Add the configuration for AD0

[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.cxx

diff --git a/EMCAL/AliEMCALRecoUtils.cxx b/EMCAL/AliEMCALRecoUtils.cxx
index 705746eb818c95c5ce9ef6cd5cecac91c198ac87..6e512eeb407b2bb0a10f49a6660cf41cc065e2c0 100644 (file)
--- a/EMCAL/AliEMCALRecoUtils.cxx
+++ b/EMCAL/AliEMCALRecoUtils.cxx
@@ -83,7 +83,7 @@ AliEMCALRecoUtils::AliEMCALRecoUtils():
   fCutMinNClusterITS(0),                  fCutMaxChi2PerClusterTPC(0),            fCutMaxChi2PerClusterITS(0),
   fCutRequireTPCRefit(kFALSE),            fCutRequireITSRefit(kFALSE),            fCutAcceptKinkDaughters(kFALSE),
   fCutMaxDCAToVertexXY(0),                fCutMaxDCAToVertexZ(0),                 fCutDCAToVertex2D(kFALSE),
-  fCutRequireITSStandAlone(kFALSE),       fCutRequireITSpureSA(kFALSE) 
+  fCutRequireITSStandAlone(kFALSE),       fCutRequireITSpureSA(kFALSE)
 {
 //
   // Constructor.
@@ -136,7 +136,7 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
   fCutRequireTPCRefit(reco.fCutRequireTPCRefit),             fCutRequireITSRefit(reco.fCutRequireITSRefit),
   fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),     fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),    
   fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),             fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
-  fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone),   fCutRequireITSpureSA(reco.fCutRequireITSpureSA) 
+  fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone),   fCutRequireITSpureSA(reco.fCutRequireITSpureSA)
 {
   //Copy ctor
   
@@ -222,7 +222,8 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   fCutMaxDCAToVertexZ        = reco.fCutMaxDCAToVertexZ;
   fCutDCAToVertex2D          = reco.fCutDCAToVertex2D;
   fCutRequireITSStandAlone   = reco.fCutRequireITSStandAlone; 
-  fCutRequireITSpureSA       = reco.fCutRequireITSpureSA; 
+  fCutRequireITSpureSA       = reco.fCutRequireITSpureSA;
+
   if (reco.fResidualEta) {
     // assign or copy construct
     if (fResidualEta) { 
@@ -270,7 +271,7 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
     if (fMatchedClusterIndex) delete fMatchedClusterIndex;
     fMatchedClusterIndex = 0;
   }
-   
+  
   return *this;
 }
 
@@ -377,6 +378,7 @@ Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom,
   if (iSM < 0 || iphi < 0 || ieta < 0 ) {
     AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
                   iSM,ieta,iphi));
+    return kFALSE; // trick coverity
   }
   
   //Check rows/phi
@@ -695,6 +697,36 @@ Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster)
       break;
     }
       
+    case kSDMv5:
+    {
+      //Based on fit to the MC/data using kNoCorrection on the data - utilizes symmetric decay method and kPi0MCv5(MC) - 28 Oct 2013
+      //fNonLinearityParams[0] =  1.0;
+      //fNonLinearityParams[1] =  6.64778e-02;
+      //fNonLinearityParams[2] =  1.570;
+      //fNonLinearityParams[3] =  9.67998e-02;
+      //fNonLinearityParams[4] =  2.19381e+02;
+      //fNonLinearityParams[5] =  6.31604e+01;
+      //fNonLinearityParams[6] =  1.01286;
+      energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5])))) * (0.964 + exp(-3.132-0.435*energy*2.0));
+      
+      break;
+    }
+      
+    case kPi0MCv5:
+    {
+      //Based on comparing MC truth information to the reconstructed energy of clusters. 
+      //fNonLinearityParams[0] =  1.0;
+      //fNonLinearityParams[1] =  6.64778e-02;
+      //fNonLinearityParams[2] =  1.570;
+      //fNonLinearityParams[3] =  9.67998e-02;
+      //fNonLinearityParams[4] =  2.19381e+02;
+      //fNonLinearityParams[5] =  6.31604e+01;
+      //fNonLinearityParams[6] =  1.01286;
+      energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
+      
+      break;
+    }
+      
     case kNoCorrection:
       AliDebug(2,"No correction on the energy\n");
       break;
@@ -782,6 +814,27 @@ void AliEMCALRecoUtils::InitNonLinearityParam()
     fNonLinearityParams[5] =  31.5028;
     fNonLinearityParams[6] =  0.968;
   }
