Adding a reminder for coders
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALReconstructor.cxx
index e41a446..b3fba39 100644 (file)
 #include "AliEMCALRawUtils.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCALClusterizerv1.h"
+#include "AliEMCALClusterizerv2.h"
 #include "AliEMCALClusterizerNxN.h"
 #include "AliEMCALRecPoint.h"
 #include "AliEMCALPID.h"
-#include "AliEMCALTrigger.h"
 #include "AliRawReader.h"
 #include "AliCDBEntry.h"
 #include "AliCDBManager.h"
@@ -67,7 +67,7 @@
 #include "AliEMCALTriggerTypes.h"
 
 ClassImp(AliEMCALReconstructor) 
-
+  
 const AliEMCALRecParam*     AliEMCALReconstructor::fgkRecParam        = 0;   // EMCAL rec. parameters
 AliEMCALRawUtils*           AliEMCALReconstructor::fgRawUtils         = 0;   // EMCAL raw utilities class
 AliEMCALClusterizer*        AliEMCALReconstructor::fgClusterizer      = 0;   // EMCAL clusterizer class
@@ -77,10 +77,12 @@ TClonesArray*               AliEMCALReconstructor::fgTriggerDigits    = 0;   //
 AliEMCALTriggerElectronics* AliEMCALReconstructor::fgTriggerProcessor = 0x0;
 //____________________________________________________________________________
 AliEMCALReconstructor::AliEMCALReconstructor() 
-  : fGeom(0),fCalibData(0),fPedestalData(0),fTriggerData(0x0) 
+  : fGeom(0),fCalibData(0),fPedestalData(0),fTriggerData(0x0), fMatches(0x0)
 {
   // ctor
 
+  // AliDebug(2, "Mark.");  
+
   fgRawUtils = new AliEMCALRawUtils;
   
   //To make sure we match with the geometry in a simulation file,
@@ -133,13 +135,19 @@ AliEMCALReconstructor::AliEMCALReconstructor()
   fgDigitsArr     = new TClonesArray("AliEMCALDigit",1000);
   fgClustersArr   = new TObjArray(1000);
   fgTriggerDigits = new TClonesArray("AliEMCALTriggerRawDigit",1000);  
+
+  //Track matching
+  fMatches = new TList();
+  fMatches->SetOwner(kTRUE);
 } 
 
 //____________________________________________________________________________
 AliEMCALReconstructor::~AliEMCALReconstructor()
 {
   // dtor
-  
+
+  //AliDebug(2, "Mark.");
+
   if(fGeom)              delete fGeom;
   
   //No need to delete, recovered from OCDB
@@ -165,6 +173,8 @@ AliEMCALReconstructor::~AliEMCALReconstructor()
   if(fgClusterizer)      delete fgClusterizer;
   if(fgTriggerProcessor) delete fgTriggerProcessor;
   
+  if(fMatches) { fMatches->Delete(); delete fMatches; fMatches = 0;}
+  
   AliCodeTimer::Instance()->Print();
 } 
 
@@ -194,24 +204,34 @@ void AliEMCALReconstructor::InitClusterizer() const
       //printf("ReCreate clusterizer? Clusterizer set <%d>, Clusterizer in use <%s>\n",
       //     clusterizerType, fgClusterizer->Version());
       
-      if     (clusterizerType == AliEMCALRecParam::kClusterizerv1 && !strcmp(fgClusterizer->Version(),"clu-v1")) return;
+      if     (clusterizerType == AliEMCALRecParam::kClusterizerv1  && !strcmp(fgClusterizer->Version(),"clu-v1"))  return;
       
       else if(clusterizerType == AliEMCALRecParam::kClusterizerNxN && !strcmp(fgClusterizer->Version(),"clu-NxN")) return;
       
+      else if(clusterizerType == AliEMCALRecParam::kClusterizerv2  && !strcmp(fgClusterizer->Version(),"clu-v2"))  return;
+      
       //Need to create new clusterizer, the one set previously is not the correct one     
       delete fgClusterizer;
     }
     else return;
   }
   
-  if (clusterizerType  == AliEMCALRecParam::kClusterizerv1)
+  if      (clusterizerType  == AliEMCALRecParam::kClusterizerv1)
     {
-      fgClusterizer = new AliEMCALClusterizerv1(fGeom, fCalibData,fPedestalData);
+      fgClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData,fPedestalData);
     }
