Modified interface to oft so as to get back the full information on the hits of the tracked gammas

The mean free paths in mgt are now interpolated between the points listed in the mgt_ge_[photo|comp|pair].dat

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1139 170316e4-aea8-4b27-aad4-0380ec0519c9

WinCtest/femul.vcproj

@@ -46,7 +46,7 @@
 				Optimization="0"
 				EnableIntrinsicFunctions="true"
 				AdditionalIncludeDirectories="..\WinCtest;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Preprocessing\includePrePSA;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\filters\PostPSA;..\builders;..\filters\Global;..\..\include;..\PRISMA\src\lib_prisma\include;"C:\Program Files (x86)\boost\boost_1_40";C:\root\include"
-				PreprocessorDefinitions="WIN32;_DEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;NRV_TYPE=NRV_OFFLINE;TF_ROOTTREE"
+				PreprocessorDefinitions="WIN32;_DEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;NRV_TYPE=NRV_OFFLINE;TRF_ROOTTREE"
 				MinimalRebuild="true"
 				BasicRuntimeChecks="3"
 				RuntimeLibrary="3"
@@ -202,8 +202,8 @@
 				Name="VCCLCompilerTool"
 				Optimization="2"
 				EnableIntrinsicFunctions="true"
-				AdditionalIncludeDirectories="..\PRISMA\src\lib_prisma\include;"C:\Program Files (x86)\boost\boost_1_40";..\..\agaprodep\adf;..\WinCtest;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Preprocessing\includePrePSA;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\builders;C:\root\include;..\filters\PostPSA;..\filters\Global"
-				PreprocessorDefinitions="WIN32;NDEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;NRV_TYPE=NRV_OFFLINE;TF_ROOTTREE"
+				AdditionalIncludeDirectories="..\PRISMA\src\lib_prisma\include;"C:\Program Files (x86)\boost\boost_1_40";..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Preprocessing\includePrePSA;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\builders;C:\root\include;..\filters\PostPSA;..\filters\Global;..\..\include"
+				PreprocessorDefinitions="WIN32;NDEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;NRV_TYPE=NRV_OFFLINE;TRF_ROOTTREE"
 				RuntimeLibrary="2"
 				EnableFunctionLevelLinking="true"
 				UsePrecompiledHeader="0"

common/commonDefs.h

@@ -179,7 +179,7 @@ const int     defTriggerSample  =  10;  // 60-10=50 samples passed to the PSA
 ////// TrackingFilter //////
 ////////////////////////////
 
-#define TRF_ROOTTREE          // enable the root tree in this actor (ENABLED IN THE MAKEFILE)
+//#define TRF_ROOTTREE          // enable the root tree in this actor (ENABLED IN THE MAKEFILE)
 //#define TRF_WRITE_TRACKED     // enable writing the tracked gammas into the Oft_TrackedEnergies ascii file
 //#define TRF_WRITE_GSORT       // enable writing input and tracked data in gasp format
 #define TFR_WRITE_HITS_MGT    // enable writing the input hits in mgt format

filters/Tracking/TrackingFilter.cpp

@@ -923,7 +923,7 @@ Int_t TrackingFilter::SetInput()
     evnumberEvt  = evnumberCry;
   }
 