+
+  if (fNonLinearityFunction == kSDMv5) {
+    fNonLinearityParams[0] =  1.0;
+    fNonLinearityParams[1] =  6.64778e-02;
+    fNonLinearityParams[2] =  1.570;
+    fNonLinearityParams[3] =  9.67998e-02;
+    fNonLinearityParams[4] =  2.19381e+02;
+    fNonLinearityParams[5] =  6.31604e+01;
+    fNonLinearityParams[6] =  1.01286;
+  }
+
+  if (fNonLinearityFunction == kPi0MCv5) {
+    fNonLinearityParams[0] =  1.0;
+    fNonLinearityParams[1] =  6.64778e-02;
+    fNonLinearityParams[2] =  1.570;
+    fNonLinearityParams[3] =  9.67998e-02;
+    fNonLinearityParams[4] =  2.19381e+02;
+    fNonLinearityParams[5] =  6.31604e+01;
+    fNonLinearityParams[6] =  1.01286;
+  }
+
 }
 
 //_________________________________________________________
@@ -927,7 +980,7 @@ void AliEMCALRecoUtils::InitParameters()
   
   fExoticCellFraction     = 0.97;
   fExoticCellDiffTime     = 1e6;
-  fExoticCellMinAmplitude = 0.5;
+  fExoticCellMinAmplitude = 4.0;
   
   fAODFilterMask    = 128;
   fAODHybridTracks  = kFALSE;
@@ -1756,12 +1809,13 @@ void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
     //If the input event is AOD, the starting point for extrapolation is at vertex
     //AOD tracks are selected according to its filterbit.
     else if (aodevent) {
-      aodTrack = aodevent->GetTrack(itr);
+      aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(itr));
+      if(!aodTrack) AliFatal("Not a standard AOD");
       if (!aodTrack) continue;
             
       if (fAODTPCOnlyTracks) { // Match with TPC only tracks, default from May 2013, before filter bit 32
         //printf("Match with TPC only tracks, accept? %d, test bit 128 <%d> \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128));
-        if (!aodTrack->IsTPCOnly()) continue ;
+        if (!aodTrack->IsTPCConstrained()) continue ;
       } else if (fAODHybridTracks) { // Match with hybrid tracks
         //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal());
         if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
@@ -1948,13 +2002,15 @@ Int_t  AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackP
 
 //------------------------------------------------------------------------------------
 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliVTrack *track,
-                                                        Double_t emcalR, Double_t mass, Double_t step)
+                                                         Double_t emcalR, Double_t mass,
+                                                         Double_t step, Double_t minpt,
+                                                         Bool_t useMassForTracking)
 { 
   // Extrapolate track to EMCAL surface
 
   track->SetTrackPhiEtaPtOnEMCal(-999, -999, -999);
 
-  if (track->Pt()<0.350) //todo: discuss with Marta, everywhere else we use pT=0
+  if (track->Pt()<minpt)
     return kFALSE;
 
   Double_t phi = track->Phi()*TMath::RadToDeg();
@@ -1969,14 +2025,22 @@ Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliVTrack *track,
       return kFALSE;
   }
 
-  if (mass<0) {
+  // Select the mass hypothesis
+  if ( mass < 0 )
+  {
     Bool_t onlyTPC = kFALSE;
-    if (mass==-99)
-      onlyTPC=kTRUE;
+    if ( mass == -99 ) onlyTPC=kTRUE;
+    
     if (esdt)
-      mass = esdt->GetMass(onlyTPC);
-    else 
-      mass = aodt->M();
+    {
+      if ( useMassForTracking ) mass = esdt->GetMassForTracking();
+      else                      mass = esdt->GetMass(onlyTPC);
+    }
+    else
+    {
+      if ( useMassForTracking ) mass = aodt->GetMassForTracking();
+      else                      mass = aodt->M();
+    }
   }
 
   AliExternalTrackParam *trackParam = 0;
@@ -2487,8 +2551,13 @@ void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
     for (Int_t iTrk=0; iTrk<nTracks; ++iTrk) 
     {
       AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
-      if (iTrk == matchTrackIndex) continue;
-      if (track->GetEMCALcluster() == iClus) {
+      
+      if( !track ) continue;
+      
+      if ( iTrk == matchTrackIndex ) continue;
+      
+      if ( track->GetEMCALcluster() == iClus )
+      {
         arrayTrackMatched[nMatched] = iTrk;
         ++nMatched;
       }