-  else
+  else if (clusterizerType  == AliEMCALRecParam::kClusterizerNxN)
     {
       fgClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData,fPedestalData);
     }
+  else if (clusterizerType  == AliEMCALRecParam::kClusterizerv2)
+  {
+    fgClusterizer = new AliEMCALClusterizerv2   (fGeom, fCalibData,fPedestalData);
+  }
+  else 
+  {
+    AliFatal(Form("Unknown clusterizer %d ", clusterizerType));
+  }
 }
 
 //____________________________________________________________________________
@@ -249,6 +269,9 @@ void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree)
   }//not a LED event
   
   clustersTree->Fill();        
+
+  // Deleting the recpoints at the end of the reconstruction call
+  fgClusterizer->DeleteRecPoints();
 }
 
 //____________________________________________________________________________
@@ -279,7 +302,7 @@ void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
     fgRawUtils->SetOption(GetOption());
     
     //  fgRawUtils->SetRawFormatHighLowGainFactor(GetRecParam()->GetHighLowGainFactor());
-  
+    
     //   fgRawUtils->SetRawFormatOrder(GetRecParam()->GetOrderParameter());
     //    fgRawUtils->SetRawFormatTau(GetRecParam()->GetTau());
     fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
@@ -323,7 +346,7 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   //########################################
   
   static int saveOnce = 0;
-       
+  
   Int_t v0M[2] = {0, 0};
   
   AliESDVZERO* esdV0 = esd->GetVZEROData();
@@ -383,22 +406,22 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
       trgESD->SetL1Threshold(0, fTriggerData->GetL1GammaThreshold());
       
       trgESD->SetL1Threshold(1, fTriggerData->GetL1JetThreshold()  );
-
-         Int_t v0[2];
-         fTriggerData->GetL1V0(v0);
-               
-         trgESD->SetL1V0(v0);  
-         trgESD->SetL1FrameMask(fTriggerData->GetL1FrameMask());            
-
-         if (!saveOnce) 
-         {
-                 int type[8] = {0};
-                 fTriggerData->GetL1TriggerType(type);
-
-                 esd->SetCaloTriggerType(type);
-                 
-                 saveOnce = 1;
-         }
+      
+      Int_t v0[2];
+      fTriggerData->GetL1V0(v0);
+      
+      trgESD->SetL1V0(v0);     
+      trgESD->SetL1FrameMask(fTriggerData->GetL1FrameMask());            
+      
+      if (!saveOnce && fTriggerData->GetL1DataDecoded()) 
+       {
+         int type[8] = {0};
+         fTriggerData->GetL1TriggerType(type);
+         
+         esd->SetCaloTriggerType(type);
+         
+         saveOnce = 1;
+       }
     }
   
   // Resetting
@@ -413,18 +436,20 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   
   Int_t nDigits = fgDigitsArr->GetEntries(), idignew = 0 ;
   AliDebug(1,Form("%d digits",nDigits));
-  
   AliESDCaloCells &emcCells = *(esd->GetEMCALCells());
   emcCells.CreateContainer(nDigits);
   emcCells.SetType(AliVCaloCells::kEMCALCell);
   Float_t energy = 0;
