Fix issue with free() index. Added min track length

Author: rigel0586

EsbGeometry/EsbSuperFGD/EsbConfig/fgdconfig

@@ -11,7 +11,7 @@ number_cubes_Z=100
 
 #magnetic field [1 kGauss = 0.1 T]
 magField_X=0    #in [kG] 
-magField_Y=20    #in [kG] 
+magField_Y=10    #in [kG] 
 magField_Z=0    #in [kG] 
 
 #visualization
@@ -37,17 +37,20 @@ searchneighborhood=0
 # setting for lower photon limit
 errPhotoLimit=0   
 
+#exclude tracks below minimum track Lenght
+min_track_lenght=15
+
 # The belog settings for graph traversal
 
 # The grad in angle is calculated as the angle between two vectors.
 #   1st vector is the vector distance between 2 cubes at gradDist apart
 #   2nd vector is the vector distance between 2 cubes shifter gradIntervalDist number of cubes from the 1st one
-gradDist=7              # Distance in number of cubes for track to calculate the grad
-gradIntervalDist=4      # Interval between two distance vectors to calculate grad
+gradDist=6              # Distance in number of cubes for track to calculate the grad
+gradIntervalDist=3      # Interval between two distance vectors to calculate grad
 gradDiff=25             # Allowable difference between grad to consider a different track
 
 # Calculate momentum from the segment with length 'momTrackSegment' by averaging the calculated momentum from 3 points (the 1st of the segment,
 # the second and for each point in between). Since the are energy losses and the momentum decreases, take the particle momentum from as the average from the 
 # first'avgTrackMomentum' points
-momTrackSegment=15
-avgTrackMomentum=3
+momTrackSegment=11
+avgTrackMomentum=5

EsbGeometry/EsbSuperFGD/EsbFgdDetectorParameters.h

@@ -47,6 +47,8 @@ namespace DP
     static const string FGD_MC_LEPTON_RECONSTRUCTION_BRANCH = "FgdMCLeptonStatsBranch";
     static const string FGD_MC_LEPTON_RECONSTRUCTION_ROOT_FILE = "FgdMCLeptonStatsRootFile";
 
+    static const string FGD_MIN_TRACK_LENGTH = "min_track_lenght";
+
     static const string FGD_GRAD_DIST = "gradDist";
     static const string FGD_GRAD_INTERVAL_DIST = "gradIntervalDist";
     static const string FGD_GRAD_ALLOWABLE_DIFF = "gradDiff";

EsbReconstruction/EsbSuperFGD/FgdGraphStats.cxx

@@ -83,6 +83,7 @@ FgdGraphStats::FgdGraphStats(const char* name
                           , double debugLlv) :
   FairTask(name, verbose) 
   , fgdConstructor(geoConfigFile)
