]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGDQ/dielectron/AliDielectronVarManager.h
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
-sync with gsi repository
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronVarManager.h
diff --git a/PWGDQ/dielectron/AliDielectronVarManager.h b/PWGDQ/dielectron/AliDielectronVarManager.h
index 0ca74e81e7e9214f2489727bfb6e067add26e40b..78a35dfee48e448713201fcf0a53a3026a9f5c69 100644 (file)
--- a/PWGDQ/dielectron/AliDielectronVarManager.h
+++ b/PWGDQ/dielectron/AliDielectronVarManager.h
@@ -24,6 +24,7 @@
 #include <TProfile.h>
 #include <TProfile2D.h>
 #include <TH3D.h>
+#include <THnBase.h>
 #include <TFile.h>
 #include <TDatabasePDG.h>
 #include <TKey.h>
@@ -161,6 +162,7 @@ public:
     kTOFPIDBit,              // TOF PID bit (1:set, 0:TOF not available)a
     kTOFmismProb,               // and mismatchPorbability as explain in TOF-twiki
        
+    kTPCnSigmaEleRaw,        // raw number of sigmas to the dE/dx electron line in the TPC
     kTPCnSigmaEle,           // number of sigmas to the dE/dx electron line in the TPC
     kTPCnSigmaPio,           // number of sigmas to the dE/dx pion line in the TPC
     kTPCnSigmaMuo,           // number of sigmas to the dE/dx muon line in the TPC
@@ -181,6 +183,8 @@ public:
     kEMCALM20,               // M20 showershape parameter
     kEMCALDispersion,        // Dispersion paramter
     
+    kLegEff,                 // single electron efficiency
+    kOneOverLegEff,          // 1 / single electron efficiency (correction factor)
     kV0Index0,               // v0 index 0
     kKinkIndex0,             // kink index 0
       
@@ -245,6 +249,8 @@ public:
     kTRDpidEffPair,          // TRD pid efficieny from conversion electrons
     kMomAsymDau1,            // momentum fraction of daughter1
     kMomAsymDau2,            // momentum fraction of daughter2
+    kPairEff,                // pair efficiency
+    kOneOverPairEff,         // 1 / pair efficiency (correction factor)
     kPairMax,                 //
   // Event specific variables
     kXvPrim=kPairMax,        // prim vertex
@@ -396,6 +402,7 @@ public:
   static void InitAODpidUtil(Int_t type=0);
   static void InitEstimatorAvg(const Char_t* filename);
   static void InitTRDpidEffHistograms(const Char_t* filename);
+  static void InitEffMap(const Char_t* filename);
   static void SetVZEROCalibrationFile(const Char_t* filename) {fgVZEROCalibrationFile = filename;}
   
   static void SetVZERORecenteringFile(const Char_t* filename) {fgVZERORecenteringFile = filename;}
@@ -410,6 +417,7 @@ public:
 
   static TProfile* GetEstimatorHistogram(Int_t period, Int_t type) {return fgMultEstimatorAvg[period][type];}
   static Double_t GetTRDpidEfficiency(Int_t runNo, Double_t centrality, Double_t eta, Double_t trdPhi, Double_t pout, Double_t& effErr);
+  static Double_t GetSingleLegEff(Double_t * const values);
 
   static const AliKFVertex* GetKFVertex() {return fgKFVertex;}
   
@@ -451,6 +459,7 @@ private:
   static TProfile        *fgMultEstimatorAvg[4][9];  // multiplicity estimator averages (4 periods x 9 estimators)
   static Double_t         fgTRDpidEffCentRanges[10][4];   // centrality ranges for the TRD pid efficiency histograms
   static TH3D            *fgTRDpidEff[10][4];   // TRD pid efficiencies from conversion electrons
+  static THnBase         *fgEffMap;             // single electron efficiencies
   static TString          fgVZEROCalibrationFile;  // file with VZERO channel-by-channel calibrations
   static TString          fgVZERORecenteringFile;  // file with VZERO Q-vector averages needed for event plane recentering
   static TProfile2D      *fgVZEROCalib[64];           // 1 histogram per VZERO channel