+  Float_t time   = 0;
   for (Int_t idig = 0 ; idig < nDigits ; idig++) {
     const AliEMCALDigit * dig = (const AliEMCALDigit*)fgDigitsArr->At(idig);
-    if(dig->GetAmplitude() > 0 ){
-      energy = fgClusterizer->Calibrate(dig->GetAmplitude(),dig->GetTime(),dig->GetId()); //TimeR or Time?
-      if(energy > 0){ //Digits tagged as bad (dead, hot, not alive) are set to 0 in calibrate, remove them     
-       emcCells.SetCell(idignew,dig->GetId(),energy, dig->GetTime());   
-       idignew++;
+    time   = dig->GetTime();      // Time already calibrated in clusterizer
+    energy = dig->GetAmplitude(); // energy calibrated in clusterizer
+    if(energy > 0 ){
+      fgClusterizer->Calibrate(energy,time,dig->GetId()); //Digits already calibrated in clusterizers
+      if(energy > 0){ //Digits tagged as bad (dead, hot, not alive) are set to 0 in calibrate, remove them
+        emcCells.SetCell(idignew,dig->GetId(),energy, time);   
+        idignew++;
       }
     }
   }
@@ -444,31 +469,14 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   
   Int_t nClusters = fgClustersArr->GetEntries(),  nClustersNew=0;
   AliDebug(1,Form("%d clusters",nClusters));
-  
-  //######################################################
-  //#######################TRACK MATCHING###############
-  //######################################################
-  //Fill list of integers, each one is index of track to which the cluster belongs.
-  
-  // step 1 - initialize array of matched track indexes
-  Int_t *matchedTrack = new Int_t[nClusters];
-  for (Int_t iclus = 0; iclus < nClusters; iclus++)
-    matchedTrack[iclus] = -1;  // neg. index --> no matched track
-  
-  // step 2, change the flag for all matched clusters found in tracks
-  Int_t iemcalMatch = -1;
-  Int_t endtpc = esd->GetNumberOfTracks();
-  for (Int_t itrack = 0; itrack < endtpc; itrack++) {
-    AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
-    iemcalMatch = track->GetEMCALcluster();
-    if(iemcalMatch >= 0) matchedTrack[iemcalMatch] = itrack;
-  } 
+
   
   //########################################
   //##############Fill CaloClusters#############
   //########################################
   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
     const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)fgClustersArr->At(iClust);
+    if(!clust) continue;
     //if(clust->GetClusterType()== AliVCluster::kEMCALClusterv1) nRP++; else nPC++;
     // clust->Print(); //For debugging
     // Get information from EMCAL reconstruction points
@@ -490,16 +498,16 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
     Int_t newCellMult = 0;
     for (Int_t iCell=0; iCell<cellMult; iCell++) {
       if (amplFloat[iCell] > 0) {
-       absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
-       //Calculate Fraction
-       if(emcCells.GetCellAmplitude(digitInts[iCell])>0 && GetRecParam()->GetUnfold()){
-         fracList[newCellMult] = amplFloat[iCell]/(emcCells.GetCellAmplitude(digitInts[iCell]));//get cell calibration value 
-         
-       }
-       else{
-         fracList[newCellMult] = 0; 
-       }
-       newCellMult++;
+        absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
+        //Calculate Fraction
+        if(emcCells.GetCellAmplitude(digitInts[iCell])>0 && GetRecParam()->GetUnfold()){
+          fracList[newCellMult] = amplFloat[iCell]/(emcCells.GetCellAmplitude(digitInts[iCell]));//get cell calibration value 
+          
+        }
+        else{
+          fracList[newCellMult] = 0; 
+        }
+        newCellMult++;
       }
     }
     
@@ -536,22 +544,74 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
       ec->SetM20(elipAxis[1]*elipAxis[1]) ;
       ec->SetTOF(clust->GetTime()) ; //time-of-fligh
       ec->SetNExMax(clust->GetNExMax());          //number of local maxima
-      TArrayI arrayTrackMatched(1);// Only one track, temporal solution.
-      arrayTrackMatched[0]= matchedTrack[iClust];
-      ec->AddTracksMatched(arrayTrackMatched);
+  
       
       TArrayI arrayParents(parentMult,parentList);
       ec->AddLabels(arrayParents);
-      
-      // add the cluster to the esd object
+      //
+      //Track matching
+      //
+      fMatches->Clear();
+      Int_t nTracks = esd->GetNumberOfTracks();
+      for (Int_t itrack = 0; itrack < nTracks; itrack++)
+       {
+         AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
+         if(track->GetEMCALcluster()==iClust)
+           {
+             Double_t dEta=-999, dPhi=-999;
+             Bool_t isMatch =  CalculateResidual(track, ec, dEta, dPhi);
+             if(!isMatch) 
+               {
+                 // AliDebug(10, "Not good");
+                 continue;
+               }
+             AliEMCALMatch *match = new AliEMCALMatch();
+             match->SetIndexT(itrack);
+             match->SetDistance(TMath::Sqrt(dEta*dEta+dPhi*dPhi));
+             match->SetdEta(dEta);
+             match->SetdPhi(dPhi);
+             fMatches->Add(match);
+           }
+       } 
+      fMatches->Sort(kSortAscending); //Sort matched tracks from closest to furthest
+      Int_t nMatch = fMatches->GetEntries();
+      TArrayI arrayTrackMatched(nMatch);
+      for(Int_t imatch=0; imatch<nMatch; imatch++)
+       {
+         AliEMCALMatch *match = (AliEMCALMatch*)fMatches->At(imatch);
+         arrayTrackMatched[imatch] = match->GetIndexT();
+         if(imatch==0)
+           {
+             ec->SetTrackDistance(match->GetdPhi(), match->GetdEta());
+           }
+       }
+      ec->AddTracksMatched(arrayTrackMatched);
+    
+      //add the cluster to the esd object
       esd->AddCaloCluster(ec);
+
       delete ec;
       delete [] newAbsIdList ;
       delete [] newFracList ;
     }
   } // cycle on clusters
