]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALReconstructor.cxx
do not delete fGeom in ~AliEMCALReconstructor since it is owned by AliRun. When delet...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALReconstructor.cxx
index e88708e1fa04be1a8d0c219051c9a225abf87402..9d3dc7239c166d0076c101d38dd2c6f495e0b02e 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 "AliEMCALRecoUtils.h"
 #include "AliRawReader.h"
 #include "AliCDBEntry.h"
 #include "AliCDBManager.h"
@@ -80,7 +81,9 @@ AliEMCALReconstructor::AliEMCALReconstructor()
   : 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,
@@ -100,7 +103,7 @@ AliEMCALReconstructor::AliEMCALReconstructor()
   if(!fCalibData)
     {
       AliCDBEntry *entry = (AliCDBEntry*) 
-       AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+      AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
     }
   
@@ -111,7 +114,7 @@ AliEMCALReconstructor::AliEMCALReconstructor()
   if(!fPedestalData)
     {
       AliCDBEntry *entry = (AliCDBEntry*) 
-       AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
+      AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
       if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
     }
   
@@ -132,8 +135,10 @@ AliEMCALReconstructor::AliEMCALReconstructor()
   //Init temporary list of digits
   fgDigitsArr     = new TClonesArray("AliEMCALDigit",1000);
   fgClustersArr   = new TObjArray(1000);
-  fgTriggerDigits = new TClonesArray("AliEMCALTriggerRawDigit",1000);  
 
+  const int kNTRU = fGeom->GetNTotalTRU();
+  fgTriggerDigits = new TClonesArray("AliEMCALTriggerRawDigit", kNTRU * 96);   
+       
   //Track matching
   fMatches = new TList();
   fMatches->SetOwner(kTRUE);
@@ -143,31 +148,34 @@ AliEMCALReconstructor::AliEMCALReconstructor()
 AliEMCALReconstructor::~AliEMCALReconstructor()
 {
   // dtor
-  
-  if(fGeom)              delete fGeom;
+
+  //AliDebug(2, "Mark.");
+
+  ////RS  if(fGeom)              delete fGeom;
   
   //No need to delete, recovered from OCDB
   //if(fCalibData)         delete fCalibData;
   //if(fPedestalData)      delete fPedestalData;
   
-  if(fgDigitsArr){
-    fgDigitsArr->Clear("C");
-    delete fgDigitsArr; 
-  }
+  if(fgDigitsArr) fgDigitsArr->Clear("C");
+  delete fgDigitsArr; 
+  fgDigitsArr = 0;
   
-  if(fgClustersArr){
-    fgClustersArr->Clear();
-    delete fgClustersArr; 
-  }
+  if(fgClustersArr) fgClustersArr->Clear();
+  delete fgClustersArr; 
+  fgClustersArr = 0;
   
-  if(fgTriggerDigits){
-    fgTriggerDigits->Clear();
-    delete fgTriggerDigits; 
-  }
+  if(fgTriggerDigits) fgTriggerDigits->Clear();
+  delete fgTriggerDigits; 
+  fgTriggerDigits = 0;
+  
+  delete fgRawUtils;
+  fgRawUtils = 0;
+  delete fgClusterizer;
+  fgClusterizer = 0;
   
-  if(fgRawUtils)         delete fgRawUtils;
-  if(fgClusterizer)      delete fgClusterizer;
-  if(fgTriggerProcessor) delete fgTriggerProcessor;
+  delete fgTriggerProcessor;
+  fgTriggerProcessor = 0;
   
   if(fMatches) { fMatches->Delete(); delete fMatches; fMatches = 0;}
   
@@ -200,24 +208,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));
+  }
 }
 
 //____________________________________________________________________________
@@ -255,6 +273,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();
 }
 
 //____________________________________________________________________________
@@ -271,7 +292,8 @@ void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
   
   if(fgDigitsArr) fgDigitsArr->Clear("C");
   
-  TClonesArray *digitsTrg = new TClonesArray("AliEMCALTriggerRawDigit", 32 * 96);
+  const int kNTRU = fGeom->GetNTotalTRU();
+  TClonesArray *digitsTrg = new TClonesArray("AliEMCALTriggerRawDigit", kNTRU * 96);
   
   Int_t bufsize = 32000;
   digitsTree->Branch("EMCAL", &fgDigitsArr, bufsize);
@@ -291,9 +313,10 @@ void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
     fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
     fgRawUtils->SetNPedSamples(GetRecParam()->GetNPedSamples());
     fgRawUtils->SetRemoveBadChannels(GetRecParam()->GetRemoveBadChannels());