+  , fMinTrackLenght(0)
 { 
     fpdgDB = make_shared<TDatabasePDG>();
     fParams.LoadPartParams(geoConfigFile);
@@ -94,7 +95,11 @@ FgdGraphStats::FgdGraphStats(const char* name
 // -----   Destructor   ----------------------------------------------------
 FgdGraphStats::~FgdGraphStats() 
 {
-    
+    if(fHitArray) 
+    {
+      fHitArray->Delete();
+    }
+    delete fHitArray;
 }
 // -------------------------------------------------------------------------
 
@@ -120,6 +125,8 @@ InitStatus FgdGraphStats::Init()
   f_total_Y = f_step_Y * f_bin_Y;
   f_total_Z = f_step_Z * f_bin_Z;
 
+  fMinTrackLenght = fParams.ParamAsInt(esbroot::geometry::superfgd::DP::FGD_MIN_TRACK_LENGTH);
+
   // Get RootManager
   FairRootManager* manager = FairRootManager::Instance();
   if ( !manager ) {
@@ -159,14 +166,17 @@ void FgdGraphStats::Exec(Option_t* opt)
 
         std::vector<ReconHit> allGRhits;
         std::vector<std::vector<ReconHit*>> foundGRTracks;
+
         rc = rc && GetGraphHits(allGRhits);
 
+
         if(!rc) LOG(info) << " Could not get graph hits";
 
         if(rc)
         { 
             LOG(debug) <<" MC Hits to retrieve stats from " << allMChits.size();
-            SplitMCTrack(allMChits, foundMCTracks);
+            //SplitMCTrack(allMChits, foundMCTracks);
+            SplitMCTrack(allGRhits, foundMCTracks);
         }
 
 
@@ -199,46 +209,45 @@ void FgdGraphStats::Exec(Option_t* opt)
 
 void FgdGraphStats::CompareTracks(std::vector<std::vector<ReconHit>>& mcTracks, std::vector<std::vector<ReconHit*>>& grTracks)
 {
-    static int TRACK_MIN_SIZE = 3;
     // 1. Calculate number of trakcs in Graph
-    int numOfGRTracksBigThan3 = 0;
+    int numOfGRTracksBigThanLimit = 0;
     for(Int_t grInd = 0 ; grInd < grTracks.size(); ++grInd)
     {
         std::vector<ReconHit*>& track = grTracks[grInd];
-        if(track.size() > TRACK_MIN_SIZE) numOfGRTracksBigThan3++;
+        if(track.size() > fMinTrackLenght) numOfGRTracksBigThanLimit++;
     }
 
     // 2. Calculate number of tracks in MC (charged particles only)
-    int numOfMCTracksBigThan3 = 0;
+    int numOfMCTracksBigThanLimit = 0;
     for(Int_t mcInd = 0 ; mcInd < mcTracks.size(); ++mcInd)
     {
         std::vector<ReconHit> track = mcTracks[mcInd];
-        if(track.size() > TRACK_MIN_SIZE
+        if(track.size() > fMinTrackLenght
             && IsChargedParticle(track[0])) 
         {
-           numOfMCTracksBigThan3++; 
+           numOfMCTracksBigThanLimit++; 
         }
     }
 
     LOG(info) << "+++++++++++++++++++++++++++++++++";
-    LOG(info) << " Graph tracks [3 cubes at least] = " << numOfGRTracksBigThan3;
-    LOG(info) << " MC tracks [3 cubes at least/ charge is not zero of the particle] = " << numOfMCTracksBigThan3;
-    LOG(info) << " Number of Graph tracks / Number of MC tracks = " << (1.0 * numOfGRTracksBigThan3) / numOfMCTracksBigThan3;
+    LOG(info) << " Graph tracks ["<< fMinTrackLenght <<" cubes at least] = " << numOfGRTracksBigThanLimit;
+    LOG(info) << " MC tracks ["<< fMinTrackLenght <<" cubes at least/ charge is not zero of the particle] = " << numOfMCTracksBigThanLimit;
+    LOG(info) << " Number of Graph tracks / Number of MC tracks = " << (1.0 * numOfGRTracksBigThanLimit) / numOfMCTracksBigThanLimit;
     LOG(info) << "+++++++++++++++++++++++++++++++++";
 
     // 3. Calculate the best match for GR track to MC track
     for(Int_t grInd = 0 ; grInd < grTracks.size(); ++grInd)
     {
         std::vector<ReconHit*>& grtrack = grTracks[grInd];
-        if(grtrack.size() <= TRACK_MIN_SIZE) continue;
+        if(grtrack.size() <= fMinTrackLenght) continue;
 
         Double_t maxFit = 0;
         Int_t indOfMcMaxFit = 0;
 
         for(Int_t mcInd = 0 ; mcInd < mcTracks.size(); ++mcInd)
         {
             std::vector<ReconHit> mctrack = mcTracks[mcInd];
-            if(mctrack.size() > TRACK_MIN_SIZE
+            if(mctrack.size() > fMinTrackLenght
                 && IsChargedParticle(mctrack[0])) 
             {
                 Double_t res = CmpGrToMCTrack(mctrack, grtrack);
@@ -347,6 +356,7 @@ Bool_t FgdGraphStats::GetGraphHits(std::vector<ReconHit>& allHits)
   for(Int_t i =0; i <  fHitArray->GetEntries() ; i++)
   {
     data::superfgd::FgdHit* hit = (data::superfgd::FgdHit*)fHitArray->At(i);
+
     TVector3  photoE = hit->GetPhotoE();    
     TVector3  mppcLoc = hit->GetMppcLoc();  
 
@@ -363,11 +373,6 @@ Bool_t FgdGraphStats::GetGraphHits(std::vector<ReconHit>& allHits)
     if(visited[ind])
     {
       // If already exists, add the photons
-      ReconHit toFind;
-      toFind.fmppcLoc = mppcLoc;
-      std::vector<ReconHit>::iterator recHit = find(allHits.begin(), allHits.end(), toFind);
-      ReconHit& foundHit = *recHit;
-      foundHit.fphotons = foundHit.fphotons + photoE;
       continue;
     }
     visited[ind] = true;

EsbReconstruction/EsbSuperFGD/FgdGraphStats.h

@@ -121,12 +121,8 @@ class FgdGraphStats : public FairTask
     /** Path to the used media.geo file - containing definitions of materials **/
     std::string fmediaFile;//!<!
 
-    /** Path to the events file containing Monte carlo simulation data - pdg codes, momentums etc. **/
-    std::string feventFile;//!<!
+    Int_t fMinTrackLenght;
 
-    Int_t fminGenFitInterations;
-    Int_t fmaxGenFitIterations;
-    Int_t fminHits;
 
     /** Are materials already defined **/
     bool isDefinedMaterials;