@@ -699,7 +708,9 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle
   // nsigma to Electron band
   // TODO: for the moment we set the bethe bloch parameters manually
   //       this should be changed in future!
+  values[AliDielectronVarManager::kTPCnSigmaEleRaw]=fgPIDResponse->NumberOfSigmasTPC(particle,AliPID::kElectron);
   values[AliDielectronVarManager::kTPCnSigmaEle]=(fgPIDResponse->NumberOfSigmasTPC(particle,AliPID::kElectron)-AliDielectronPID::GetCorrVal()-AliDielectronPID::GetCntrdCorr(particle)) / AliDielectronPID::GetWdthCorr(particle);
+
   values[AliDielectronVarManager::kTPCnSigmaPio]=fgPIDResponse->NumberOfSigmasTPC(particle,AliPID::kPion);
   values[AliDielectronVarManager::kTPCnSigmaMuo]=fgPIDResponse->NumberOfSigmasTPC(particle,AliPID::kMuon);
   values[AliDielectronVarManager::kTPCnSigmaKao]=fgPIDResponse->NumberOfSigmasTPC(particle,AliPID::kKaon);
@@ -728,8 +739,9 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle
   values[AliDielectronVarManager::kEMCALM02]        = showershape[1];
   values[AliDielectronVarManager::kEMCALM20]        = showershape[2];
   values[AliDielectronVarManager::kEMCALDispersion] = showershape[3];
-  
-  
+
+  values[AliDielectronVarManager::kLegEff]        = GetSingleLegEff(values);
+  values[AliDielectronVarManager::kOneOverLegEff] = (values[AliDielectronVarManager::kLegEff]>0.0 ? 1./values[AliDielectronVarManager::kLegEff] : 0.0);
   //restore TPC signal if it was changed
   if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN());
 }
@@ -813,6 +825,7 @@ inline void AliDielectronVarManager::FillVarAODTrack(const AliAODTrack *particle
   values[AliDielectronVarManager::kTOFsignal]=0;
   //values[AliDielectronVarManager::kTOFbeta]=0;
 
+  values[AliDielectronVarManager::kTPCnSigmaEleRaw]=0;
   values[AliDielectronVarManager::kTPCnSigmaEle]=0;
   values[AliDielectronVarManager::kTPCnSigmaPio]=0;
   values[AliDielectronVarManager::kTPCnSigmaMuo]=0;
@@ -853,7 +866,7 @@ inline void AliDielectronVarManager::FillVarAODTrack(const AliAODTrack *particle
     
     values[AliDielectronVarManager::kPIn]=mom;
     values[AliDielectronVarManager::kTPCsignal]=pid->GetTPCsignal();
-
+    values[AliDielectronVarManager::kTPCnSigmaEleRaw]=fgPIDResponse->NumberOfSigmasTPC(particle,AliPID::kElectron);
     values[AliDielectronVarManager::kTPCnSigmaEle]=tpcNsigmaEle;
     values[AliDielectronVarManager::kTPCnSigmaPio]=tpcNsigmaPio;
     values[AliDielectronVarManager::kTPCnSigmaMuo]=tpcNsigmaMuo;
@@ -948,6 +961,9 @@ inline void AliDielectronVarManager::FillVarAODTrack(const AliAODTrack *particle
   } //if(mc->HasMC())
   
   values[AliDielectronVarManager::kTOFPIDBit]=(particle->GetStatus()&AliESDtrack::kTOFpid? 1: 0);
+  values[AliDielectronVarManager::kLegEff] = GetSingleLegEff(values);
+  values[AliDielectronVarManager::kOneOverLegEff] = (values[AliDielectronVarManager::kLegEff]>0.0 ? 1./values[AliDielectronVarManager::kLegEff] : 0.0);
+
 }
 
 inline void AliDielectronVarManager::FillVarMCParticle(const AliMCParticle *particle, Double_t * const values)
@@ -983,6 +999,7 @@ inline void AliDielectronVarManager::FillVarMCParticle(const AliMCParticle *part
   values[AliDielectronVarManager::kTPCsignal]     = 0;
   values[AliDielectronVarManager::kTOFsignal]     = 0;
   values[AliDielectronVarManager::kTOFbeta]       = 0;
+  values[AliDielectronVarManager::kTPCnSigmaEleRaw]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaEle]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaPio]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaMuo]  = 0;
@@ -1051,6 +1068,7 @@ inline void AliDielectronVarManager::FillVarMCParticle2(const AliVParticle *p1,
   values[AliDielectronVarManager::kPIn]           = 0;
   values[AliDielectronVarManager::kYsignedIn]     = 0;
   values[AliDielectronVarManager::kTPCsignal]     = 0;
+  values[AliDielectronVarManager::kTPCnSigmaEleRaw]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaEle]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaPio]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaMuo]  = 0;
@@ -1151,6 +1169,7 @@ inline void AliDielectronVarManager::FillVarAODMCParticle(const AliAODMCParticle
   values[AliDielectronVarManager::kPIn]           = 0;
   values[AliDielectronVarManager::kYsignedIn]     = 0;
   values[AliDielectronVarManager::kTPCsignal]     = 0;
+  values[AliDielectronVarManager::kTPCnSigmaEleRaw]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaEle]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaPio]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaMuo]  = 0;
@@ -1308,7 +1327,7 @@ inline void AliDielectronVarManager::FillVarDielectronPair(const AliDielectronPa
        values[AliDielectronVarManager::kOpeningAngle] =  lv1.Angle(lv2.Vect());
 
        values[AliDielectronVarManager::kOneOverPt] = (values[AliDielectronVarManager::kPt]>0. ? 1./values[AliDielectronVarManager::kPt] : -9999.);
-       values[AliDielectronVarManager::kPhi]       = (lv1+lv2).Phi();
+       values[AliDielectronVarManager::kPhi]       = (lv1+lv2).Phi()+TMath::Pi(); // change interval to [0,+2pi]
        values[AliDielectronVarManager::kEta]       = (lv1+lv2).Eta();
 
        values[AliDielectronVarManager::kY]       = (lv1+lv2).Rapidity();
@@ -1408,7 +1427,7 @@ inline void AliDielectronVarManager::FillVarDielectronPair(const AliDielectronPa
     values[AliDielectronVarManager::kDeltaPhiTPCrpH2] += TMath::TwoPi(); 
 
   //angle between ee plane and Mag/Reaction plane
-  values[AliDielectronVarManager::kPairPlaneAngle] = pair->PairPlaneAngle(values[AliDielectronVarManager::kv0CrpH2]);
+  values[AliDielectronVarManager::kPairPlaneAngle] = pair->PairPlaneAngle();\r
   //ee plane vector
   Double_t RotPairx = 0;
   Double_t RotPairy = 0;
@@ -1468,13 +1487,23 @@ inline void AliDielectronVarManager::FillVarDielectronPair(const AliDielectronPa
   AliVParticle* leg1 = pair->GetFirstDaughter();
   AliVParticle* leg2 = pair->GetSecondDaughter();
   if (leg1)
-       values[AliDielectronVarManager::kMomAsymDau1] = (values[AliDielectronVarManager::kP] != 0)? leg1->P()  / values[AliDielectronVarManager::kP]: 0;
+    values[AliDielectronVarManager::kMomAsymDau1] = (values[AliDielectronVarManager::kP] != 0)? leg1->P()  / values[AliDielectronVarManager::kP]: 0;
   else 
-       values[AliDielectronVarManager::kMomAsymDau1] = -9999.;
+    values[AliDielectronVarManager::kMomAsymDau1] = -9999.;
   if (leg2)
-       values[AliDielectronVarManager::kMomAsymDau2] = (values[AliDielectronVarManager::kP] != 0)? leg2->P()  / values[AliDielectronVarManager::kP]: 0;
+    values[AliDielectronVarManager::kMomAsymDau2] = (values[AliDielectronVarManager::kP] != 0)? leg2->P()  / values[AliDielectronVarManager::kP]: 0;
   else 
-       values[AliDielectronVarManager::kMomAsymDau2] = -9999.;
+    values[AliDielectronVarManager::kMomAsymDau2] = -9999.;
+
+  Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
+  Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
+  if (leg1 && leg2 && fgEffMap) {
+    Fill(leg1, valuesLeg1);
+    Fill(leg2, valuesLeg2);
+    values[AliDielectronVarManager::kPairEff] = valuesLeg1[AliDielectronVarManager::kLegEff] *valuesLeg2[AliDielectronVarManager::kLegEff];
+    values[AliDielectronVarManager::kOneOverPairEff] = (values[AliDielectronVarManager::kPairEff]>0.0 ? 1./values[AliDielectronVarManager::kPairEff] : 0.0);
+
+  }
 }
 
 inline void AliDielectronVarManager::FillVarKFParticle(const AliKFParticle *particle, Double_t * const values)
@@ -1527,6 +1556,7 @@ inline void AliDielectronVarManager::FillVarKFParticle(const AliKFParticle *part
   values[AliDielectronVarManager::kTPCsignal]     = 0;
   values[AliDielectronVarManager::kTOFsignal]     = 0;
   values[AliDielectronVarManager::kTOFbeta]       = 0;
+  values[AliDielectronVarManager::kTPCnSigmaEleRaw]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaEle]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaPio]  = 0;
   values[AliDielectronVarManager::kTPCnSigmaMuo]  = 0;
@@ -2109,6 +2139,42 @@ inline void AliDielectronVarManager::InitTRDpidEffHistograms(const Char_t* filen
   }
 }
 
+inline void AliDielectronVarManager::InitEffMap(const Char_t* filename) {
+  //
+  // init an efficiency object for on-the-fly correction calculations
+  //
+  fgEffMap=0x0;
+  TFile* file=TFile::Open(filename);
+  THnBase *hGen = (THnBase*) file->Get("hGenerated");
+  THnBase *hFnd = (THnBase*) file->Get("hFound");
+  if(!hFnd || !hGen) return;
+
+  fgEffMap  = (THnBase*) hFnd->Clone("effMap");
+  fgEffMap->Reset();
+  fgEffMap->Sumw2();
+  fgEffMap->Divide(hFnd, hGen);//, 1, 1, "");  //assume uncorrelated err, otherwise give option "B"
+  if(fgEffMap) printf("[I] AliDielectronVarManager::InitEffMap efficiency maps loaded! \n");
+}
+
+inline Double_t AliDielectronVarManager::GetSingleLegEff(Double_t * const values) {
+  //
+  // get the single leg efficiency for a given particle
+  //
+  if(!fgEffMap) return -1.;
+
+  Int_t dim=fgEffMap->GetNdimensions();
+  Int_t idx[dim];
+  for(Int_t idim=0; idim<dim; idim++) {
+    UInt_t var = GetValueType(fgEffMap->GetAxis(idim)->GetName());
+    idx[idim] = fgEffMap->GetAxis(idim)->FindBin(values[var]);
+    if(idx[idim] < 0 || idx[idim]>fgEffMap->GetAxis(idim)->GetNbins()) 
+      printf(" [E] AliDielectronVarManager::GetSingleLegEff values %f for %s not found in axis range \n",values[var],fgEffMap->GetAxis(idim)->GetName());
+    //    printf(" (%d,%f,%s) \t",idx[idim],values[var],fgEffMap->GetAxis(idim)->GetName());
+  }
+  //  printf(" bin content %f+-%f \n",fgEffMap->GetBinContent(idx), fgEffMap->GetBinError(idx));
+  return (fgEffMap->GetBinContent(idx));
+}
+
 
 inline void AliDielectronVarManager::InitVZEROCalibrationHistograms(Int_t runNo) {
   //