-    fgRawUtils->SetFittingAlgorithm(GetRecParam()->GetFittingAlgorithm());
+    if (!fgRawUtils->GetFittingAlgorithm()) fgRawUtils->SetFittingAlgorithm(GetRecParam()->GetFittingAlgorithm());
     fgRawUtils->SetFALTROUsage(GetRecParam()->UseFALTRO());
-    
+    //  fgRawUtils->SetFALTROUsage(0);
     //fgRawUtils->SetTimeMin(GetRecParam()->GetTimeMin());
     //fgRawUtils->SetTimeMax(GetRecParam()->GetTimeMax());
     
@@ -336,18 +359,15 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   
   if (esdV0) 
     {
-      for (Int_t i = 0; i < 32; i++)
-       {
-         v0M[0] += (Int_t)esdV0->GetAdcV0C(i);
-         v0M[1] += (Int_t)esdV0->GetAdcV0A(i);
-       }
+               v0M[0] = esdV0->GetTriggerChargeC();
+               v0M[1] = esdV0->GetTriggerChargeA();
     }
   else
     {
-      AliWarning("Cannot retrieve V0 ESD! Run w/ null V0 charges");
+      AliWarning("No V0 ESD! Run trigger processor w/ null V0 charges");
     }
   
-  if (fgTriggerDigits) fgTriggerDigits->Clear();
+  if (fgTriggerDigits && fgTriggerDigits->GetEntriesFast()) fgTriggerDigits->Delete();
   
   TBranch *branchtrg = digitsTree->GetBranch("EMTRG");
   
@@ -373,7 +393,8 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
       for (Int_t i = 0; i < fgTriggerDigits->GetEntriesFast(); i++)
        {         
          AliEMCALTriggerRawDigit* rdig = (AliEMCALTriggerRawDigit*)fgTriggerDigits->At(i);
-         
+         if (AliDebugLevel() > 999) rdig->Print("");
+               
          Int_t px, py;
          if (fGeom->GetPositionInEMCALFromAbsFastORIndex(rdig->GetId(), px, py))
            {
@@ -381,14 +402,15 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
              
              rdig->GetMaximum(a, t);
              rdig->GetL0Times(times);
-             
+                       
              trgESD->Add(px, py, a, t, times, rdig->GetNL0Times(), rdig->GetL1TimeSum(), rdig->GetTriggerBits());
            }
        }
       
-      trgESD->SetL1Threshold(0, fTriggerData->GetL1GammaThreshold());
-      
-      trgESD->SetL1Threshold(1, fTriggerData->GetL1JetThreshold()  );
+               for (int i = 0; i < 2; i++) {
+                       trgESD->SetL1Threshold(2 * i    , fTriggerData->GetL1JetThreshold(  i));
+                       trgESD->SetL1Threshold(2 * i + 1, fTriggerData->GetL1GammaThreshold(i));
+               }
       
       Int_t v0[2];
       fTriggerData->GetL1V0(v0);
@@ -398,7 +420,7 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
       
       if (!saveOnce && fTriggerData->GetL1DataDecoded()) 
        {
-         int type[8] = {0};
+         int type[15] = {0};
          fTriggerData->GetL1TriggerType(type);
          
          esd->SetCaloTriggerType(type);
@@ -423,17 +445,50 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   AliESDCaloCells &emcCells = *(esd->GetEMCALCells());
   emcCells.CreateContainer(nDigits);
   emcCells.SetType(AliVCaloCells::kEMCALCell);
+  
   Float_t energy = 0;
-  for (Int_t idig = 0 ; idig < nDigits ; idig++) {
+  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());   
+    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
+      { 
+        // Only for MC
+        // Get the label of the primary particle that generated the cell
+        // Assign the particle that deposited more energy
+        Int_t   nprimaries = dig->GetNprimary() ;
+        Int_t   digLabel   =-1 ;
+        Float_t edep       =-1.;
+        if ( nprimaries > 0 )
+        {
+          Int_t jndex ;
+          for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit
+            if(edep < dig->GetDEParent(jndex+1))
+            {
+              digLabel = dig->GetIparent (jndex+1);
+              edep     = dig->GetDEParent(jndex+1);
+            }
+                   
+          } // all primaries in digit      
+        } // select primary label
+        
+        Bool_t highGain = kFALSE;
+        if( dig->GetType() == AliEMCALDigit::kHG ) highGain = kTRUE;
+        
+        emcCells.SetCell(idignew,dig->GetId(),energy, time,digLabel,0.,highGain);
         idignew++;
       }
     }
   }
+  
   emcCells.SetNumberOfCells(idignew);
   emcCells.Sort();
   