-  // correct the time information of the hits to refer to timestampEvt
+  // correct the time information of the hits to refer to timestampEvt (?????)
   exyzHit *pt = pEXYZ;
   for(UShort_t nh = 0; nh < number_of_hits; nh++, pt++) {
     pt->T += cryst[pt->Ind].tstamp - timestampEvt;
@@ -1048,22 +1048,20 @@ void TrackingFilter::GetDataPSA()
 
   // loop over number of interaction points in fFramePSA
   float esumhits = 0;
-  float averageT = 0;
   int   hit1st   = number_of_hits;
   exyzHit *pt = pEXYZ + number_of_hits;
   for(UShort_t nh = 0; nh < nbhits; nh++, pt++) {
     const PSAHit *ahit = dynamic_cast<PSAHit*>(data->GetHit(nh));
-    //pt->Det = crystal_id;         // needed in previous versions
-    pt->Det = ahit->GetID(1);
-    pt->Seg = ahit->GetID(0);
-    pt->Ind = number_of_crystals;   // used to recover crystal-level infos (tstamp, ...)
-    pt->E  = (float)ahit->GetE(); // passed in units of keV
-    pt->X = (float)ahit->GetX();
+    pt->E = (float)ahit->GetE();    // keV
+    pt->X = (float)ahit->GetX();    // mm
     pt->Y = (float)ahit->GetY();
     pt->Z = (float)ahit->GetZ();
     //pt->T = (float)ahit->GetT();  // not clear what this is --> replaced with common time
     pt->T = CoreT[0];               // same time for all hits; this should be in units of samples
-    averageT += pt->T* pt->E;       // nonsense as long as all times are the same
+    pt->Det = ahit->GetID(1);
+    pt->Seg = ahit->GetID(0);
+    pt->Ind = number_of_crystals;   // used to recover crystal-level infos (tstamp, ...)
+    pt->Gam = -1;                   // not yet assiged to a gamma
     esumhits += pt->E;              // build the sum-of-segments energy
 
 #ifdef TRF_FromGRU      
@@ -1110,10 +1108,6 @@ void TrackingFilter::GetDataPSA()
   //}
 
   number_of_hits += nbhits;
-  if(esumhits > 0)
-    averageT /= esumhits;
-  else
-    averageT  = 0;
 
 #ifdef TRF_FromGRU
   CE_Tracking->Fill(esumhits*10);
@@ -1227,22 +1221,18 @@ Int_t TrackingFilter::SetOutput()
     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
-    exyzHit *pP = pEXYZ + pG->I1;           // original PSA hit of the the first interaction point of the tracked gamma
-    for(Int_t nh = 0; nh < pG->nH; nh++) {  // Decision on how many hits to write taken later in  if(nh==0) and if(pG-type!=2)
-      Hit *pH = ahit->NewHit();             // BEST OF ALL WOULD BE TO WRITE ALL HITS IN THE PROPER SEQUENCE!! (OR TO USE A FLAG)
-      pH->SetXYZ(pP->X, pP->Y, pP->Z);      // info taken from the original PSA hit
+    ahit->SetT(pG->T);                      // time of the tracked gamma
+    // 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->SetT (pP->T);                     // info taken from the original PSA hit
+      pH->SetXYZ(pP->X, pP->Y, pP->Z);
+      pH->SetT (pP->T);
       pH->SetID(pP->Seg,0);
       pH->SetID(pP->Det,1);
-      if(nh==0) {
-        ahit->SetT(pP->T);                  // time of the tracked gamma taken from its first hit
-        //ahit->SetID(pP->Seg, 0);    // there is no such field in TrackedHit
-        //ahit->SetID(pP->Ind, 1);    //  ""
-        if(pG->type!=2)                     // Photoelectric and Pair Production write only on hit. When reading back the
-          break;                            // adf file they differ in Egamma_photo = Ehit vs Egamma_PP ~ Ehit + 1022
-        pP = pEXYZ + pG->I2;
-      }
+      if(pG->type != 2) break;            // pair-production writes only on point (like photo, of course) 
+      if(nh == 1) break;                  // limit compton scattering to firse two points
     }
   }
 
@@ -1683,7 +1673,7 @@ void TrackingFilter::InitGenStructures()
     TrackSpec_ES->setComment("sum-energy spectra (gain 1 keV/chan)");
     hGroup.add(TrackSpec_ES);
 
-    TrackSpec_EE = new MultiHist<unsigned int>(2, 10, specLenE);
+    TrackSpec_EE = new MultiHist<unsigned int>(2, 15, specLenE);
     TrackSpec_EE->setFileName(fOdirPrefix+"Track?EE.spec");
     TrackSpec_EE->setComment("input and tracked spectra");
     hGroup.add(TrackSpec_EE);
@@ -1923,14 +1913,8 @@ bool TrackingFilter::OpenInputHits()
   // 6  time
   // 7  interaction    (rightmost character)  (for us this is timeStamp+eventNumber)
   fprintf(fp_hits, "OUTPUT_MASK 111001");
-  if(fWriteMgtTiming)
-    fprintf(fp_hits, "1");
-  else
-    fprintf(fp_hits, "0");
-  if(fWriteMgtExtended)
-    fprintf(fp_hits, "1");
-  else
-    fprintf(fp_hits, "0");
+  fprintf(fp_hits, fWriteMgtTiming   ? "1" : "0" );
+  fprintf(fp_hits, fWriteMgtExtended ? "1" : "0" );
   fprintf(fp_hits, "\n");
 
   fprintf(fp_hits, "DISTFACTOR 1\n");
@@ -2046,6 +2030,43 @@ int TrackingFilter::PostProcessEvent()
       TrackSpec_Anc->Incr(0, nn, (int)rawBuf[nn]);
     }
   }