+
+  //
+  //Reset the index of matched cluster for tracks
+  //to the one in CaloCluster array
+  Int_t ncls = esd->GetNumberOfCaloClusters();
+  for(Int_t icl=0; icl<ncls; icl++)
+    {
+      AliESDCaloCluster *cluster = esd->GetCaloCluster(icl);
+      if(!cluster || !cluster->IsEMCAL()) continue;
+      TArrayI *trackIndex = cluster->GetTracksMatched();
+      for(Int_t itr=0; itr<trackIndex->GetSize(); itr++)
+       {
+         AliESDtrack *track = esd->GetTrack(trackIndex->At(itr));
+         track->SetEMCALcluster(cluster->GetID());
+       }
+    }
   
-  delete [] matchedTrack;
   
   //Fill ESDCaloCluster with PID weights
   AliEMCALPID *pid = new AliEMCALPID;
@@ -629,4 +689,97 @@ void AliEMCALReconstructor::ReadDigitsArrayFromTree(TTree *digitsTree) const
   branch->GetEntry(0);
 }
 
+//==================================================================================
+Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCluster *cluster, Double_t &dEta, Double_t &dPhi)const
+{
+  //
+  // calculate the residual between track and cluster
+  //
+
+  // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
+  // Otherwise use the TPCInner point
+  AliExternalTrackParam *trkParam = 0;
+  const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
+  if(friendTrack && friendTrack->GetTPCOut())
+    trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
+  else
+    trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
+  if(!trkParam) return kFALSE;
+
+  //Perform extrapolation
+  Double_t trkPos[3];
+  Float_t  clsPos[3];
+
+  AliExternalTrackParam trkParamTmp (*trkParam);
+  cluster->GetPosition(clsPos);
+  TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
+  Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
+  //Rotate to proper local coordinate
+  vec.RotateZ(-alpha); 
+  trkParamTmp.Rotate(alpha); 
+  //extrapolation is done here
+  if(!AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, vec.X(), track->GetMass(), GetRecParam()->GetExtrapolateStep(), kFALSE, 0.8, -1)) 
+    return kFALSE; 
+
+  //Calculate the residuals
+  trkParamTmp.GetXYZ(trkPos); 
+   
+  TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
+  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+      
+  Double_t clsPhi = clsPosVec.Phi();
+  if(clsPhi<0) clsPhi+=2*TMath::Pi();
+  Double_t trkPhi = trkPosVec.Phi();
+  if(trkPhi<0) trkPhi+=2*TMath::Pi();
+
+  dPhi = clsPhi-trkPhi;
+  dEta = clsPosVec.Eta()-trkPosVec.Eta();
+
+  return kTRUE;
+}
+
+//
+//==================================================================================
+//
+AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch() 
+  : TObject(),  
+    fIndexT(-1), 
+    fDistance(-999.),
+    fdEta(-999.),
+    fdPhi(-999.)
+{
+  //default constructor
 
+}
+
+//
+//==================================================================================
+//
+AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
+  : TObject(),
+    fIndexT(copy.fIndexT),
+    fDistance(copy.fDistance),
+    fdEta(copy.fdEta),
+    fdPhi(copy.fdPhi)
+{
+  //copy ctor
+}
+
+//
+//==================================================================================
+//
+Int_t AliEMCALReconstructor::AliEMCALMatch::Compare(const TObject *obj) const 
+{
+  //
+  // Compare wrt the residual
+  //
+       
+  AliEMCALReconstructor::AliEMCALMatch *that = (AliEMCALReconstructor::AliEMCALMatch*)obj;
+       
+  Double_t thisDist = fDistance;//fDistance;
+  Double_t thatDist = that->fDistance;//that->GetDistance();
+       
+  if (thisDist > thatDist) return 1;
+  else if (thisDist < thatDist) return -1;
+  return 0;
+}