@@ -539,12 +594,12 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
          AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
          if(track->GetEMCALcluster()==iClust)
            {
-             Double_t dEta=-999, dPhi=-999;
+             Float_t dEta=-999, dPhi=-999;
              Bool_t isMatch =  CalculateResidual(track, ec, dEta, dPhi);
              if(!isMatch) 
                {
+                 // AliDebug(10, "Not good");
                  continue;
-                 cout<<"Not good"<<endl;
                }
              AliEMCALMatch *match = new AliEMCALMatch();
              match->SetIndexT(itrack);
@@ -625,10 +680,25 @@ void AliEMCALReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{
   const Int_t bufsize = 255;
   char path[bufsize] ;
   TGeoHMatrix * m = 0x0;
+  Int_t tmpType = -1;
+  Int_t SMOrder = 0;
+  TString SMName;
   for(Int_t sm = 0; sm < fGeom->GetNumberOfSuperModules(); sm++){
-    snprintf(path,bufsize,"/ALIC_1/XEN1_1/SMOD_%d",sm+1) ; //In Geometry modules numbered 1,2,.,5
-    if(sm >= 10) snprintf(path,bufsize,"/ALIC_1/XEN1_1/SM10_%d",sm-10+1) ;
-    
+    if(fGeom->GetSMType(sm) == AliEMCALGeometry::kEMCAL_Standard )      SMName = "SMOD";
+    else if(fGeom->GetSMType(sm) == AliEMCALGeometry::kEMCAL_Half )     SMName = "SM10";
+    else if(fGeom->GetSMType(sm) == AliEMCALGeometry::kEMCAL_3rd )      SMName = "SM3rd";
+    else if( fGeom->GetSMType(sm) == AliEMCALGeometry::kDCAL_Standard ) SMName = "DCSM";
+    else if( fGeom->GetSMType(sm) == AliEMCALGeometry::kDCAL_Ext )      SMName = "DCEXT";
+    else AliError("Unkown SM Type!!");
+
+    if(fGeom->GetSMType(sm) == tmpType) {
+      SMOrder++;
+    } else {
+      tmpType = fGeom->GetSMType(sm);
+      SMOrder = 1;
+    }
+    snprintf(path,bufsize,"/ALIC_1/XEN1_1/%s_%d", SMName.Data(), SMOrder) ;
+
     if (gGeoManager->CheckPath(path)){
       gGeoManager->cd(path);
       m = gGeoManager->GetCurrentMatrix() ;
@@ -671,7 +741,7 @@ void AliEMCALReconstructor::ReadDigitsArrayFromTree(TTree *digitsTree) const
 }
 
 //==================================================================================
-Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCluster *cluster, Double_t &dEta, Double_t &dPhi)const
+Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCluster *cluster, Float_t &dEta, Float_t &dPhi)const
 {
   //
   // calculate the residual between track and cluster
@@ -679,44 +749,30 @@ Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCl
 
   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
   // Otherwise use the TPCInner point
-  AliExternalTrackParam *trkParam;
+
+  dEta = -999, dPhi = -999;
+  Bool_t ITSTrackSA = 0;
+
+  AliExternalTrackParam *trkParam = 0;
+  
   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
   if(friendTrack && friendTrack->GetTPCOut())
     trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
-  else
+  else if(track->GetInnerParam())
     trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
+  else{
+    trkParam = new AliExternalTrackParam(*track); //If there is ITSSa track 
+    ITSTrackSA = 1;    
+  }
   if(!trkParam) return kFALSE;
+  
+  AliExternalTrackParam trkParamTmp (*trkParam);
+  if(!AliEMCALRecoUtils::ExtrapolateTrackToCluster(&trkParamTmp, cluster, track->GetMass(kTRUE), GetRecParam()->GetExtrapolateStep(), dEta, dPhi)){
+       if(ITSTrackSA) delete trkParam;
+       return kFALSE;
+  }
 
-  //Perform extrapolation
-  Double_t trkPos[3];
-  Float_t  clsPos[3];
-
-  AliExternalTrackParam *trkParamTmp = new AliExternalTrackParam(*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);     
-  delete trkParamTmp;
-   
-  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();
-
+  if(ITSTrackSA) delete trkParam;
   return kTRUE;
 }
 
@@ -746,7 +802,14 @@ AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
 {
   //copy ctor
 }
+//_____________________________________________________________________
+AliEMCALReconstructor::AliEMCALMatch& AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch::operator = (const AliEMCALMatch &source)
+{ // assignment operator; use copy ctor
+  if (&source == this) return *this;
 
+  new (this) AliEMCALMatch(source);
+  return *this;
+}
 //
 //==================================================================================
 //