+  // input and unused hits
+  float eDetUnused[180]; memset(eDetUnused, 0, sizeof(eDetUnused));
+  int   dUnusedFirst = 179, dUnusedLast = 0;
+  float eSumUnused  = 0;
+  int   countUnused = 0;
+  exyzHit * pt = pEXYZ;
+  for(int nn = 0; nn < number_of_hits; nn++, pt++) {
+    float ehit= pt->E;
+    TrackSpec_EE->Incr  (0, 10, int(fEnergyGain*ehit));                     // [0][A] input hits
+    if(hasRecoil)
+      TrackSpec_EE->Incr(1, 10, int(fEnergyGain*ehit));                     // [1][A] input hits
+    if(pt->Gam <0) {
+      countUnused++;
+      TrackSpec_EE->Incr  (0, 11, int(fEnergyGain*ehit));                   // [0][B] unused hits
+      if(hasRecoil)
+        TrackSpec_EE->Incr(1, 11, int(fEnergyGain*ehit));                   // [1][B] unused hits
+      eSumUnused += ehit;
+      int detId = pt->Ind;
+      eDetUnused[detId] += ehit;
+      dUnusedFirst = min(dUnusedFirst, detId);
+      dUnusedLast  = min(dUnusedLast,  detId);
+    }
+  }
+  // sumspec of unused hits
+  if(countUnused) {
+    TrackSpec_EE->Incr  (0, 12, int(fEnergyGain*eSumUnused));               // [0][C] unused hits sum spectrum
+    if(hasRecoil)
+      TrackSpec_EE->Incr(1, 12, int(fEnergyGain*eSumUnused));               // [1][C] unused hits sum spectrum
+    for(int nn = dUnusedFirst; nn <= dUnusedLast; nn++) {
+      float eDet = eDetUnused[nn];
+      if(eDet) {
+        TrackSpec_EE->Incr  (0, 13, int(fEnergyGain*eDet));                 // [0][D] unused common spectrum of crystals
+        if(hasRecoil)
+          TrackSpec_EE->Incr(1, 13, int(fEnergyGain*eDet));                 // [1][D] unused common spectrum of crystals
+      }
+    }
+  }
 #endif
 
   if(!number_of_gammas)
@@ -2073,12 +2094,12 @@ int TrackingFilter::PostProcessEvent()
 #endif
 #ifdef TRF_MULTIHIST
     if(TrackSpec_EE) {
-      TrackSpec_EE->Incr(0, 2,                   int(fEnergyGain*pg->E));   // 0-2 Tracked
-      TrackSpec_EE->Incr(0, 2 + pg->type,        int(fEnergyGain*pg->E));   // 0-3 Photo=1   0-4 Compt=2 0-5 pair=3
+      TrackSpec_EE->Incr(0, 2,                   int(fEnergyGain*pg->E));   // [0][2] Tracked
+      TrackSpec_EE->Incr(0, 2 + pg->type,        int(fEnergyGain*pg->E));   // [0][3] Photo  [0][4] Compt   [0][5] Pair
       eSumTracked += pg->E;
       if(hasRecoil) {
-        TrackSpec_EE->Incr(1, 2,                 int(fEnergyGain*pg->E));   // 1-2 Tracked in coinc with ancillary
-        TrackSpec_EE->Incr(1, 2 + pg->type,      int(fEnergyGain*pg->E));   // 1-3 Photo=1   1-4 Compt=2 1-5 pair=3
+        TrackSpec_EE->Incr(1, 2,                 int(fEnergyGain*pg->E));   // [1][2] Tracked in coinc with ancillary
+        TrackSpec_EE->Incr(1, 2 + pg->type,      int(fEnergyGain*pg->E));   // [1][3] Photo  [1][4] Compt   [1][5] Pair
         eSumTrackedA += pg->E;
       }
     }
