Commit 4df97edf authored by Jérémie Dudouet's avatar Jérémie Dudouet Committed by Jérémie Dudouet
Browse files

Add a pre tracking time window to build blocks of hits to be sent to the tracking

parent ec86ba3e
......@@ -89,6 +89,7 @@ TrackingFilter::TrackingFilter() :
rotTr = NULL;
pGams = NULL;
number_of_hits = 0;
first_hit = 0;
number_of_gammas = 0;
number_of_tracked_gammas = 0;
......@@ -122,6 +123,8 @@ TrackingFilter::TrackingFilter() :
fTsGeAnc_min = 0x80000000;
fTsGeAnc_max = 0x7FFFFFFF;
fTrackTimeWindow=0.;
//fForceSGtoCC = false;
fWriteMgtData = false;
......@@ -721,7 +724,86 @@ UInt_t TrackingFilter::ProcessBlock(ADF::FrameBlock &inBlock, ADF::FrameBlock &o
errStatus = TrackingFilter::Process(); // accept everything as individual gammas
}
else {
errStatus = Process(); // normal daughter::Process()
// if the tracking hits window has been define, blocks of hits in this time window are sent to the tracking
if(fTrackTimeWindow>=0){
// saving the initial block of hits
int saved_number_of_hits = number_of_hits;
memcpy(saved_pEXYZ,pEXYZ,saved_number_of_hits*sizeof(exyzHit));
// sorting of hits in time
for(ushort t=1 ; t<number_of_hits ; t++)
for(int i=0 ; i<number_of_hits-1 ; i++)
if(pEXYZ[i].T > pEXYZ[i+1].T){
swap(pEXYZ[i],pEXYZ[i+1]);
}
// reference time for the next builded block
uint index=0;
int tot_number_of_gammas=0;
first_hit = 0;
while(index<saved_number_of_hits) {
// reference time for the next builded block
Float_t RefTime = saved_pEXYZ[index].T;
// value sent to the tracking
number_of_hits=0;
// cout<<"New block (TrackingTimeWindow="<<fTrackTimeWindow<<"):"<<endl;
while(saved_pEXYZ[index].T-RefTime<=fTrackTimeWindow && index<saved_number_of_hits){
pEXYZ[number_of_hits] = saved_pEXYZ[index];
// cout<<setw(2)<<pEXYZ[number_of_hits].Det<<" (seg "<<setw(2)<<pEXYZ[number_of_hits].Seg<<") dt="<<setw(8)<<pEXYZ[number_of_hits].T-RefT<<" ; T="<<setw(8)<<pEXYZ[number_of_hits].T<<" E="<<setw(8)<<pEXYZ[number_of_hits].E<<endl;
number_of_hits++;
index++;
}
// cout<<"nhits in block: "<<number_of_hits<<endl;
//Process tracking for the selected block of gammas
Process();
// shift first_hit by the number of already processed hits
first_hit += number_of_hits;
// Return to first position in pGams
pGams = &pGams[0];
// add the new tracked gamma rays to the list saved_pGams
memcpy(&saved_pGams[tot_number_of_gammas],pGams,number_of_gammas*sizeof(trGamma));
tot_number_of_gammas += number_of_gammas;
// cout<<number_of_gammas<<" tracked gammas"<<endl;
// for(int ii=0 ; ii<number_of_gammas ; ii++){
// cout<<pGams[ii].E<<" "<<pGams[ii].nH<<"(";
// for(int iii=0 ; iii<pGams[ii].nH ; iii++) cout<<pGams[ii].iH[iii]<<" ";
// cout<<")"<<endl;
// }
}
// recover the real total number of hits
number_of_hits = saved_number_of_hits;
memcpy(pEXYZ,saved_pEXYZ,number_of_hits*sizeof(exyzHit));
// recover the real total number of gammas
number_of_gammas = tot_number_of_gammas;
memcpy(pGams,saved_pGams,number_of_gammas*sizeof(trGamma));
// cout<<"Final traked result:"<<endl;
// cout<<number_of_gammas<<" tracked gammas"<<endl;
// for(int ii=0 ; ii<number_of_gammas ; ii++){
// cout<<pGams[ii].E<<" "<<pGams[ii].nH<<"(";
// for(int iii=0 ; iii<pGams[ii].nH ; iii++) {
// exyzHit *pP = &pEXYZ[pGams[ii].iH[iii]];
// cout<<pGams[ii].iH[iii]<<":"<<pP->E;
// }
// cout<<")"<<endl;
// }
}
else
errStatus = Process(); // normal daughter::Process()
}
if( errStatus ) {
error_code = 2;
......@@ -1181,11 +1263,15 @@ Int_t TrackingFilter::SetOutput()
number_of_tracked_gammas += number_of_gammas;
Float_t sum1=0.;
Float_t sum2=0.;
// loop on the tracked gammas
trGamma *pG = pGams;
int numUsed = 0;
for (Int_t ii = 0; ii < number_of_gammas; ii++, pG++) {
TrackedHit *ahit = trackdata->NewGamma();
sum1+=pG->E;
ahit->SetStatus(0);
ahit->SetE(pG->E); // energy, passed over in units of keV
ahit->SetXYZ(pG->X1, pG->Y1, pG->Z1); // 3D position of the tracked gamma is the first hit
......@@ -1201,11 +1287,13 @@ Int_t TrackingFilter::SetOutput()
pH->SetT(pP->T);
pH->SetID(pP->Seg, 0);
pH->SetID(pP->Det, 1);
sum2+=pP->E;
//if(pG->type != 2) break; // pair-production writes only one point (like photo, of course)
//if(nh == 1) break; // limit compton scattering to firse two points
}
numUsed += pG->nH;
}
if (!fDiscardEmpty && numUsed < number_of_hits) {
// write the unused hits in a gamma with E=0
TrackedHit *ahit = trackdata->NewGamma();
......@@ -1513,6 +1601,9 @@ UInt_t TrackingFilter::GetParameters(const std::string& confFile, Bool_t doList)
conf.Add("TimeWindowGeAnc", "time gate for ge-ancillary"
, &fTsGeAnc_min, &fTsGeAnc_max);
conf.Add("TrackingTimeWindow", "only successive packs hits in this time window will be sent to the tracking (useful for isomer)"
, &fTrackTimeWindow);
conf.Add("RotoTranslations", "file to transform coordinates from intrinsic to space"
, &fWhichRotoTranslations);
conf.Add("GeometrySummary", "for the header of the Mgt_Hits.txt file of input hits"
......@@ -1798,7 +1889,10 @@ void TrackingFilter::InitGenStructures()
for(int nn = 0; nn < 180; nn++) rotTr[nn].rot[0][0] = rotTr[nn].rot[1][1] = rotTr[nn].rot[2][2] = 1;
pEXYZ = (exyzHit *)calloc(intmax, sizeof(exyzHit));
saved_pEXYZ = (exyzHit *)calloc(intmax, sizeof(exyzHit));
pGams = (trGamma *)calloc(intmax, sizeof(trGamma));
saved_pGams = (trGamma *)calloc(intmax, sizeof(trGamma));
number_of_hits = 0;
number_of_gammas = 0;
......
......@@ -101,7 +101,7 @@ protected:
float Ecorr;
float X1, Y1, Z1; // same as original PSA hit iH[0] (should be removed)
float T; // timing of gamma
int type; // 0-photo 1-compt 2-pair
int type; // 1-photo 2-compt 3-pair
double FOM; // figure of merit of cluster
int nH; // number of hits of this gamma
int iH[64]; // there should be a global parameter for the maximum number of hits in the gammas (oft::kmax=7 mgt::MAXACCEPT=6)
......@@ -131,11 +131,14 @@ protected:
bool KeepCrystal[180]; // to manage time order of crystals
exyzHit *pEXYZ;
exyzHit *saved_pEXYZ; // used for tracking time gate
int number_of_hits;
int first_hit; // first hit of a time window pack
rotTr3D *rotTr; // placement of detectors
trGamma *pGams;
trGamma *saved_pGams;
int number_of_gammas;
int number_of_tracked_gammas; // total number of tracked gammas
......@@ -204,7 +207,7 @@ protected:
UInt_t evnumberEvt; // this should have the same value for all data frames in the event
Int_t fTsGeAnc_min, fTsGeAnc_max; // the time gate on Tstamp(ge)-Tstamp(ancillary)
Float_t fTrackTimeWindow;
Bool_t stop_called;
bool hasAncillary; // event has ancillary component
......
......@@ -104,7 +104,7 @@ void TrackingFilterMGT::FillMGTStructures()
pc->ns = pt->Seg; // segnum
pc->nd = pt->Det; // detnum
pc->ng = -1; // not assigned to a gamma
pc->pid = nn; // id of hit
pc->pid = first_hit+nn; // id of hit
eseen += pt->E;
if(fVerbose) {
if(nn==0)
......
......@@ -107,12 +107,12 @@ void TrackingFilterOFT::FillOFTStructures()
local_ptr->source.z = 0.1*fSource[2];
local_ptr->nb_int = number_of_hits;
exyzHit * pt = pEXYZ;
for(int nn = 0; nn < number_of_hits; nn++, pt++) {
local_ptr->eh[nn] = pt->E * 0.001; // MeV
for(int nn = 0 ; nn < number_of_hits; nn++, pt++) {
local_ptr->eh[nn] = pt->E * 0.001; // MeV
local_ptr->ipos[nn].x = pt->X * 0.1; // cm
local_ptr->ipos[nn].y = pt->Y * 0.1; // cm
local_ptr->ipos[nn].z = pt->Z * 0.1; // cm
local_ptr->ipos[nn].iH = pt->Ind = nn; // which hit is this;
local_ptr->ipos[nn].iH = pt->Ind = first_hit+nn; // which hit is this;
local_ptr->DetId[nn] = pt->Det;
if(fVerbose) {
if(nn==0)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment