Commit af3bf92d authored by Jérémie Dudouet's avatar Jérémie Dudouet
Browse files

Merge branch 'TrackBuilder' into preprod

Resolved conflicts related to implementation of CMake Modern
parents b9738fbe 776c763f
......@@ -6,5 +6,3 @@ if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/build/T${AGAPRO_BUILDING_TYPE}/CMakeLists
else()
message(AUTHOR_WARNING "Building type ${AGAPRO_BUILDING_TYPE} unknown in ${CMAKE_CURRENT_SOURCE_DIR}/CMakeLists.txt")
endif()
......@@ -45,7 +45,9 @@ endif()
if (ROOT_FOUND)
SET(HAS_ROOT ON)
SET(WITH_ROOT ON) # For Vamos lib
if (VAMOS_FOUND)
SET(WITH_ROOT ON) # For Vamos lib
endif()
message("[AGAPRO] + Configured with ROOT enabled")
else()
message("[AGAPRO] + Configured with ROOT disabled")
......@@ -53,7 +55,9 @@ endif()
if (MFM_FOUND)
SET(HAS_MFM ON)
SET(WITH_MFM ON) # For Vamos lib
if (VAMOS_FOUND)
SET(WITH_MFM ON) # For Vamos lib
endif()
message("[AGAPRO] + Configured with MFM enabled")
else()
message("[AGAPRO] + Configured with MFM disabled")
......
......@@ -44,7 +44,9 @@ endif()
# ROOT related
if (ROOT_FOUND)
set(HAS_ROOT ON)
set(WITH_ROOT ON)
if (VAMOS_FOUND)
set(WITH_ROOT ON) # For Vamos lib
endif()
message(STATUS " --> with ROOT enabled")
else()
message(WARNING " --> with ROOT disabled")
......@@ -52,7 +54,9 @@ endif()
# MFM related
if (MFM_FOUND)
set(HAS_MFM ON)
set(WITH_MFM ON)
if (VAMOS_FOUND)
set(WITH_MFM ON) # For Vamos lib
endif()
message(STATUS " --> with MFM enabled")
else()
message(AUTHOR_WARNING " --> with MFM disabled")
......
......@@ -89,6 +89,8 @@ void TB_AGATA_Builder::InitTree()
fTree->Branch("nbHitsperCry", nbHitsperCry, "nbHitsperCry[nbCores]/I");
fTree->Branch("coreId", coreId, "coreId[nbCores]/I");
fTree->Branch("coreE0", coreE0, "coreE0[nbCores]/F");
fTree->Branch("coreE1", coreE1, "coreE1[nbCores]/F");
fTree->Branch("coreT0", coreT0, "coreT0[nbCores]/F");
fTree->Branch("coreTS", coreTS, "coreTS[nbCores]/l");
fTree->Branch("nbAdd", &nbAdd, "nbAdd/I");
......@@ -212,7 +214,9 @@ void TB_AGATA_Builder::Process(Int_t idet)
// Core like branches
coreId[i] = data->GetUID();
coreE0[i] = data->GetE();
coreE0[i] = data->GetE(0);
coreE1[i] = data->GetE(1);
coreT0[i] = data->GetT(0);
coreTS[i] = ((AgataKey *)fFrame->GetFrame()->GetKey())->GetTimeStamp();
// Loop on the hit to get over the hits with the max of the energy
......@@ -473,6 +477,11 @@ void TB_AGATA_Tracking::InitTree()
fTree->Branch("trackX1", trackX1, "trackX1[nbTrack]/F");
fTree->Branch("trackY1", trackY1, "trackY1[nbTrack]/F");
fTree->Branch("trackZ1", trackZ1, "trackZ1[nbTrack]/F");
fTree->Branch("trackX2", trackX2, "trackX2[nbTrack]/F");
fTree->Branch("trackY2", trackY2, "trackY2[nbTrack]/F");
fTree->Branch("trackZ2", trackZ2, "trackZ2[nbTrack]/F");
fTree->Branch("trackType", trackType, "trackType[nbTrack]/I");
fTree->Branch("trackFOM", trackFOM, "trackFOM[nbTrack]/F");
fTree->Branch("trackT", trackT, "trackT[nbTrack]/F");
fTree->Branch("trackCrystalID", trackCrystalID, "trackCrystalID[nbTrack]/I");
fTree->Branch("TStrack", &TSTrack, "TStrack/l");
......@@ -521,6 +530,9 @@ void TB_AGATA_Tracking::Process(Int_t idet)
trackE[i] = gamma1->GetE();
trackType[i] = gamma1->GetType();
trackFOM[i] = gamma1->GetFOM();
trackT[i] = hit1->GetT();
trackCrystalID[i] = hit1->GetID(1);
......@@ -528,6 +540,19 @@ void TB_AGATA_Tracking::Process(Int_t idet)
trackY1[i] = gamma1->GetY();
trackZ1[i] = gamma1->GetZ();
//if tracked gamma is photo or pair, the second int point is a copy of the first one (because not useful)
if(trackType[i]==2){
trackX2[i] = gamma1->GetX();
trackY2[i] = gamma1->GetY();
trackZ2[i] = gamma1->GetZ();
}
else{
trackX2[i] = trackX1[i];
trackY2[i] = trackY1[i];
trackZ2[i] = trackZ1[i];
}
if(trackE[i]==0)
nbTrack--;
}
......
......@@ -64,11 +64,13 @@ protected:
Int_t nbCores; // for more than 1 frame (built event)
Int_t nbHitsperCry[MaxCore]; // nb of hits epr crystal
Int_t coreId[MaxCore]; // for each crystal fired its ID
Float_t coreE0[MaxCore]; // Core energie, high gain
Float_t coreE0[MaxCore]; // Core energy, high gain
Float_t coreE1[MaxCore]; // Core energy, low gain
Float_t coreT0[MaxCore]; // Core time, high gain
Int_t nbAdd; // Addback multiplicity
Int_t AddId[MaxCore]; // Addback crystal Id of the hit with highest energy
Float_t AddE[MaxCore]; // Addback Energy
Float_t AddE[MaxCore]; // Addback energy
Float_t AddX[MaxCore]; // XYZ Position of the hit with highest energy
Float_t AddY[MaxCore];
Float_t AddZ[MaxCore];
......@@ -113,6 +115,12 @@ protected:
Float_t trackX1[MaxTrackedGamma]; // X position of the first interaction point
Float_t trackY1[MaxTrackedGamma]; // Y position of the first interaction point
Float_t trackZ1[MaxTrackedGamma]; // Z position of the first interaction point
Float_t trackX2[MaxTrackedGamma]; // X position of the Second interaction point
Float_t trackY2[MaxTrackedGamma]; // Y position of the Second interaction point
Float_t trackZ2[MaxTrackedGamma]; // Z position of the Second interaction point
Int_t trackType[MaxTrackedGamma]; // Interaction type (1-photo 2-compt 3-pair)
Float_t trackFOM[MaxTrackedGamma]; // figure of merit
Float_t trackT[MaxTrackedGamma]; // time of the track gamma
Int_t trackCrystalID[MaxTrackedGamma]; // CrystalID of the first interaction point of the track
......
......@@ -240,7 +240,7 @@ UInt_t TreeBuilder::ProcessBlock(ADF::FrameBlock &inBlock)
fReadEvts++;
if(fReadEvts%fAutoSave == 0) {
fTree->AutoSave("SaveSelf");
fTree->AutoSave();
}
}
......
......@@ -58,7 +58,7 @@ protected:
protected:
static const Long64_t fMaxFileSize = 2*1024*1024*1024LL; //! the output file number is incremented all fMaxFileSize (Here 2GB)
static const Long64_t fAutoSave = 1*100*1000; //! the tree is saved all fAutoSave events (Here, 1e5 evts)
static const Long64_t fAutoSave = 1*1000*1000; //! the tree is saved all fAutoSave events (Here, 1e6 evts)
string fSaveDir = "";
string fNamePatern = "";
......
......@@ -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=-1.;
//fForceSGtoCC = false;
fWriteMgtData = false;
......@@ -714,6 +717,16 @@ UInt_t TrackingFilter::ProcessBlock(ADF::FrameBlock &inBlock, ADF::FrameBlock &o
fTrackEvents[0]++;
fTrackEvents[1] += number_of_hits;
}
if(fVerbose) {
cout<<endl<<endl;
cout<<"List of hits raw:"<<endl;
for(ushort t=0 ; t<number_of_hits ; t++){
cout<<" DetId: "<<pEXYZ[t].Ind<<" - hitE: "<<pEXYZ[t].E<<" - hitT: "<<pEXYZ[t].T<<endl;
}
cout<<endl;
}
if(fExcludeTracking) {
errStatus = 0; // OK if tracking excluded
}
......@@ -721,7 +734,94 @@ 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));
uint index=0; // position in the initial array of hits
int tot_number_of_gammas=0; // cumulative number of tracked gamma rays
first_hit = 0; // position of the first hit of a block, to be set in the active Tracking class
int NBlocks=0; // Number of blocks of hits
// loop on all the hits
while(index<saved_number_of_hits) {
if(fVerbose) cout<<"New block (TrackingTimeWindow="<<fTrackTimeWindow<<"):"<<endl;
NBlocks++;
// reference time of the current block
Float_t RefTime = pEXYZ[index].T;
// number of hits in the current block
number_of_hits=0;
// loop on all the hits within the tracking time window
while(pEXYZ[index].T-RefTime<=fTrackTimeWindow && index<saved_number_of_hits){
pEXYZ[number_of_hits] = pEXYZ[index];
number_of_hits++;
index++;
if(fVerbose) {
cout<<"DetId:"<<setw(2)<<pEXYZ[number_of_hits].Det<<" (seg "<<setw(2)<<pEXYZ[number_of_hits].Seg<<")";
cout<<" - hiT:"<<setw(7)<<pEXYZ[number_of_hits].T;
cout<<" - dt in block="<<setw(8)<<pEXYZ[number_of_hits].T-RefTime;
cout<<" - dt Global="<<setw(8)<<pEXYZ[number_of_hits].T-pEXYZ[0].T;
cout<<" - hitE="<<setw(8)<<pEXYZ[number_of_hits].E<<endl;
}
}
if(fVerbose) 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;
if(fVerbose) {
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));
if(fVerbose) {
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;
......@@ -856,7 +956,7 @@ Int_t TrackingFilter::SetInput()
return 18; // ??
}
else {
trigType[nTriggerEventData_10]++; // event:data:psa WITHOUT data:rancN
trigType[nTriggerEventData_10]++; // eve]nt:data:psa WITHOUT data:rancN
}
}
......@@ -945,6 +1045,13 @@ Int_t TrackingFilter::SetInput()
for(UShort_t nh = 0; nh < number_of_hits; nh++, pt++) {
pt->T += cryst[pt->Ind].tstamp - timestampEvt;
}
// ordering the hist in time rather than in detector id
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]);
}
}
if(hasAncillary)
timestampAnc += timestampAnc-timestampEvt;
......@@ -1181,6 +1288,8 @@ Int_t TrackingFilter::SetOutput()
number_of_tracked_gammas += number_of_gammas;
if(fVerbose) cout<<"Number of tracked gammas: "<<number_of_gammas<<endl;
// loop on the tracked gammas
trGamma *pG = pGams;
int numUsed = 0;
......@@ -1192,20 +1301,31 @@ Int_t TrackingFilter::SetOutput()
ahit->SetT(pG->T); // time of the tracked gamma
ahit->SetType((UShort_t)pG->type);
ahit->SetFOM((Double_t)pG->FOM);
Float_t ESum=0.;
if(fVerbose) cout<<"gamma "<<ii<<": "<<pG->T<<" E: "<<pG->E<<endl;
// the interaction points of the gamma; still writing only first 2 for compton and only one for pair-production
for (Int_t nh = 0; nh < pG->nH; nh++) { // Decision on how many hits to write taken later in (pG-type!=2) and if(nh==1)
exyzHit *pP = pEXYZ + pG->iH[nh]; // Original hit of this interaction point; all info taken from here
Hit *pH = ahit->NewHit();
pH->SetE(pP->E);
pH->SetXYZ(pP->X, pP->Y, pP->Z);
pH->SetT(pP->T);
pH->SetT(pP->T-pEXYZ[0].T); // Hit T related to the first hit time value
pH->SetID(pP->Seg, 0);
pH->SetID(pP->Det, 1);
ESum+=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
if(fVerbose) cout<<"hit "<<nh<<": "<<pP->T<<" E:"<<pP->E<<endl;
}
if(fVerbose) cout<<" ==> ESum= "<<ESum<<endl;
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 +1633,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 +1921,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;
......@@ -2004,7 +2130,7 @@ void TrackingFilter::WriteInputHits()
fprintf(fp_hits, "\n");
}
fprintf(fp_hits, "-1 %.2f 0 0 0 %d\n", etot, evnumberCry); // evnumberCry ?? seems meaningless; should become timestampEvt
fprintf(fp_hits, "-1 %.2f 0 0 0 %llu\n", etot, timestampEvt);
if(trueHits<=0)
return;
......
......@@ -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