@@ -2096,9 +2117,9 @@ int TrackingFilter::PostProcessEvent()
     if(Ecorr > 0) {
 #ifdef TRF_MULTIHIST
       if(TrackSpec_EE) {
-        TrackSpec_EE->Incr  (0, 6, int(fEnergyGain*Ecorr)); // 0-6 Doppler corrected
+        TrackSpec_EE->Incr  (0, 6, int(fEnergyGain*Ecorr));                 // [0][6] Doppler corrected
         if(hasRecoil)
-          TrackSpec_EE->Incr(1, 6, int(fEnergyGain*Ecorr)); // 1-6 Doppler corrected in coincidence with ancillary
+          TrackSpec_EE->Incr(1, 6, int(fEnergyGain*Ecorr));                 // [1][6] Doppler corrected in coincidence with ancillary
       }
 #endif  //TRF_MULTIHIST
 #ifdef TRF_FromGRU
@@ -2108,12 +2129,12 @@ int TrackingFilter::PostProcessEvent()
 
   }  // loop over the gammas
 
+#ifdef TRF_MULTIHIST
   // sum spectra
-  TrackSpec_EE->Incr  (0, 7, int(fEnergyGain*eSumTracked));   // 0-7 sum energy
+  TrackSpec_EE->Incr  (0, 7, int(fEnergyGain*eSumTracked));                 // [0][7] tracked sum energy
   if(hasRecoil)
-    TrackSpec_EE->Incr(1, 7, int(fEnergyGain*eSumTrackedA));  // 1-7 sum energy
+    TrackSpec_EE->Incr(1, 7, int(fEnergyGain*eSumTrackedA));                // [1][7] tracked sum energy
 
-#ifdef TRF_MULTIHIST
   // time of the individual gammas and time between gammas using only info from the output
   if(TrackSpec_TG) {
     for(int i1 = 0; i1 < number_of_gammas-1; i1++) {
@@ -2121,7 +2142,7 @@ int TrackingFilter::PostProcessEvent()
       TrackSpec_TG->Incr(0, int(gamTime1));   // time of the individual gammas (re global time reference ??) 
       for(int i2 = i1+1; i2 < number_of_gammas; i2++) {
         Double_t gamTime2 = pGams[i2].T;
-        TrackSpec_TG->Incr(1, int(gamTime2-gamTime1 + 2*specLenT1));  // time between gammas
+        TrackSpec_TG->Incr(1, int(gamTime2-gamTime1 + 2*specLenT1));        // time between gammas
       }
     }
   }
@@ -2160,15 +2181,15 @@ bool TrackingFilter::AnalysisOfCrystals()
     int id1 = SM180(pCr1->id);           // off range should return fNumgeDets-1
     esumCC += pCr1->CCE0;
     if(TrackSpec_EC) {
-      TrackSpec_EC->Incr(0, id1, int(fEnergyGain*pCr1->CCE0));     // energy of the core sorted by detector
-      TrackSpec_EC->Incr(1, id1, int(fEnergyGain*pCr1->Esum));     // sum energy of the segments
+      TrackSpec_EC->Incr(0, id1, int(fEnergyGain*pCr1->CCE0));              // energy of the core sorted by detector
+      TrackSpec_EC->Incr(1, id1, int(fEnergyGain*pCr1->Esum));              // sum energy of the segments
     }
     if(TrackSpec_EE) {
-      TrackSpec_EE->Incr(0,   0, int(fEnergyGain*pCr1->CCE0));     // 0-0 energy of the core common spectrum
-      TrackSpec_EE->Incr(0,   1, int(fEnergyGain*pCr1->Esum));     // 0-1 energy of the sumsegs common spectrum
+      TrackSpec_EE->Incr(0,   0, int(fEnergyGain*pCr1->CCE0));              // [0][0] energy of core common spectrum
+      TrackSpec_EE->Incr(0,   1, int(fEnergyGain*pCr1->Esum));              // [0][1] energy of sumsegs common spectrum
     }
     if(TrackSpec_TC)
-      TrackSpec_TC->Incr(   id1, int(tScale*pCr1->CCT0));         // trigger point of the input traces from PreprocessingFilterPSA
+      TrackSpec_TC->Incr(   id1, int(tScale*pCr1->CCT0));                   // trigger point of the input traces from PreprocessingFilterPSA
     //if(TrackMatr_ETC)
     //  TrackMatr_ETC->Incr(int(pCr1->CCE0/10), id1*1024+int(tScale*pCr1->CCT0));     // energy-trigger point from preprocessing
     int seg1 = 36*id1 + pEXYZ[pCr1->hit1st  ].Seg;
@@ -2207,9 +2228,9 @@ bool TrackingFilter::AnalysisOfCrystals()
     }
   }
   if(TrackSpec_EE) {
-    TrackSpec_EE->Incr(0,   8, int(fEnergyGain*esumCC));        // 0-8 input sum energy
+    TrackSpec_EE->Incr(0,   8, int(fEnergyGain*esumCC));                    // [0][8] input sum energy
     if(number_of_crystals==1)
-      TrackSpec_EE->Incr(0,   9, int(fEnergyGain*esumCC));      // 0-9 input single detector
+      TrackSpec_EE->Incr(0,   9, int(fEnergyGain*esumCC));                  // [0][9] input single crystal
   }
   if(TrackSpec_ES && number_of_crystals > 0 && number_of_crystals < fNumGeDets) {
     TrackSpec_ES->Incr(0,                  0, int(esumCC));     // total sum energy
@@ -2282,13 +2303,13 @@ bool TrackingFilter::AnalysisGeAncillary()
     // ge-Ancillary in the time window
     if(isInTime) {
       if(TrackSpec_EA)
-        TrackSpec_EA->Incr(1, ii1, eG);                 // coinc. energy of the core
+        TrackSpec_EA->Incr(1, ii1, eG);                                     // coinc. energy of the core
       if(TrackSpec_TA)
-        TrackSpec_TA->Incr(1, ii1, dT);                 // tstamp difference ge-Anc
+        TrackSpec_TA->Incr(1, ii1, dT);                                     // tstamp difference ge-Anc
       //if( rVc ) {
         if(TrackSpec_EE) {
-          TrackSpec_EE->Incr(1, 0, eG);                 // 1-0 total energy spectrum of the core
-          TrackSpec_EE->Incr(1, 1, eS);                 // 1-1 energy of the sumsegs common spectrum
+          TrackSpec_EE->Incr(1, 0, eG);                                     // [1][0] total energy spectrum of the core
+          TrackSpec_EE->Incr(1, 1, eS);                                     // [1][1] energy of sumsegs common spectrum
         }
       //}
       if(TrackMatr_AG) {
@@ -2447,6 +2468,9 @@ void TrackingFilter::openTTree ()
     t->Branch("hitTstamp",hitTstamp,"hitTstamp[number_of_hits]/l"); // original timestamps from the PSA
     t->Branch("hitId",hitId,"hitId[number_of_hits]/I");
     t->Branch("hitSg",hitSg,"hitSg[number_of_hits]/I");
+    // branches to be added
+    // add core as 37 (Benedikt)
+    // a mark of which gamma got the hit (-1 if not assigned)
   }
 
 }
@@ -2460,7 +2484,7 @@ void TrackingFilter::CopyHits2Tree()
     hitY[nn]  = (Float_t)pt->Y;
     hitZ[nn]  = (Float_t)pt->Z;
     hitT[nn]  = (Float_t)cryst[pt->Ind].CCT0;
-    hitTstamp[nn] = cryst[pt->Ind].tstamp;
+    hitTstamp[nn] =      cryst[pt->Ind].tstamp;
     hitId[nn] = (Int_t)pt->Det;
     hitSg[nn] = (Int_t)pt->Seg;
   }
@@ -2468,6 +2492,7 @@ void TrackingFilter::CopyHits2Tree()
 
 void TrackingFilter::CopyGamma2Tree (trGamma *pg, Int_t number_of_gammas)
 {
+  // should be generalized to write all interaction points of the gamma
   number_of_gammas = number_of_gammas > intmax ? intmax : number_of_gammas;
   for(Int_t i = 0; i < number_of_gammas; i++, pg++ ) {
     gammaE[i]  = (Float_t)pg->E;
@@ -2477,9 +2502,14 @@ void TrackingFilter::CopyGamma2Tree (trGamma *pg, Int_t number_of_gammas)
     gammaX1[i] = (Float_t)pg->X1;
     gammaY1[i] = (Float_t)pg->Y1;
     gammaZ1[i] = (Float_t)pg->Z1;
-    gammaX2[i] = (Float_t)pg->X2;
-    gammaY2[i] = (Float_t)pg->Y2;
-    gammaZ2[i] = (Float_t)pg->Z2;
+    if(pg->type == 2) {
+     gammaX2[i] = pEXYZ[pg->iH[1]].X;
+     gammaY2[i] = pEXYZ[pg->iH[1]].Y;
+     gammaZ2[i] = pEXYZ[pg->iH[1]].Z;
+    }
+    else {
+      gammaX2[i] = gammaY2[i] = gammaZ2[i] = 0;
+    }
   }
 }
 

filters/Tracking/TrackingFilter.h

@@ -42,6 +42,7 @@
   #include <TH3.h>
 #endif  // TRF_FromGRU
 
+// enable this here (to avoid full rebuild when changing them in the projet)
 //#define TRF_ROOTTREE
 //#define TRF_WRITE_GSORT
 
@@ -100,29 +101,24 @@ protected:
     float X;
     float Y;
     float Z;
-    float T;          // time of hit (samples) referred to the common reference timestampEvt
-    UInt_t Det;        // detector ID
-    UInt_t Seg;        // segment number
-    UInt_t Ind;        // sequential crystal number (from which to take all needed info)
+    float T;            // time of hit (samples) referred to the common reference timestampEvt ???
+    UInt_t Det;         // detector ID
+    UInt_t Seg;         // segment number
+    UInt_t Ind;         // sequential crystal number (from which to take all needed info)
+    Int_t  Gam;         // sequential number of the gamma to which this hit has been assigned (or -1)
   } exyzHit;
 
   // Output gammas from tracking
-  // The only information needed is actually the identifier of the original PSA hit (I1 and I2)
+  // The only information needed is actually the identifier of the original PSA hit iH[nH]
   // the rest (xyz) is redundant and could be removed
   typedef struct {
     float E;
     float Ecorr;
-    int   nH;     // number of hits of this gamma
-    int   type;   // 0-photo 1-compt 2-pair
-    int   I1;     // original PSA hit of the 1st interaction point
-    float X1;     // info recoverable from the original PSA hit I1
-    float Y1;
-    float Z1;
-    int   I2;     // original PSA hit of the 2nd interaction point
-    float X2;     // info recoverable from the original PSA hit I2
-    float Y2;
-    float Z2;
-    float T;
+    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   nH;           // number of hits of this gamma
+    int   iH[8];        // there should be a global for the maximum number of hits in the gammas (oft:kmax=7 mgt:MAXACCEPT=6)
   } trGamma;
 
   typedef struct {

filters/Tracking/TrackingFilterMGT.cpp

@@ -102,7 +102,8 @@ void TrackingFilterMGT::FillMGTStructures()
     pc->tt = pt->T*10;  // sample --> ns
     pc->ns = pt->Seg;   // segnum
     pc->nd = pt->Det;   // detnum
-    pc->pid = nn;
+    pc->ng = -1;        // not assigned to a gamma
+    pc->pid = nn;       // id of hit
     eseen += pt->E;
     if(fVerbose) {
       if(nn==0)
@@ -118,58 +119,31 @@ void TrackingFilterMGT::FillMGTStructures()
 
 void TrackingFilterMGT::MoveMGTStructures()
 {
-  char tag1[] = "gam#";
-  char tag2[] = "    ";
-  char * tag = tag1;
+  // gammas are assigned the time of their first hit
+  // it would be better to average over the interaction points, or take the time of largest energy
   number_of_gammas = local_ptr->ngoodclusters;
   trGamma * pg = pGams;
   clust   * pc = local_ptr->pgoodclusters[0];
   for(Int_t nn = 0; nn < number_of_gammas; nn++, pg++, pc++) {
-    pg->E  = pg->Ecorr = float(pc->eclust);
-    pg->nH = pc->npoints;
-    pg->X1 = float(pc->ppc[0]->xx);
-    pg->Y1 = float(pc->ppc[0]->yy);
-    pg->Z1 = float(pc->ppc[0]->zz);
-    pg->I1 = pc->ppc[0]->pid;
-    if(pg->nH > 1) {
-      pg->X2 = float(pc->ppc[1]->xx);
-      pg->Y2 = float(pc->ppc[1]->yy);
-      pg->Z2 = float(pc->ppc[1]->zz);
-      pg->I2 = pc->ppc[1]->pid;
-    }
-    else {
-      pg->X2 =  pg->Y2 = pg->Z2 = 0;
+    pg->E    = pg->Ecorr = float(pc->eclust);   // energy of gamma
+    pg->nH   = pc->npoints;                     // number of interactions points of the gamma
+    pg->type = pc->mechanism;                   // type of reconstruction (0 reserved for none)
+    int iH1  = pc->ppc[0]->pid;                 // info on first hit saved on gamma for convenience
+    pg->X1 = pEXYZ[iH1].X;                      // float(pc->ppc[0]->xx);
+    pg->Y1 = pEXYZ[iH1].Y;                      // float(pc->ppc[0]->yy);
+    pg->Z1 = pEXYZ[iH1].Z;                      // float(pc->ppc[0]->zz);
+    pg->T  = pEXYZ[iH1].T;                      // timing of gamma, taken from the first hit
+    for(Int_t kk = 0; kk < pg->nH; kk++) {
+      pg->iH[kk] = pc->ppc[kk]->pid;            // original PSA hit of interaction points
+      pEXYZ[pg->iH[kk]].Gam = pc->ppc[kk]->ng;  // gamma to which it has been assigned
     }
-    pg->type = pc->mechanism; // 0 reserved for none
-    pg->T    = pEXYZ[pg->I1].T;
-
     if(fVerbose) {
-      if(nn==1) tag = tag2;
-      printf("%s%3d        E x y z x y z = %7.1f %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f\n",
-        tag, nn, pg->E, pg->X1, pg->Y1, pg->Z1, pg->X2, pg->Y2, pg->Z2);
+      if(nn == 0) printf("gam#");
+      else        printf("    ");
+      printf("%3d %5d  E x y z x y z = %7.1f", nn, pg->nH, pg->E); 
+      for(int kk = 0; kk < min(2, pg->nH); kk++)
+        printf(" %8.1f %8.1f %8.1f  ", pEXYZ[pg->iH[kk]].X, pEXYZ[pg->iH[kk]].Y, pEXYZ[pg->iH[kk]].Z);
+      printf("\n");
     }
-    //        pg->t      = (Float_t)cryst[icry].CCT0; //OFT is not giving me the information on which is the detector of the first interaction
-    //        pg->tstamp = cryst[icry].tstamp; //OFT is not giving me the information on which is the detector of the first interaction
-    //int ipgX1 = (int)(pg->X1+0.5f);
-    //int ipgY1 = (int)(pg->Y1+0.5f);
-    //int ipgZ1 = (int)(pg->Z1+0.5f);
-    //bool found = false;
-    //for (int ihit = 0; ihit < number_of_hits; ihit++) {
-    //  if ((ipgX1 == (int)(pEXYZ[ihit].X+0.5f))&&
-    //      (ipgY1 == (int)(pEXYZ[ihit].Y+0.5f))&&
-    //      (ipgZ1 == (int)(pEXYZ[ihit].Z+0.5f))) {
-    //    for (int icry = 0; icry < number_of_crystals; icry++) {
-    //      if (cryst[icry].id == pEXYZ[ihit].Id) {
-    //        pg->t = (Float_t)cryst[icry].CCT0;
-    //        //pg->t += (cryst[icry].tstamp-timestampEvt)*10;  // relative to the smallest timestamp in event (10 ns/sample should be a symbol)
-    //        pg->tstamp = cryst[icry].tstamp;
-    //        found = true;
-    //        break;
-    //      }
-    //    }
-    //    if(found)
-    //      break;
-    //  }
-    //}
   }
 }

filters/Tracking/TrackingFilterOFT.cpp

@@ -11,18 +11,18 @@ using namespace ADF;
 
 #ifdef _MSC_VER  // *.c stuff compiled using C
 extern "C" {
-  void process_event(OftStruct *mydata, unsigned int *error_code);
-  void initialize_oft(double probLimit1, double probLimit2);
+  void oft_process_event(OftStruct *mydata, unsigned int *error_code);
+  void oft_initialize(double probLimit1, double probLimit2, double inner_r, double outer_r);
 }
 #else 
 	#if DO_OFT != ON // internal cmake force C++ compilation
 	extern "C" {
-		void process_event(OftStruct *mydata, unsigned int *error_code);
-		void initialize_oft(double probLimit1, double probLimit2);
+		void oft_process_event(OftStruct *mydata, unsigned int *error_code);
+		void oft_initialize(double probLimit1, double probLimit2, double inner_r, double outer_r);
 	}
 	#else // *.c stuff compiled using C++
-		void process_event(OftStruct *mydata, unsigned int *error_code);
-		void initialize_oft(double probLimit1, double probLimit2);
+		void oft_process_event(OftStruct *mydata, unsigned int *error_code);
+		void oft_initialize(double probLimit1, double probLimit2, double inner_r, double outer_r);
 	#endif
 #endif
 
@@ -44,7 +44,7 @@ Int_t TrackingFilterOFT::Process()
 
   FillOFTStructures();
 
-  process_event(local_ptr, &error_code);
+  oft_process_event(local_ptr, &error_code);
 
   MoveOFTStructures();
 
@@ -57,59 +57,30 @@ Int_t TrackingFilterOFT::AlgoSpecificInitialise()
 
   cout << endl;
 
-  local_ptr = (OftStruct *) malloc( sizeof(OftStruct) );
+  local_ptr = (OftStruct *) calloc(1, sizeof(OftStruct) );
 
-  local_ptr->etot = (double *)calloc(intmax, sizeof(double));
-
-  local_ptr->pos    = (posi3D *)calloc(intmax, sizeof(posi3D));
-  local_ptr->first  = (posi3D *)calloc(intmax, sizeof(posi3D));
-  local_ptr->second = (posi3D *)calloc(intmax, sizeof(posi3D));
-
-  local_ptr->probtot = (double *)calloc(intmax, sizeof(double));
-
-  local_ptr->flagu = (int *)calloc(intmax, sizeof(int));
-  local_ptr->nbtot = (int *)calloc(intmax, sizeof(int));
-
-  local_ptr->e    = (double *)calloc(intmax, sizeof(double));
-  local_ptr->et   = (double *)calloc(intmax, sizeof(double));
-
-  local_ptr->r    = (double *)calloc(intmax*intmax, sizeof(double));
-  local_ptr->r_ge = (double *)calloc(intmax*intmax, sizeof(double));
-
-  local_ptr->sn          = (int *)calloc(kmax*intmax, sizeof(int));
-  local_ptr->interaction = (int *)calloc(kmax*intmax, sizeof(int));
-
-  //local_ptr->seg = (int *)calloc(intmax, sizeof(int));
-  //local_ptr->det = (int *)calloc(intmax, sizeof(int));
-
-  local_ptr->angtheta = (double *)calloc(intmax, sizeof(double));
-  local_ptr->angphi   = (double *)calloc(intmax, sizeof(double));
-  local_ptr->angn     = (double *)calloc(intmax, sizeof(double));
-
-  local_ptr->numn = (int *)calloc(intmax, sizeof(int));
-
-  local_ptr->source.x = fSource[0]/10.;
-  local_ptr->source.y = fSource[1]/10;
-  local_ptr->source.z = fSource[2]/10;
+  local_ptr->source.x = 0.1*fSource[0];
+  local_ptr->source.y = 0.1*fSource[1];
+  local_ptr->source.z = 0.1*fSource[2];
 
   istringstream dd(fGeometrySummary);
-  double inner_r = 0, outer_r = 0;
-  dd >> inner_r >> outer_r;
-  if(!dd.good() || inner_r<=0 || outer_r<=0 || outer_r<=inner_r) {
+  double r_inner = 0, r_outer = 0;
+  dd >> r_inner >> r_outer;
+  if(!dd.good() || r_inner<=0 || r_outer<=0 || r_outer<=r_inner) {
     cout << "Error reading the AGATA summary ( " << fGeometrySummary << " ) " << endl;
     cout << "Back to default values" << endl;
-    inner_r = 235.;
-    outer_r = 329.;
+    r_inner = 235.;
+    r_outer = 329.;
   }
-  local_ptr->inner_r = inner_r/10.;  // inner radius of AGATA or AGATA demonstrator 
-  local_ptr->outer_r = outer_r/10.;  // outer radius of AGATA or AGATA demonstrator 
+  r_inner /= 10.;  // inner radius of AGATA or AGATA demonstrator (cm)
+  r_outer /= 10.;  // outer radius of AGATA or AGATA demonstrator (cm) 
 
-  cout << "AGATA inner and outer radius : " << local_ptr->inner_r << "  " << local_ptr->outer_r <<  " (cm)" << endl;
+  cout << "AGATA inner and outer radius : " << r_inner << "  " << r_outer <<  " (cm)" << endl;
 
-  local_ptr->mult   = 0;
+  local_ptr->nb_gam = 0;
   local_ptr->nb_int = 0;
 
-  initialize_oft(fAcceptanceValue1, fAcceptanceValue2);
+  oft_initialize(fAcceptanceValue1, fAcceptanceValue2, r_inner, r_outer);
 
   return 0;
 }
@@ -120,14 +91,17 @@ void TrackingFilterOFT::FillOFTStructures()
   char tag2[] = "    ";
   char * tag = tag1;
 
+  local_ptr->source.x = 0.1*fSource[0];
+  local_ptr->source.y = 0.1*fSource[1];
+  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++) {