Add flow tasks to the compilation, update others
authorrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 6 Mar 2012 10:05:06 +0000 (10:05 +0000)
committerrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 6 Mar 2012 10:05:06 +0000 (10:05 +0000)
26 files changed:
PWGHF/CMakelibPWGHFhfe.pkg
PWGHF/PWGHFhfeLinkDef.h
PWGHF/hfe/AliAnalysisTaskHFE.cxx
PWGHF/hfe/AliAnalysisTaskHFE.h
PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx
PWGHF/hfe/AliAnalysisTaskHFEFlow.h
PWGHF/hfe/AliHFEVZEROEventPlane.cxx [new file with mode: 0644]
PWGHF/hfe/AliHFEVZEROEventPlane.h [new file with mode: 0644]
PWGHF/hfe/AliHFEcollection.cxx
PWGHF/hfe/AliHFEcollection.h
PWGHF/hfe/AliHFEdebugTreeTask.cxx
PWGHF/hfe/AliHFEdebugTreeTask.h
PWGHF/hfe/AliHFEemcalPIDqa.cxx
PWGHF/hfe/AliHFEextraCuts.h
PWGHF/hfe/AliHFEmcQA.cxx
PWGHF/hfe/AliHFEmcQA.h
PWGHF/hfe/AliHFEpidEMCAL.cxx
PWGHF/hfe/AliHFEpidQA.cxx
PWGHF/hfe/AliHFEpidQA.h
PWGHF/hfe/AliHFEsignalCuts.cxx
PWGHF/hfe/AliHFEspectrum.cxx
PWGHF/hfe/AliHFEspectrum.h
PWGHF/hfe/AliHFEtools.cxx
PWGHF/hfe/AliHFEtools.h
PWGHF/hfe/AliHFEtrdPIDqa.cxx
PWGHF/hfe/AliHFEtrdPIDqa.h

index a8c1964..c6f8284 100644 (file)
@@ -79,10 +79,12 @@ set (SRCS
     hfe/AliHFEpidObject.cxx
     hfe/AliAnalysisTaskElecHadronCorrel.cxx
     hfe/AliHFEdebugTreeTask.cxx
+    hfe/AliHFEVZEROEventPlane.cxx
+    hfe/AliAnalysisTaskHFEFlow.cxx
     )
 
 string (REPLACE ".cxx" ".h" HDRS "${SRCS}")
 
 set ( DHDR  PWGHFhfeLinkDef.h)
 
-set ( EINCLUDE  PWGHF/hfe TPC CORRFW STEER/STEERBase)
+set ( EINCLUDE  PWGHF/hfe TPC CORRFW STEER/STEERBase PWG/FLOW/Base PWG/FLOW/Tasks)
index 5cc2f59..7bdd0f0 100644 (file)
@@ -70,4 +70,7 @@
 
 #pragma link C++ class  AliHFEdebugTreeTask+;
 
+#pragma link C++ class  AliHFEVZEROEventPlane+;
+#pragma link C++ class  AliAnalysisTaskHFEFlow+;
+
 #endif
index c5386ab..c9e6eb0 100644 (file)
@@ -13,7 +13,6 @@
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
 //
-//
 // The analysis task:
 // Filling an AliCFContainer with the quantities pt, eta and phi
 // for tracks which survivied the particle cuts (MC resp. ESD tracks)
@@ -432,6 +431,9 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
     fHistMCQA->SetOwner();
     if(IsPbPb()) fMCQA->SetPbPb();
     if(fisppMultiBin) fMCQA->SetPPMultiBin();
+    if(TestBit(kTreeStream)){
+      fMCQA->EnableDebugStreamer();
+    }
     fMCQA->CreatDefaultHistograms(fHistMCQA);
     fMCQA->SetBackgroundWeightFactor(fElecBackgroundFactor[0][0][0],fBinLimit);
     fQA->Add(fHistMCQA);
@@ -711,19 +713,9 @@ void AliAnalysisTaskHFE::ProcessMC(){
           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
-          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0); // no accept cut
-          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
-          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 0); // no accept cut
-          if (TMath::Abs(mcpart->Eta()) < 0.9) {
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
-          }
-          if (TMath::Abs(AliHFEtools::GetRapidity(mcpart)) < 0.5) {
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-            fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
-          }
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG); // no accept cut
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG); // no accept cut
+          fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG); // no accept cut
         }
         //fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
         //fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
@@ -830,6 +822,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
   memset(container, 0, sizeof(Double_t) * 10);
   // container for the output THnSparse
   Double_t dataE[6]; // [pT, eta, Phi, type, 'C' or 'B']
+  Double_t dataDca[7]; // [source, pT, dca, dcaSig, centrality]
   Int_t nElectronCandidates = 0;
   AliESDtrack *track = NULL, *htrack = NULL;
   AliMCParticle *mctrack = NULL;
@@ -938,10 +931,22 @@ void AliAnalysisTaskHFE::ProcessESD(){
     if(HasMCData()){
       FillProductionVertex(track);
 
-      if(fMCQA){
+      if(fMCQA && signal){
         fMCQA->SetCentrality(fCentralityF);
         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){
          Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.};
+         Double_t hfeimpactRtmp=0., hfeimpactnsigmaRtmp=0.;
+         fExtraCuts->GetHFEImpactParameters(track, hfeimpactRtmp, hfeimpactnsigmaRtmp);
+         UChar_t itsPixel = track->GetITSClusterMap();
+         Double_t ilyrhit=0, ilyrstat=0;
+         for(Int_t ilyr=0; ilyr<6; ilyr++){
+           if(TESTBIT(itsPixel, ilyr)) ilyrhit += TMath::Power(2,ilyr);
+           if(fExtraCuts->CheckITSstatus(fExtraCuts->GetITSstatus(track,ilyr))) ilyrstat += TMath::Power(2,ilyr);
+         }
+         fMCQA->SetITSInfo(ilyrhit,ilyrstat);
+         fMCQA->SetHFEImpactParameters(hfeimpactRtmp, hfeimpactnsigmaRtmp);
+         fMCQA->SetTrkKine(track->Pt(),track->Eta(), track->Phi());
+         fMCQA->SetContainerStep(3);
          for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
            weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE 
            if(!fisNonHFEsystematics)break;   
@@ -958,8 +963,12 @@ void AliAnalysisTaskHFE::ProcessESD(){
              if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
              if(elecSource == iSource){
                for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
-                 if(weightElecBgV0[iLevel]>0){ fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);}
-                 else if(weightElecBgV0[iLevel]<0){ fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);}
+                 if(weightElecBgV0[iLevel]>0){ 
+                   fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);
+                 } 
+                 else if(weightElecBgV0[iLevel]<0){ 
+                   fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);
+                 }
                }
                break;
              }
@@ -968,15 +977,22 @@ void AliAnalysisTaskHFE::ProcessESD(){
            }
          }
          //else{
-           if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]);
-           else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]);
+           if(weightElecBgV0[0]>0) {
+            fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]);
+            fVarManager->FillContainer(fContainer, "conversionElecs", 4, kTRUE, weightElecBgV0[0]);
+          }
+           else if(weightElecBgV0[0]<0) {
+            fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]);
+            fVarManager->FillContainer(fContainer, "mesonElecs", 4, kTRUE, -1*weightElecBgV0[0]);
+          }
            //}
         }
       }
     }
 
     if(TMath::Abs(track->Eta()) < 0.5){
-      fQACollection->Fill("TPCdEdxBeforePID", track->P(), track->GetTPCsignal());
+      if(track->GetInnerParam())
+        fQACollection->Fill("TPCdEdxBeforePID", track->GetInnerParam()->P(), track->GetTPCsignal());
       fQACollection->Fill("TPCnSigmaBeforePID", track->P(), fInputHandler->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron));
     }
 
@@ -1002,11 +1018,11 @@ void AliAnalysisTaskHFE::ProcessESD(){
            if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
 
             Double_t itsChi2[7] = {track->Pt(),track->Eta(), track->Phi(),
-                                  fCentralityF,track->GetTPCsignalN(), sharebit, itschi2percluster};
+                                  fCentralityF,track->GetTPCsignalN(), sharebit,itschi2percluster};
             fQACollection->Fill("fChi2perITScluster", itsChi2);
     }
     else{
-      
+
       Double_t itschi2percluster = 0.0;
       Double_t itsnbcls = static_cast<Double_t>(track->GetNcls(0));
       if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
@@ -1024,9 +1040,9 @@ void AliAnalysisTaskHFE::ProcessESD(){
         Int_t glabel=TMath::Abs(mctrack->GetMother());
         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
           if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==321)
-            fQACollection->Fill("Ke3Kecorr",mctrackmother->Pt(),mctrack->Pt());
+            fQACollection->Fill("Ke3Kecorr",mctrack->Pt(),mctrackmother->Pt());
           else if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==130)
-            fQACollection->Fill("Ke3K0Lecorr",mctrackmother->Pt(),mctrack->Pt());
+            fQACollection->Fill("Ke3K0Lecorr",mctrack->Pt(),mctrackmother->Pt());
         }
       }
     }
@@ -1107,42 +1123,41 @@ void AliAnalysisTaskHFE::ProcessESD(){
     } // end of electron background analysis
 
 
-
+    Int_t sourceDca =-1;
     if (GetPlugin(kDEstep)) { 
       Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.,};
       Int_t elecSource = 0;
       // minjung for IP QA(temporary ~ 2weeks)
       Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
       fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
-      fQACollection->Fill("dataDca",track->Pt(),hfeimpactR);
-      fQACollection->Fill("dataDcaSig",track->Pt(),hfeimpactnsigmaR);
-      fQACollection->Fill("dataDcaSigDca",hfeimpactR,hfeimpactnsigmaR);
-      if(HasMCData()){
+      sourceDca=0;
+      if(HasMCData())
+      {
         // minjung for IP QA(temporary ~ 2weeks)
-        if(fSignalCuts->IsCharmElectron(track)){
-          fQACollection->Fill("charmDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("charmDcaSig",track->Pt(),hfeimpactnsigmaR);
+         if(fSignalCuts->IsCharmElectron(track)){
+             sourceDca=1;
         }
-        else if(fSignalCuts->IsBeautyElectron(track)){
-          fQACollection->Fill("beautyDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("beautyDcaSig",track->Pt(),hfeimpactnsigmaR);
+         else if(fSignalCuts->IsBeautyElectron(track)){
+             sourceDca=2;
         }
-        else if(fSignalCuts->IsGammaElectron(track)){
-          fQACollection->Fill("conversionDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("conversionDcaSig",track->Pt(),hfeimpactnsigmaR);
-          fQACollection->Fill("conversionDcaSigDca",hfeimpactR,hfeimpactnsigmaR);
+         else if(fSignalCuts->IsGammaElectron(track)){
+             sourceDca=3;
         }
-        else if(fSignalCuts->IsNonHFElectron(track)){
-          fQACollection->Fill("nonhfeDca",track->Pt(),hfeimpactR);
-          fQACollection->Fill("nonhfeDcaSig",track->Pt(),hfeimpactnsigmaR);
+         else if(fSignalCuts->IsNonHFElectron(track)){
+             sourceDca=4;
         }
-
-        if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
+         else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
+             sourceDca=5;
           fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
           fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
         }
-        if(fMCQA) {
+         else {
+             sourceDca=6;
+       }
+
+        if(fMCQA && signal) {
           
+          fMCQA->SetContainerStep(0);
           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE 
             if(!fisNonHFEsystematics)break;        
@@ -1168,19 +1183,45 @@ void AliAnalysisTaskHFE::ProcessESD(){
             }
           }
           //else{
-          if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]);
-          else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]);
+          if(weightElecBgV0[0]>0) {
+           fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]);
+           fVarManager->FillContainer(fContainer, "conversionElecs", 5, kTRUE, weightElecBgV0[0]);
+         }
+          else if(weightElecBgV0[0]<0) {
+           fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]);
+           fVarManager->FillContainer(fContainer, "mesonElecs", 5, kTRUE, -1*weightElecBgV0[0]);
+         }  
           //}
           if(bTagged){ // bg estimation for the secondary vertex tagged signals
             if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 2, kFALSE, weightElecBgV0[0]);
             else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 2, kFALSE, -1*weightElecBgV0[0]);
           }
         }
-      }
+      } // end of MC
+
+      dataDca[0]=sourceDca;
+      dataDca[1]=track->Pt();
+      dataDca[2]=hfeimpactR;
+      dataDca[3]=hfeimpactnsigmaR;
+      dataDca[4]=fCentralityF;
+      dataDca[5] = 49;
+      Double_t xr[3]={49,49,49};
+      if(HasMCData()) {
+        mctrack->XvYvZv(xr);
+        dataDca[5] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+      } 
+      dataDca[6] = v0pid;
+
+ //     printf("Entries dca: [%.3f|%.3f|%.3f|%f|%f]\n", dataDca[0],dataDca[1],dataDca[2],dataDca[3],dataDca[4]);
+      if (!HasMCData()) fQACollection->Fill("Dca", dataDca);
+      else if(signal) fQACollection->Fill("Dca", dataDca);
+
+
       // Fill Containers for impact parameter analysis
       if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
       if(HasMCData()){
-        if(fMCQA) {
+        if(fMCQA && signal) {
+          fMCQA->SetContainerStep(1);
           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE 
             if(!fisNonHFEsystematics)break;        
@@ -1205,8 +1246,14 @@ void AliAnalysisTaskHFE::ProcessESD(){
             }
           }
           // else{
-            if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]);
-            else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]);
+            if(weightElecBgV0[0]>0) {
+             fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]);
+             fVarManager->FillContainer(fContainer, "conversionElecs", 6, kTRUE, weightElecBgV0[0]);
+           }
+            else if(weightElecBgV0[0]<0) {
+             fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]);
+             fVarManager->FillContainer(fContainer, "mesonElecs", 6, kTRUE, -1*weightElecBgV0[0]);
+            }
             //}
         }
       }
@@ -1517,8 +1564,8 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
   fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1);
 
   if(HasMCData()){
-    fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",4);
-    fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",4);
+    fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",7);
+    fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",7);
     fContainer->Sumw2("conversionElecs");
     fContainer->Sumw2("mesonElecs");
    
@@ -1588,36 +1635,55 @@ void AliAnalysisTaskHFE::InitContaminationQA(){
   // 
   // Add QA for Impact Parameter cut
   //
-  const Double_t kPtbound[2] = {0.1, 20.};
-  Int_t iBin[1];
-  iBin[0] = 44; // bins in pt
-  fQACollection->CreateTH1F("hadronsBeforeIPcut", "Hadrons before IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("hadronsAfterIPcut", "Hadrons after IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", iBin[0], kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
-
-  fQACollection->CreateTH2F("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1);
-  fQACollection->CreateTH2F("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1);
-  fQACollection->CreateTH1F("Kptspectra", "Charged Kaons: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
-  fQACollection->CreateTH1F("K0Lptspectra", "K0L: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
-
-  const Double_t kDCAbound[2] = {-5., 5.};
-  const Double_t kDCAsigbound[2] = {-50., 50.};
-
-  fQACollection->CreateTH2F("dataDcaSig", "data dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("charmDcaSig", "charm dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("beautyDcaSig", "beauty dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("conversionDcaSig", "conversion dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("nonhfeDcaSig", "nonhfe dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0);
-
-  fQACollection->CreateTH2F("dataDca", "data dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("charmDca", "charm dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("beautyDca", "beauty dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("conversionDca", "conversion dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-  fQACollection->CreateTH2F("nonhfeDca", "nonhfe dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0);
-
-  fQACollection->CreateTH2F("dataDcaSigDca", "data dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0);
-  fQACollection->CreateTH2F("conversionDcaSigDca", "conversion dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0);
+
+  TObjArray *array = fVarManager->GetVariables();
+  Int_t nvars = array->GetEntriesFast();
+  for(Int_t v = 0; v < nvars; v++) {
+    AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
+    if(!variable) continue;
+    TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
+    if(!name.CompareTo("pt")) {
+      const Int_t nBinPt  = variable->GetNumberOfBins();
+      const Double_t *kPtRange = variable->GetBinning();
+
+      fQACollection->CreateTH1Farray("hadronsBeforeIPcut", "Hadrons before IP cut", nBinPt, kPtRange);
+      fQACollection->CreateTH1Farray("hadronsAfterIPcut", "Hadrons after IP cut", nBinPt, kPtRange);
+      fQACollection->CreateTH1Farray("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", nBinPt, kPtRange);
+      fQACollection->CreateTH1Farray("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", nBinPt, kPtRange);
+
+      fQACollection->CreateTH2Farray("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
+      fQACollection->CreateTH2Farray("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
+      fQACollection->CreateTH1Farray("Kptspectra", "Charged Kaons: MC p_{t} ", nBinPt, kPtRange);
+      fQACollection->CreateTH1Farray("K0Lptspectra", "K0L: MC p_{t} ", nBinPt, kPtRange);
+
+      const Double_t kDCAbound[2] = {-5., 5.};
+      const Double_t kDCAsigbound[2] = {-50., 50.};
+
+      const Int_t nDimDca=7;
+      const Int_t nBinDca[nDimDca] = { 8, nBinPt, 2000, 2000, 12, 500, 6};
+      Double_t minimaDca[nDimDca]  = { -1., 0., kDCAbound[0], kDCAsigbound[0], -1., 0, -1};
+      Double_t maximaDca[nDimDca]  = { 7., 20., kDCAbound[1], kDCAsigbound[1], 11., 50, 5};
+
+      Double_t *sourceBins = AliHFEtools::MakeLinearBinning(nBinDca[0], minimaDca[0], maximaDca[0]);
+      Double_t *dcaBins = AliHFEtools::MakeLinearBinning(nBinDca[2], minimaDca[2], maximaDca[2]);
+      Double_t *dcaSigBins = AliHFEtools::MakeLinearBinning(nBinDca[3], minimaDca[3], maximaDca[3]);
+      Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBinDca[4], minimaDca[4], maximaDca[4]);
+      Double_t *eProdRBins = AliHFEtools::MakeLinearBinning(nBinDca[5], minimaDca[5], maximaDca[5]);
+      Double_t *v0PIDBins = AliHFEtools::MakeLinearBinning(nBinDca[6], minimaDca[6], maximaDca[6]);
+
+      fQACollection->CreateTHnSparseNoLimits("Dca", "Dca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; dcasig; centrality bin; eProdR; v0pid", nDimDca, nBinDca);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(0, sourceBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(1, kPtRange);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(2, dcaBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(3, dcaSigBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(4, centralityBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(5, eProdRBins);
+      ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(6, v0PIDBins);
+
+      break;
+    }  
+  }
+
 }
 
 //____________________________________________________________
index 9021a30..e10ee95 100644 (file)
@@ -133,13 +133,15 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{
     Bool_t ReadCentrality();
     void RejectionPileUpVertexRangeEventCut();  
     void SelectSpecialTrigger(const Char_t *trgclust, Int_t runMin = 0, Int_t runMax = 999999999); 
+    void SetDebugStreaming() {SetBit(kTreeStream);};
     
   private:
     enum{
       kHasMCdata = BIT(19),
       kAODanalysis = BIT(20),
       kBeamType = BIT(21),
-      kBackgroundInitialized = BIT(22)
+      kBackgroundInitialized = BIT(22),
+      kTreeStream = BIT(23)
     };
 
     Bool_t FillProductionVertex(const AliVParticle * const track) const;
index 437424d..e11df38 100644 (file)
 #include "AliFlowVector.h"\r
 #include "AliFlowCommonConstants.h"\r
 \r
-\r
 #include "AliHFEcuts.h"\r
 #include "AliHFEpid.h"\r
 #include "AliHFEpidQAmanager.h"\r
 #include "AliHFEtools.h"\r
-\r
-\r
+#include "AliHFEVZEROEventPlane.h"\r
 \r
 #include "AliCentrality.h"\r
 #include "AliEventplane.h"\r
@@ -62,7 +60,7 @@
 //____________________________________________________________________\r
 AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :\r
   AliAnalysisTaskSE(),\r
-  fListHist(), \r
+  fListHist(0x0), \r
   fVZEROEventPlane(kFALSE),\r
   fVZEROEventPlaneA(kFALSE),\r
   fVZEROEventPlaneC(kFALSE),\r
@@ -83,6 +81,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :
   fPrecisionPhi(0.001),\r
   fUseMCReactionPlane(kFALSE),\r
   fMCPID(kFALSE),\r
+  fNoPID(kFALSE),\r
   fDebugLevel(0),\r
   fcutsRP(0),\r
   fcutsPOI(0),\r
@@ -90,6 +89,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :
   fPID(0),\r
   fPIDqa(0),\r
   fflowEvent(0x0),\r
+  fHFEVZEROEventPlane(0x0),\r
   fHistEV(0),\r
   fEventPlane(0x0),\r
   fEventPlaneaftersubtraction(0x0),\r
@@ -100,10 +100,12 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :
   fSin2phiep(0x0),\r
   fSin2phiephiep(0x0),\r
   fCosResabc(0x0),\r
+  fSinResabc(0x0),\r
   fProfileCosResab(0x0),\r
   fProfileCosResac(0x0),\r
   fProfileCosResbc(0x0),\r
   fCosRes(0x0),\r
+  fSinRes(0x0),\r
   fProfileCosRes(0x0),\r
   fDeltaPhiMaps(0x0),\r
   fCosPhiMaps(0x0),\r
@@ -119,7 +121,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :
 //______________________________________________________________________________\r
 AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :\r
   AliAnalysisTaskSE(name),\r
-  fListHist(), \r
+  fListHist(0x0), \r
   fVZEROEventPlane(kFALSE),\r
   fVZEROEventPlaneA(kFALSE),\r
   fVZEROEventPlaneC(kFALSE),\r
@@ -140,6 +142,7 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :
   fPrecisionPhi(0.001),\r
   fUseMCReactionPlane(kFALSE),\r
   fMCPID(kFALSE),\r
+  fNoPID(kFALSE),\r
   fDebugLevel(0),\r
   fcutsRP(0),\r
   fcutsPOI(0),\r
@@ -147,6 +150,7 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :
   fPID(0),\r
   fPIDqa(0),\r
   fflowEvent(0x0),\r
+  fHFEVZEROEventPlane(0x0),\r
   fHistEV(0),\r
   fEventPlane(0x0),\r
   fEventPlaneaftersubtraction(0x0),\r
@@ -157,10 +161,12 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :
   fSin2phiep(0x0),\r
   fSin2phiephiep(0x0),\r
   fCosResabc(0x0),\r
+  fSinResabc(0x0),\r
   fProfileCosResab(0x0),\r
   fProfileCosResac(0x0),\r
   fProfileCosResbc(0x0),\r
   fCosRes(0x0),\r
+  fSinRes(0x0),\r
   fProfileCosRes(0x0),\r
   fDeltaPhiMaps(0x0),\r
   fCosPhiMaps(0x0),\r
@@ -189,6 +195,147 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :
   }\r
   \r
 }\r
+//____________________________________________________________\r
+AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref):\r
+  AliAnalysisTaskSE(ref),\r
+  fListHist(0x0), \r
+  fVZEROEventPlane(ref.fVZEROEventPlane),\r
+  fVZEROEventPlaneA(ref.fVZEROEventPlaneA),\r
+  fVZEROEventPlaneC(ref.fVZEROEventPlaneC),\r
+  fSubEtaGapTPC(ref.fSubEtaGapTPC),\r
+  fEtaGap(ref.fEtaGap),\r
+  fNbBinsCentralityQCumulant(ref.fNbBinsCentralityQCumulant),\r
+  fNbBinsPtQCumulant(ref.fNbBinsPtQCumulant),\r
+  fMinPtQCumulant(ref.fMinPtQCumulant),\r
+  fMaxPtQCumulant(ref.fMaxPtQCumulant),\r
+  fAfterBurnerOn(ref.fAfterBurnerOn),\r
+  fNonFlowNumberOfTrackClones(ref.fNonFlowNumberOfTrackClones),\r
+  fV1(ref.fV1),\r
+  fV2(ref.fV2),\r
+  fV3(ref.fV3),\r
+  fV4(ref.fV4),\r
+  fV5(ref.fV5),\r
+  fMaxNumberOfIterations(ref.fMaxNumberOfIterations),\r
+  fPrecisionPhi(ref.fPrecisionPhi),\r
+  fUseMCReactionPlane(ref.fUseMCReactionPlane),\r
+  fMCPID(ref.fMCPID),\r
+  fNoPID(ref.fNoPID),\r
+  fDebugLevel(ref.fDebugLevel),\r
+  fcutsRP(0),\r
+  fcutsPOI(0),\r
+  fHFECuts(0),\r
+  fPID(0),\r
+  fPIDqa(0),\r
+  fflowEvent(0x0),\r
+  fHFEVZEROEventPlane(0x0),\r
+  fHistEV(0),\r
+  fEventPlane(0x0),\r
+  fEventPlaneaftersubtraction(0x0),\r
+  fCosSin2phiep(0x0),\r
+  fCos2phie(0x0),\r
+  fSin2phie(0x0),\r
+  fCos2phiep(0x0),\r
+  fSin2phiep(0x0),\r
+  fSin2phiephiep(0x0),\r
+  fCosResabc(0x0),\r
+  fSinResabc(0x0),\r
+  fProfileCosResab(0x0),\r
+  fProfileCosResac(0x0),\r
+  fProfileCosResbc(0x0),\r
+  fCosRes(0x0),\r
+  fSinRes(0x0),\r
+  fProfileCosRes(0x0),\r
+  fDeltaPhiMaps(0x0),\r
+  fCosPhiMaps(0x0),\r
+  fProfileCosPhiMaps(0x0)\r
+{\r
+  //\r
+  // Copy Constructor\r
+  //\r
+  ref.Copy(*this);\r
+}\r
+\r
+//____________________________________________________________\r
+AliAnalysisTaskHFEFlow &AliAnalysisTaskHFEFlow::operator=(const AliAnalysisTaskHFEFlow &ref){\r
+  //\r
+  // Assignment operator\r
+  //\r
+  if(this == &ref) \r
+    ref.Copy(*this);\r
+  return *this;\r
+}\r
+\r
+//____________________________________________________________\r
+void AliAnalysisTaskHFEFlow::Copy(TObject &o) const {\r
+  // \r
+  // Copy into object o\r
+  //\r
+  AliAnalysisTaskHFEFlow &target = dynamic_cast<AliAnalysisTaskHFEFlow &>(o);\r
+  target.fVZEROEventPlane = fVZEROEventPlane;\r
+  target.fVZEROEventPlaneA = fVZEROEventPlaneA;\r
+  target.fVZEROEventPlaneC = fVZEROEventPlaneC;\r
+  target.fSubEtaGapTPC = fSubEtaGapTPC;\r
+  target.fEtaGap = fEtaGap;\r
+  target.fNbBinsCentralityQCumulant = fNbBinsCentralityQCumulant;\r
+  target.fNbBinsPtQCumulant = fNbBinsPtQCumulant;\r
+  target.fMinPtQCumulant = fMinPtQCumulant;\r
+  target.fMaxPtQCumulant = fMaxPtQCumulant;\r
+  target.fAfterBurnerOn = fAfterBurnerOn;\r
+  target.fNonFlowNumberOfTrackClones = fNonFlowNumberOfTrackClones;\r
+  target.fV1 = fV1;\r
+  target.fV2 = fV2;\r
+  target.fV3 = fV3;\r
+  target.fV4 = fV4;\r
+  target.fV5 = fV5;\r
+  target.fMaxNumberOfIterations = fMaxNumberOfIterations;\r
+  target.fPrecisionPhi = fPrecisionPhi;\r
+  target.fUseMCReactionPlane = fUseMCReactionPlane;\r
+  target.fMCPID = fMCPID;\r
+  target.fNoPID = fNoPID;\r
+  target.fDebugLevel = fDebugLevel;\r
+  target.fcutsRP = fcutsRP;\r
+  target.fcutsPOI = fcutsPOI;\r
+  target.fHFECuts = fHFECuts;\r
+  target.fPID = fPID;\r
+  target.fPIDqa = fPIDqa;\r
+  target.fHFEVZEROEventPlane = fHFEVZEROEventPlane;\r
\r
+}\r
+//____________________________________________________________\r
+AliAnalysisTaskHFEFlow::~AliAnalysisTaskHFEFlow(){\r
+  //\r
+  // Destructor\r
+  //\r
+  if(fListHist) delete fListHist;\r
+  if(fcutsRP) delete fcutsRP;\r
+  if(fcutsPOI) delete fcutsPOI;\r
+  if(fHFECuts) delete fHFECuts;\r
+  if(fPID) delete fPID;\r
+  if(fPIDqa) delete fPIDqa;\r
+  if(fflowEvent) delete fflowEvent;\r
+  if(fHFEVZEROEventPlane) delete fHFEVZEROEventPlane;\r
+  if(fHistEV) delete fHistEV;\r
+  if(fEventPlane) delete fEventPlane;\r
+  if(fEventPlaneaftersubtraction) delete fEventPlaneaftersubtraction;\r
+  if(fCosSin2phiep) delete fCosSin2phiep;\r
+  if(fCos2phie) delete fCos2phie;\r
+  if(fSin2phie) delete fSin2phie;\r
+  if(fCos2phiep) delete fCos2phiep;\r
+  if(fSin2phiep) delete fSin2phiep;\r
+  if(fSin2phiephiep) delete fSin2phiephiep;\r
+  if(fCosResabc) delete fCosResabc;\r
+  if(fSinResabc) delete fSinResabc;\r
+  if(fProfileCosResab) delete fProfileCosResab;\r
+  if(fProfileCosResac) delete fProfileCosResac;\r
+  if(fProfileCosResbc) delete fProfileCosResbc;\r
+  if(fCosRes) delete fCosRes;\r
+  if(fSinRes) delete fSinRes;\r
+  if(fProfileCosRes) delete fProfileCosRes;\r
+  if(fDeltaPhiMaps) delete fDeltaPhiMaps;\r
+  if(fCosPhiMaps) delete fCosPhiMaps;\r
+  if(fProfileCosPhiMaps) delete fProfileCosPhiMaps;\r
+\r
+}\r
 //________________________________________________________________________\r
 void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()\r
 {\r
@@ -391,6 +538,15 @@ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()
   fCosResabc->SetBinEdges(3,binLimCMore);\r
   fCosResabc->Sumw2();\r
 \r
+  const Int_t nDimfbiss=4;\r
+  Int_t nBinfbiss[nDimfbiss] = {nBinsCos,nBinsCos,nBinsCos,nBinsC};\r
+  fSinResabc = new THnSparseF("SinRes_abc","SinRes_abc",nDimfbiss,nBinfbiss);\r
+  fSinResabc->SetBinEdges(0,binLimCos);\r
+  fSinResabc->SetBinEdges(1,binLimCos);\r
+  fSinResabc->SetBinEdges(2,binLimCos);\r
+  fSinResabc->SetBinEdges(3,binLimC);\r
+  fSinResabc->Sumw2();\r
+\r
   // Profile cosres centrality with 3 subevents\r
   fProfileCosResab = new TProfile("ProfileCosRes_a_b","ProfileCosRes_a_b",nBinsCMore,binLimCMore);\r
   fProfileCosResab->Sumw2();\r
@@ -407,6 +563,13 @@ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()
   fCosRes->SetBinEdges(1,binLimCMore);\r
   fCosRes->Sumw2();\r
 \r
+  const Int_t nDimff=2;\r
+  Int_t nBinff[nDimff] = {nBinsCos, nBinsC};\r
+  fSinRes = new THnSparseF("SinRes","SinRes",nDimff,nBinff);\r
+  fSinRes->SetBinEdges(0,binLimCos);\r
+  fSinRes->SetBinEdges(1,binLimC);\r
+  fSinRes->Sumw2();\r
+\r
   // Profile cosres centrality\r
   fProfileCosRes = new TProfile("ProfileCosRes","ProfileCosRes",nBinsCMore,binLimCMore);\r
   fProfileCosRes->Sumw2();\r
@@ -455,9 +618,13 @@ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()
   fListHist->Add(fSin2phiephiep);\r
   fListHist->Add(fCosRes);\r
   fListHist->Add(fCosResabc);\r
+  fListHist->Add(fSinRes);\r
+  fListHist->Add(fSinResabc);\r
   fListHist->Add(fDeltaPhiMaps);\r
   fListHist->Add(fCosPhiMaps);\r
   fListHist->Add(fProfileCosPhiMaps);\r
+\r
+  if(fHFEVZEROEventPlane && (fDebugLevel > 2)) fListHist->Add(fHFEVZEROEventPlane->GetOutputList());\r
   \r
 \r
   PostData(1, fListHist);\r
@@ -487,7 +654,9 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   Double_t valuensparseabis[5];\r
   Double_t valuensparsee[4];\r
   Double_t valuensparsef[2];\r
+  Double_t valuensparsefsin[2];\r
   Double_t valuensparsefbis[4];\r
+  Double_t valuensparsefbissin[4];\r
   Double_t valuensparseg[3];\r
   Double_t valuensparseh[3];\r
   Double_t valuensparsehprofile[3];\r
@@ -542,14 +711,7 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   if((95.0< cntr) && (cntr<100.0)) binctMore = 19.5;\r
 \r
   binctLess = cntr;\r
-  \r
-  //for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) {\r
-  //if((fBinCentralityLess[bincless]< cntr) && (cntr < fBinCentralityLess[bincless+1])) {\r
-  // binctLess = bincless+0.5;\r
-  //printf("fBinCentralityLess[bincless] %f and binctLess %f\n",fBinCentralityLess[bincless],binctLess);\r
-  //}\r
-  //}\r
-    \r
+   \r
   \r
   if(binct > 11.0) return;\r
  \r
@@ -558,7 +720,9 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   valuensparseabis[1] = binct;  \r
   valuensparsee[1] = binct;    \r
   valuensparsef[1] = binctMore;  \r
+  valuensparsefsin[1] = binct;  \r
   valuensparsefbis[3] = binctMore;  \r
+  valuensparsefbissin[3] = binct;  \r
   valuensparseg[1] = binct;\r
   valuensparseh[1] = binct; \r
   valuensparsehprofile[1] = binctLess; \r
@@ -615,14 +779,41 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   TVector2 *qsub2a = 0x0;\r
 \r
   // V0\r
\r
-  eventPlaneV0 = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0", esd,2));\r
-  if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
-  eventPlaneV0A = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0A", esd,2));\r
-  if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();\r
-  eventPlaneV0C = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0C", esd,2));\r
-  if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();\r
+\r
+  if(fHFEVZEROEventPlane){\r
+    \r
+    fHFEVZEROEventPlane->ProcessEvent(esd);\r
+    \r
+    if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0A()+100) < 0.0000001) eventPlaneV0A = -100.0;\r
+    else {\r
+      eventPlaneV0A = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0A());\r
+      if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();\r
+    }\r
+    \r
+    if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0C()+100) < 0.0000001) eventPlaneV0C = -100.0;\r
+    else {\r
+      eventPlaneV0C = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0C());\r
+      if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();\r
+    }\r
+\r
+    if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0()+100) < 0.0000001) eventPlaneV0 = -100.0;\r
+    else {\r
+      eventPlaneV0 = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0());\r
+      if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
+    }\r
+    \r
+  }\r
+  else {\r
+    \r
+    eventPlaneV0 = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0", esd,2));\r
+    if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
+    eventPlaneV0A = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0A", esd,2));\r
+    if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();\r
+    eventPlaneV0C = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0C", esd,2));\r
+    if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();\r
   \r
+  }\r
+\r
   // TPC\r
 \r
   standardQ = esdEPa->GetQVector(); \r
@@ -649,21 +840,32 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   Float_t eventPlanesub1a = -100.0;\r
   Float_t eventPlanesub2a = -100.0;\r
   Double_t diffsub1sub2a = -100.0;\r
+  Double_t diffsub1sub2asin = -100.0;\r
   Double_t diffsubasubb = -100.0;\r
   Double_t diffsubasubc = -100.0;\r
   Double_t diffsubbsubc = -100.0;\r
+  Double_t diffsubasubbsin = -100.0;\r
+  Double_t diffsubasubcsin = -100.0;\r
+  Double_t diffsubbsubcsin = -100.0;\r
 \r
   //if(fVZEROEventPlane) {\r
   diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C));\r
   diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC));\r
   diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC));\r
+\r
+  diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneV0C));\r
+  diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneTPC));\r
+  diffsubbsubcsin = TMath::Sin(2.*(eventPlaneV0C - eventPlaneTPC));\r
   //}\r
   //else {\r
   qsub1a = esdEPa->GetQsub1();\r
   qsub2a = esdEPa->GetQsub2();\r
   if(qsub1a) eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.;\r
   if(qsub2a) eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.;\r
-  if(qsub1a && qsub2a) diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
+  if(qsub1a && qsub2a) {\r
+    diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
+    diffsub1sub2asin = TMath::Sin(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
+  }\r
   //}\r
   \r
 \r
@@ -736,7 +938,9 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
     if(!fVZEROEventPlane) {\r
       valuensparsef[0] = diffsub1sub2a;\r
       fCosRes->Fill(&valuensparsef[0]);\r
-      if(fDebugLevel > 0) {\r
+      valuensparsefsin[0] = diffsub1sub2asin;\r
+      fSinRes->Fill(&valuensparsefsin[0]);\r
+      if(fDebugLevel > 1) {\r
        fProfileCosRes->Fill(valuensparsef[1],valuensparsef[0]);\r
       }\r
     }\r
@@ -745,7 +949,11 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
       valuensparsefbis[1] = diffsubbsubc;\r
       valuensparsefbis[2] = diffsubasubc;\r
       fCosResabc->Fill(&valuensparsefbis[0]);\r
-      if(fDebugLevel > 0) {\r
+      valuensparsefbissin[0] = diffsubasubbsin;\r
+      valuensparsefbissin[1] = diffsubbsubcsin;\r
+      valuensparsefbissin[2] = diffsubasubcsin;\r
+      fSinResabc->Fill(&valuensparsefbissin[0]);\r
+      if(fDebugLevel > 1) {\r
        fProfileCosResab->Fill(valuensparsefbis[3],valuensparsefbis[0]);\r
        fProfileCosResac->Fill(valuensparsefbis[3],valuensparsefbis[1]);\r
        fProfileCosResbc->Fill(valuensparsefbis[3],valuensparsefbis[2]);\r
@@ -775,23 +983,26 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
     \r
     if(!survived) continue;\r
 \r
-    // Apply PID for Data\r
-    if(!fMCPID) {\r
-      AliHFEpidObject hfetrack;\r
-      hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);\r
-      hfetrack.SetRecTrack(track);\r
-      hfetrack.SetCentrality((Int_t)binct);\r
-      //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality());\r
-      hfetrack.SetPbPb();\r
-      if(!fPID->IsSelected(&hfetrack,0x0,"",fPIDqa)) {\r
-       continue;\r
+    // Apply PID\r
+    if(!fNoPID) {\r
+      // Apply PID for Data\r
+      if(!fMCPID) {\r
+       AliHFEpidObject hfetrack;\r
+       hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);\r
+       hfetrack.SetRecTrack(track);\r
+       hfetrack.SetCentrality((Int_t)binct);\r
+       //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality());\r
+       hfetrack.SetPbPb();\r
+       if(!fPID->IsSelected(&hfetrack,0x0,"",fPIDqa)) {\r
+         continue;\r
+       }\r
+      }\r
+      else {\r
+       if(!mcEvent) continue;\r
+       if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;\r
+       //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode()));\r
+       if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;\r
       }\r
-    }\r
-    else {\r
-      if(!mcEvent) continue;\r
-      if(!(mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;\r
-      //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode()));\r
-      if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue;\r
     }\r
 \r
 \r
@@ -894,7 +1105,7 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
     if(fillEventPlane) fEventPlaneaftersubtraction->Fill(&valuensparseabis[0]);\r
     \r
 \r
-    if(fDebugLevel > 1) \r
+    if(fDebugLevel > 3) \r
       {\r
        //\r
        valuensparsee[0] = TMath::Cos(2*phitrack);\r
index ebf9895..7850d00 100644 (file)
@@ -36,25 +36,31 @@ class TProfile2D;
 class THnSparse;\r
 class AliHFEpidQAmanager;\r
 class AliFlowEvent;\r
+class AliHFEVZEROEventPlane;\r
 \r
 class AliAnalysisTaskHFEFlow: public AliAnalysisTaskSE {\r
   \r
 public:\r
   AliAnalysisTaskHFEFlow();\r
   AliAnalysisTaskHFEFlow(const char *name);\r
-  //TODO: if we get the ESDpid object from somewhere else, we should not delete it anymore!!!\r
-  virtual ~AliAnalysisTaskHFEFlow(){;}\r
+  AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref);\r
+  AliAnalysisTaskHFEFlow& operator=(const AliAnalysisTaskHFEFlow &ref);\r
+  virtual void Copy(TObject &o) const;\r
+  virtual ~AliAnalysisTaskHFEFlow();\r
   \r
   virtual void  UserExec(Option_t */*option*/);\r
   virtual void  UserCreateOutputObjects();\r
 \r
   AliHFEpid *GetPID() const { return fPID; }\r
+  AliHFEpidQAmanager *GetPIDQAManager() const { return fPIDqa; }\r
+\r
   void SetHFECuts(AliHFEcuts * const cuts) { fHFECuts = cuts; };\r
   void SetSubEtaGapTPC(Bool_t  subEtaGapTPC) { fSubEtaGapTPC = subEtaGapTPC; };\r
   void SetEtaGap(Double_t  etaGap) { fEtaGap = etaGap; };\r
   void SetVZEROEventPlane(Bool_t vzeroEventPlane) { fVZEROEventPlane = vzeroEventPlane; };\r
   void SetVZEROEventPlaneA(Bool_t vzeroEventPlaneA) { fVZEROEventPlaneA = vzeroEventPlaneA; };\r
   void SetVZEROEventPlaneC(Bool_t vzeroEventPlaneC) { fVZEROEventPlaneC = vzeroEventPlaneC; };\r
+  void SetHFEVZEROEventPlane(AliHFEVZEROEventPlane *hfeVZEROEventPlane) { fHFEVZEROEventPlane = hfeVZEROEventPlane; };\r
 \r
   void SetNbBinsCentralityQCumulant(Int_t nbBinsCentralityQCumulant) { fNbBinsCentralityQCumulant = nbBinsCentralityQCumulant; };\r
   void SetBinCentralityLess(Int_t k, Float_t value)  { fBinCentralityLess[k] = value; };\r
@@ -69,6 +75,7 @@ public:
   void SetPrecisionPhi(Double_t precisionPhi) { fPrecisionPhi = precisionPhi;};\r
   void SetUseMCReactionPlane(Bool_t useMCReactionPlane) { fUseMCReactionPlane = useMCReactionPlane;};\r
   void SetMCPID(Bool_t mcPID) { fMCPID = mcPID;};\r
+  void SetNoPID(Bool_t noPID) { fNoPID = noPID;};\r
 \r
   void SetDebugLevel(Int_t debugLevel) { fDebugLevel = debugLevel;};\r
 \r
@@ -79,21 +86,21 @@ public:
   Double_t GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const;\r
   \r
 private:\r
-  TList       *fListHist;  //TH list\r
+  TList       *fListHist;       //! TH list\r
 \r
-  Bool_t    fVZEROEventPlane; // Use Event Planes from VZERO\r
+  Bool_t    fVZEROEventPlane;  // Use Event Planes from VZERO\r
   Bool_t    fVZEROEventPlaneA; // Use Event Planes from VZERO A\r
   Bool_t    fVZEROEventPlaneC; // Use Event Planes from VZERO C\r
 \r
-  Bool_t    fSubEtaGapTPC;  // bool to fill with eta gap\r
-  Double_t  fEtaGap;       // Value of the eta gap\r
+  Bool_t    fSubEtaGapTPC;    // bool to fill with eta gap\r
+  Double_t  fEtaGap;          // Value of the eta gap\r
 \r
   Int_t     fNbBinsCentralityQCumulant;  // Number of Bins Q Cumulant\r
-  Double_t  fBinCentralityLess[10]; // Centrality Bin lower value\r
-  Int_t     fNbBinsPtQCumulant; // Nbbinspt QCumulant method\r
-  Double_t  fMinPtQCumulant;   // Min pt QCumulant method\r
-  Double_t  fMaxPtQCumulant;   // Max pt QCumulant method\r
-  Bool_t    fAfterBurnerOn;    // Add flow to all tracks\r
+  Double_t  fBinCentralityLess[10];      // Centrality Bin lower value\r
+  Int_t     fNbBinsPtQCumulant;          // Nbbinspt QCumulant method\r
+  Double_t  fMinPtQCumulant;             // Min pt QCumulant method\r
+  Double_t  fMaxPtQCumulant;             // Max pt QCumulant method\r
+  Bool_t    fAfterBurnerOn;              // Add flow to all tracks\r
   Int_t     fNonFlowNumberOfTrackClones; // number of times to clone the particles (nonflow) \r
   Double_t  fV1;        // Add Flow. Must be in range [0,0.5].\r
   Double_t  fV2;        // Add Flow. Must be in range [0,0.5].\r
@@ -105,6 +112,7 @@ private:
   Bool_t    fUseMCReactionPlane; // use MC reaction plane\r
 \r
   Bool_t    fMCPID; // MC PID for electrons\r
+  Bool_t    fNoPID; // No PID for checks\r
 \r
   Int_t     fDebugLevel; // Debug Level  \r
 \r
@@ -113,53 +121,54 @@ private:
   AliFlowTrackCuts* fcutsPOI; // Particle Of Interest cut\r
   \r
   // Cuts for HFE\r
-  AliHFEcuts *fHFECuts;       // HFE cuts\r
-  AliHFEpid  *fPID;           // PID cuts \r
-  AliHFEpidQAmanager *fPIDqa;   // QA Manager\r
-  AliFlowEvent *fflowEvent;     // Flow event   \r
-  \r
+  AliHFEcuts *fHFECuts;           // HFE cuts\r
+  AliHFEpid  *fPID;               // PID cuts \r
+  AliHFEpidQAmanager *fPIDqa;     // QA Manager\r
+  AliFlowEvent *fflowEvent;       //! Flow event   \r
 \r
+  // VZERO Event plane after calibration 2010\r
+  AliHFEVZEROEventPlane *fHFEVZEROEventPlane; // VZERO event plane calibrated\r
+  \r
   // Histos\r
-  TH2D *fHistEV;               // Number of events\r
+  TH2D *fHistEV;               //! Number of events\r
   \r
   // A Event plane as function of phiepa, phiepb, phiepc, phiepd centrality \r
   // a V0A, b V0C, c TPC, d V0\r
-  THnSparseF *fEventPlane;     // Event plane\r
+  THnSparseF *fEventPlane;     //! Event plane\r
   \r
   // B Event Plane after subtraction as function of phiep, centrality \r
-  THnSparseF *fEventPlaneaftersubtraction; // Event plane\r
+  THnSparseF *fEventPlaneaftersubtraction; //! Event plane\r
 \r
   // Monitoring Event plane: cos2phi, sin2phi, centrality\r
-  THnSparseF *fCosSin2phiep;        // Cos(2phi), Sin(2phi)\r
+  THnSparseF *fCosSin2phiep;        //! Cos(2phi), Sin(2phi)\r
   \r
   // E Monitoring Event plane after subtraction of the track: cos, centrality, pt, eta\r
-  THnSparseF *fCos2phie;  // Monitoring\r
-  THnSparseF *fSin2phie;  // Monitoring\r
-  THnSparseF *fCos2phiep;  // Monitoring\r
-  THnSparseF *fSin2phiep;  // Monitoring\r
-  THnSparseF *fSin2phiephiep;  // Monitoring\r
+  THnSparseF *fCos2phie;  //! Monitoring\r
+  THnSparseF *fSin2phie;  //! Monitoring\r
+  THnSparseF *fCos2phiep;  //! Monitoring\r
+  THnSparseF *fSin2phiep;  //! Monitoring\r
+  THnSparseF *fSin2phiephiep;  //! Monitoring\r
 \r
   // Fbis Resolution as function of cosres, cosres, cosres, centrality for three subevents (V0)\r
   // a V0A, b V0C, c TPC\r
-  THnSparseF *fCosResabc; // Res\r
-  TProfile   *fProfileCosResab; // Profile Res_a_b\r
-  TProfile   *fProfileCosResac; // Profile Res_a_c\r
-  TProfile   *fProfileCosResbc; // Profile Res_b_c\r
+  THnSparseF *fCosResabc; //! Res\r
+  THnSparseF *fSinResabc; //! Res\r
+  TProfile   *fProfileCosResab; //! Profile Res_a_b\r
+  TProfile   *fProfileCosResac; //! Profile Res_a_c\r
+  TProfile   *fProfileCosResbc; //! Profile Res_b_c\r
   \r
   // F Resolution as function of cosres, centrality for two subevents (TPC)\r
-  THnSparseF *fCosRes; // Res\r
-  TProfile   *fProfileCosRes; // Profile Res\r
+  THnSparseF *fCosRes; //! Res\r
+  THnSparseF *fSinRes; //! Res\r
+  TProfile   *fProfileCosRes; //! Profile Res\r
   \r
   // G Maps delta phi as function of deltaphi, centrality, pt\r
-  THnSparseF *fDeltaPhiMaps; // Delta phi\r
+  THnSparseF *fDeltaPhiMaps; //! Delta phi\r
   \r
   // H Maps cos phi : cos, centrality, pt\r
-  THnSparseF *fCosPhiMaps;         //Cos\r
-  TProfile2D *fProfileCosPhiMaps;  // Profile Cos\r
-  \r
+  THnSparseF *fCosPhiMaps;         //! Cos\r
+  TProfile2D *fProfileCosPhiMaps;  //! Profile Cos\r
   \r
-  AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow&); // not implemented\r
-  AliAnalysisTaskHFEFlow& operator=(const AliAnalysisTaskHFEFlow&); // not implemented\r
   \r
   ClassDef(AliAnalysisTaskHFEFlow, 1); // analysisclass\r
 };\r
diff --git a/PWGHF/hfe/AliHFEVZEROEventPlane.cxx b/PWGHF/hfe/AliHFEVZEROEventPlane.cxx
new file mode 100644 (file)
index 0000000..3580da9
--- /dev/null
@@ -0,0 +1,407 @@
+#include "AliHFEVZEROEventPlane.h"
+
+// AliRoot includes
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliESDtrack.h"
+#include "AliCentrality.h"
+#include "AliESDVZERO.h"
+#include "TFile.h"
+#include "AliOADBContainer.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TList.h"
+#include "TChain.h"
+#include "THnSparse.h"
+#include "TString.h"
+#include "TVector2.h"
+#include "AliEventplane.h"
+
+
+// STL includes
+#include <iostream>
+using namespace std;
+
+
+//_____________________________________________________________________________
+AliHFEVZEROEventPlane::AliHFEVZEROEventPlane():
+  TNamed(),
+  fEventPlaneV0A(-100.0),
+  fEventPlaneV0C(-100.0),
+  fEventPlaneV0(-100.0),
+  fESD(0),
+  fRun(-1),
+  fMultV0(0x0),
+  fV0Cpol(100),
+  fV0Apol(100),
+  fnamefile(TString("")),
+  fOutputList(0x0),
+  fMultV0Before(0x0),
+  fMultV0After(0x0)
+{
+  //
+  // Default constructor (should not be used)
+  //
+
+  for(Int_t k = 0; k < fgknCentrBin; k++) {
+    for(Int_t iside = 0; iside < 2; iside++) {
+      for(Int_t icoord = 0; icoord < 2; icoord++) {
+       fQBefore[k][iside][icoord] = 0x0;
+       fQAfter[k][iside][icoord] = 0x0;
+      }
+    }
+  }
+  
+}
+
+//______________________________________________________________________________
+AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const char *name, const Char_t *title):
+  TNamed(name,title),
+  fEventPlaneV0A(-100.0),
+  fEventPlaneV0C(-100.0),
+  fEventPlaneV0(-100.0),
+  fESD(0),
+  fRun(-1),
+  fMultV0(0x0),
+  fV0Cpol(100),
+  fV0Apol(100),
+  fnamefile(TString("")),
+  fOutputList(0x0),
+  fMultV0Before(0x0),
+  fMultV0After(0x0)
+{
+  //
+  // Constructor
+  //
+
+  char namelist[100];
+  sprintf(namelist,"QA_hist_%s",name);
+  fOutputList = new TList();
+  fOutputList->SetName((const char*)namelist);
+  fOutputList->SetOwner(kTRUE);
+
+  // Multiplicity
+  fMultV0Before = new TProfile("MultV0Before","",64,0,64);
+  fMultV0Before->Sumw2();
+  fMultV0After = new TProfile("MultV0After","",64,0,64);
+  fMultV0After->Sumw2();
+  fOutputList->Add(fMultV0Before);
+  fOutputList->Add(fMultV0After);
+
+  // Recentering
+  char namecontbefore[100];
+  char namecontafter[100];    
+  for(Int_t k = 0; k < fgknCentrBin; k++) {
+    for(Int_t iside = 0; iside < 2; iside++) {
+      for(Int_t icoord = 0; icoord < 2; icoord++) {
+       
+       if(iside==0 && icoord==0)
+         sprintf(namecontbefore,"hQxc2_%i",k);
+       else if(iside==1 && icoord==0)
+         sprintf(namecontbefore,"hQxa2_%i",k);
+       else if(iside==0 && icoord==1)
+         sprintf(namecontbefore,"hQyc2_%i",k);
+       else if(iside==1 && icoord==1)
+         sprintf(namecontbefore,"hQya2_%i",k);
+       //
+       if(iside==0 && icoord==0)
+         sprintf(namecontafter,"hQxc2_%i_after",k);
+       else if(iside==1 && icoord==0)
+         sprintf(namecontafter,"hQxa2_%i_after",k);
+       else if(iside==0 && icoord==1)
+         sprintf(namecontafter,"hQyc2_%i_after",k);
+       else if(iside==1 && icoord==1)
+         sprintf(namecontafter,"hQya2_%i_after",k);
+       
+       //
+       fQBefore[k][iside][icoord] = new TH1F(((const char*)namecontbefore),"",800,-400.0,400.0);
+       fQBefore[k][iside][icoord]->Sumw2();
+       fOutputList->Add(fQBefore[k][iside][icoord]);
+       fQAfter[k][iside][icoord] = new TH1F(((const char*)namecontafter),"",800,-400.0,400.0);
+       fQAfter[k][iside][icoord]->Sumw2();
+       fOutputList->Add(fQAfter[k][iside][icoord]);
+       //
+      }
+    }
+  }
+  
+}
+//____________________________________________________________
+AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const AliHFEVZEROEventPlane &ref):
+  TNamed(ref),
+  fEventPlaneV0A(ref.fEventPlaneV0A),
+  fEventPlaneV0C(ref.fEventPlaneV0C),
+  fEventPlaneV0(ref.fEventPlaneV0),
+  fESD(0),
+  fRun(-1),
+  fMultV0(0x0),
+  fV0Cpol(100),
+  fV0Apol(100),
+  fnamefile(ref.fnamefile),
+  fOutputList(0x0),
+  fMultV0Before(0x0),
+  fMultV0After(0x0)
+{
+  //
+  // Copy Constructor
+  //
+  ref.Copy(*this);
+}
+
+//____________________________________________________________
+AliHFEVZEROEventPlane &AliHFEVZEROEventPlane::operator=(const AliHFEVZEROEventPlane &ref){
+  //
+  // Assignment operator
+  //
+  if(this == &ref) 
+    ref.Copy(*this);
+  return *this;
+}
+
+//____________________________________________________________
+void AliHFEVZEROEventPlane::Copy(TObject &o) const {
+  // 
+  // Copy into object o
+  //
+  AliHFEVZEROEventPlane &target = dynamic_cast<AliHFEVZEROEventPlane &>(o);
+
+  target.fEventPlaneV0A = fEventPlaneV0A;
+  target.fEventPlaneV0C = fEventPlaneV0C;
+  target.fEventPlaneV0 = fEventPlaneV0;
+  target.fnamefile = fnamefile;
+  
+}
+
+//__________________________________________________________________
+AliHFEVZEROEventPlane::~AliHFEVZEROEventPlane(){
+  //
+  // Destruktor
+  //
+  if(fOutputList){
+    fOutputList->Delete();
+    delete fOutputList;
+  }
+  fOutputList = 0x0;
+  if(fMultV0Before) delete fMultV0Before;
+  if(fMultV0After) delete fMultV0After;
+
+  for(Int_t k = 0; k < fgknCentrBin; k++) {
+    for(Int_t iside = 0; iside < 2; iside++) {
+      for(Int_t icoord = 0; icoord < 2; icoord++) {
+       if(fQBefore[k][iside][icoord]) delete fQBefore[k][iside][icoord];
+       if(fQAfter[k][iside][icoord]) delete fQAfter[k][iside][icoord];
+      }
+    }
+  }
+
+}
+//______________________________________________________________________________
+void AliHFEVZEROEventPlane::ProcessEvent(AliESDEvent *event) 
+{
+  //
+  // Process the event
+  // 
+  
+  // Reset
+  fEventPlaneV0A = -100.0;
+  fEventPlaneV0C = -100.0;
+  fEventPlaneV0  = -100.0;
+  
+  //
+  Int_t run = event->GetRunNumber();
+  Bool_t ok = kTRUE;
+  if(run != fRun){
+    if(fMultV0) delete fMultV0;
+    fMultV0 = 0x0;
+    // Load the calibrations run dependent
+    if(!OpenInfoCalbration(run)) ok = kFALSE;
+    fRun=run;
+  }
+
+  if(ok) Analyze(event);
+
+}
+
+//________________________________________________________________________
+void AliHFEVZEROEventPlane::Analyze(AliESDEvent* esdEvent)
+{  
+  //
+  // Do VZERO calibration + centering
+  //
+
+  if(!fMultV0) {
+    //printf("Did not find calibration VZERO\n");
+    return;
+  }
+  //printf("HERE!!!\n");
+
+  //Centrality
+  Float_t v0Centr  = -10.;
+  AliCentrality* centrality = esdEvent->GetCentrality();
+  if (centrality){
+    v0Centr  = centrality->GetCentralityPercentile("V0M");
+  }
+  
+  // Analyse only for 0-80% PbPb collisions
+  Int_t iC = -1;    
+  if (v0Centr >0 && v0Centr < 80){
+    if(v0Centr < 5) iC = 0;
+    else if(v0Centr < 10) iC = 1;
+    else if(v0Centr < 20) iC = 2;
+    else if(v0Centr < 30) iC = 3;
+    else if(v0Centr < 40) iC = 4;
+    else if(v0Centr < 50) iC = 5;
+    else if(v0Centr < 60) iC = 6;
+    else if(v0Centr < 70) iC = 7;
+    else iC = 8;
+    
+    //general  
+    Double_t qxa2 = 0, qya2 = 0;
+    Double_t qxc2 = 0, qyc2 = 0;
+    
+    //V0 info    
+    AliESDVZERO* esdV0 = esdEvent->GetVZEROData();
+    
+    for (Int_t iv0 = 0; iv0 < 64; iv0++) {
+      Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
+      Float_t multv0 = esdV0->GetMultiplicity(iv0);
+      if (iv0 < 32){
+       //printf("Value %f\n",fMultV0->GetBinContent(iv0+1));
+       if(fMultV0->GetBinContent(iv0+1)>0.0) {
+         qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+         qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
+         fMultV0After->Fill(iv0,multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1));
+       }
+      } else {
+       if(fMultV0->GetBinContent(iv0+1)>0.0) {
+         qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+         qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
+         fMultV0After->Fill(iv0,multv0*fV0Apol/fMultV0->GetBinContent(iv0+1));
+       }
+      }
+      fMultV0Before->Fill(iv0,multv0);
+    }
+
+    // Fill histos
+    fQBefore[iC][1][0] -> Fill(qxa2);
+    fQBefore[iC][1][1] -> Fill(qya2);
+    fQBefore[iC][0][0] -> Fill(qxc2);
+    fQBefore[iC][0][1] -> Fill(qyc2);
+    
+    //grab for each centrality the proper histo with the Qx and Qy to do the recentering
+    Double_t qxamean2 = fMeanQ[iC][1][0];
+    Double_t qxarms2  = fWidthQ[iC][1][0];
+    Double_t qyamean2 = fMeanQ[iC][1][1];
+    Double_t qyarms2  = fWidthQ[iC][1][1];
+    
+    Double_t qxcmean2 = fMeanQ[iC][0][0];
+    Double_t qxcrms2  = fWidthQ[iC][0][0];
+    Double_t qycmean2 = fMeanQ[iC][0][1];
+    Double_t qycrms2  = fWidthQ[iC][0][1];     
+
+    if((TMath::Abs(qxarms2) < 0.00000001) || (TMath::Abs(qyarms2) < 0.00000001) || (TMath::Abs(qxcrms2) < 0.00000001) || (TMath::Abs(qycrms2) < 0.00000001)) return;    
+
+    Double_t qxaCor2 = (qxa2 - qxamean2)/qxarms2;
+    Double_t qyaCor2 = (qya2 - qyamean2)/qyarms2;
+    Double_t qxcCor2 = (qxc2 - qxcmean2)/qxcrms2;
+    Double_t qycCor2 = (qyc2 - qycmean2)/qycrms2;
+
+    Double_t qxCor2 = qxaCor2 + qxcCor2;
+    Double_t qyCor2 = qyaCor2 + qycCor2;
+
+    // Fill histos
+    fQAfter[iC][1][0] -> Fill(qxaCor2);
+    fQAfter[iC][1][1] -> Fill(qyaCor2);
+    fQAfter[iC][0][0] -> Fill(qxcCor2);
+    fQAfter[iC][0][1] -> Fill(qycCor2);
+    
+    Double_t evPlAngV0ACor2 = TVector2::Phi_0_2pi(TMath::ATan2(qyaCor2, qxaCor2))/2.;
+    Double_t evPlAngV0CCor2 = TVector2::Phi_0_2pi(TMath::ATan2(qycCor2, qxcCor2))/2.;
+    Double_t evPlAngV0Cor2  = TVector2::Phi_0_2pi(TMath::ATan2(qyCor2, qxCor2))/2.;
+
+    fEventPlaneV0A = evPlAngV0ACor2;
+    fEventPlaneV0C = evPlAngV0CCor2;
+    fEventPlaneV0  = evPlAngV0Cor2;
+    
+    //printf("Eventplane V0A %f, V0C %f\n",fEventPlaneV0A,fEventPlaneV0C);
+
+  }
+}
+//_____________________________________________________________________________
+Bool_t AliHFEVZEROEventPlane::OpenInfoCalbration(Int_t run)
+{
+  //
+  // Take the calibration coefficients
+  //  
+
+  //printf("Name of the file %s\n",(const char*)fnamefile);
+  TFile *foadb = TFile::Open(fnamefile.Data());
+  if(!foadb){
+    printf("OADB file %s cannot be opened\n",fnamefile.Data());
+    return kFALSE;
+  }
+
+  //printf("test\n");
+
+  AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
+  if(!cont){
+    printf("OADB object hMultV0BefCorr is not available in the file\n");
+    return kFALSE;     
+  }
+
+  if(!(cont->GetObject(run))){
+    printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
+    return kFALSE;     
+  }
+  //TProfile *multV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
+  TProfile *multV0 = ((TProfile *) cont->GetObject(run));
+  if(!multV0) return kFALSE;
+  fMultV0 = (TProfile *) multV0->Clone();
+  fMultV0->SetDirectory(0);
+  
+  TF1 *fpol0 = new TF1("fpol0","pol0"); 
+  fMultV0->Fit(fpol0,"Q0","",0,31);
+  fV0Cpol = fpol0->GetParameter(0);
+  fMultV0->Fit(fpol0,"Q0","",32,64);
+  fV0Apol = fpol0->GetParameter(0);
+  
+  for(Int_t iside=0;iside<2;iside++){
+    for(Int_t icoord=0;icoord<2;icoord++){
+      for(Int_t i=0;i  < fgknCentrBin;i++){
+       char namecont[100];
+       if(iside==0 && icoord==0)
+         sprintf(namecont,"hQxc2_%i",i);
+       else if(iside==1 && icoord==0)
+         sprintf(namecont,"hQxa2_%i",i);
+       else if(iside==0 && icoord==1)
+         sprintf(namecont,"hQyc2_%i",i);
+       else if(iside==1 && icoord==1)
+         sprintf(namecont,"hQya2_%i",i);
+       
+       cont = (AliOADBContainer*) foadb->Get(namecont);
+       if(!cont){
+         printf("OADB object %s is not available in the file\n",namecont);
+         delete fpol0;
+         return kFALSE;        
+       }
+       
+       if(!(cont->GetObject(run))){
+         printf("OADB object %s is not available for run %i\n",namecont,run);
+         delete fpol0;
+         return kFALSE;        
+       }
+       fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
+       fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
+      }
+    }
+  }
+  
+  delete fpol0;
+  foadb->Close();
+  
+  return kTRUE;
+
+}
diff --git a/PWGHF/hfe/AliHFEVZEROEventPlane.h b/PWGHF/hfe/AliHFEVZEROEventPlane.h
new file mode 100644 (file)
index 0000000..e0f1ac3
--- /dev/null
@@ -0,0 +1,59 @@
+//
+// VZERO event plane task for 2010
+//
+#ifndef AliHFEVZEROEVENTPLANE_H
+#define AliHFEVZEROEVENTPLANE_H
+
+#include "TProfile.h"
+#include <AliESDEvent.h>
+#include "TString.h"
+#include "TH1F.h"
+#include "TList.h"
+
+
+class AliHFEVZEROEventPlane : public TNamed {
+ public:
+  AliHFEVZEROEventPlane();
+  AliHFEVZEROEventPlane(const char *name, const Char_t *title);
+  AliHFEVZEROEventPlane(const AliHFEVZEROEventPlane &ref);
+  AliHFEVZEROEventPlane& operator=(const AliHFEVZEROEventPlane &ref);
+  virtual void Copy(TObject &o) const;
+  ~AliHFEVZEROEventPlane();
+
+  void ProcessEvent(AliESDEvent *event);
+  
+  void  SetNameFile(TString namefile) {fnamefile = namefile;};
+  Bool_t OpenInfoCalbration(Int_t run);
+
+  Double_t GetEventPlaneV0A() const {return fEventPlaneV0A;};
+  Double_t GetEventPlaneV0C() const {return fEventPlaneV0C;};
+  Double_t GetEventPlaneV0()  const {return fEventPlaneV0;};
+
+  TList *GetOutputList() const {return fOutputList;};
+
+ private:
+  virtual void Analyze(AliESDEvent* esdEvent); 
+
+  Double_t fEventPlaneV0A;          // Corrected event plane V0A
+  Double_t fEventPlaneV0C;          // Corrected event plane V0C
+  Double_t fEventPlaneV0;           // Corrected event plane V0
+
+  AliESDEvent* fESD;                //! ESD object
+  Int_t        fRun;                // Run number
+  TProfile *fMultV0;                //! fMultiplicityV0
+  Float_t fV0Cpol,fV0Apol;          // fV0Cpol, fV0Apol
+  static const Int_t fgknCentrBin = 9; // Centrality bins
+  Float_t fMeanQ[fgknCentrBin][2][2];  // mean for centering
+  Float_t fWidthQ[fgknCentrBin][2][2]; // rms for centering
+  TString fnamefile;                // name of the file with the coefficient
+  TList       *fOutputList;         //! Output list
+  TProfile    *fMultV0Before;       //! fMultiplicityV0 Before
+  TProfile    *fMultV0After;        //! fMultiplicityV0 After
+  TH1F *fQBefore[fgknCentrBin][2][2]; //! Q centering before
+  TH1F *fQAfter[fgknCentrBin][2][2];  //! Q centering after
+
+  ClassDef(AliHFEVZEROEventPlane, 1);    //Analysis task for high pt analysis 
+};
+
+#endif
index 1f2a481..75f93b7 100644 (file)
@@ -150,6 +150,23 @@ Bool_t AliHFEcollection::CreateTH1Farray(const char* name, const char* title, In
 }
 
 //___________________________________________________________________
+Bool_t AliHFEcollection::CreateTH2Farray(const char* name, const char* title, Int_t nBin, const Double_t* xbins, Int_t nBinY, Float_t nMinY, Float_t nMaxY){
+
+  //
+  // Creates a TH1F histogram for the collection 2nd version 
+  //
+
+  if(!fList){
+    AliError("No TList pointer ! ");
+    return kFALSE;
+  }
+  else{
+    fList->Add(new TH2F(name, title, nBin, xbins, nBinY, nMinY, nMaxY));
+    return CheckObject(name);
+  }
+}
+
+//___________________________________________________________________
 Bool_t AliHFEcollection::CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t logAxis){
 
   //
@@ -303,6 +320,21 @@ Bool_t AliHFEcollection::CreateTHnSparse(const char* name, const char* title, In
 
 }
 //___________________________________________________________________
+Bool_t AliHFEcollection::CreateTHnSparseNoLimits(const char* name, const char* title, Int_t dim, const Int_t* nbins){
+
+  //
+  // create 'dim' dimensional THnSparse without limits
+  //
+
+  if(!fList){
+    AliError("No TList pointer ! ");
+    return kFALSE;
+  }
+  fList->Add(new THnSparseF(name, title, dim, nbins));
+  return CheckObject(name);
+
+}
+//___________________________________________________________________
 TObject* AliHFEcollection::Get(const char* name){ 
 
   //
index 5de56ce..1464f0f 100644 (file)
@@ -55,6 +55,7 @@ class AliHFEcollection : public TNamed{
   Bool_t CreateTH1Farray(const char* name, const char* title, Int_t nBin, const Double_t* xbins);
 
   Bool_t CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t logAxis = -1);
+  Bool_t CreateTH2Farray(const char* name, const char* title, Int_t nBin, const Double_t* xbins, Int_t nBinY, Float_t nMinY, Float_t nMaxY);
   Bool_t CreateTH3F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t nBinZ, Float_t minZ, Float_t maxZ, Int_t logAxis = -1);
 
   Bool_t CreateTH1Fvector1(Int_t X, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax, Int_t logAxis = -1);
@@ -62,6 +63,7 @@ class AliHFEcollection : public TNamed{
   Bool_t CreateTH2Fvector1(Int_t X, const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t logAxis = -1);
   Bool_t CreateProfile(const char* name, const char* title, Int_t nbins, Double_t xmin, Double_t xmax);
   Bool_t CreateTHnSparse(const char* name, const char* title, Int_t dim, const Int_t* nbins, const Double_t* xmin, const Double_t* xmax);
+  Bool_t CreateTHnSparseNoLimits(const char* name, const char* title, Int_t dim, const Int_t* nbins);
 
   Bool_t BinLogAxis(const char* name, Int_t dim);
   Bool_t Sumw2(const char*name);
index 9edd39d..234cf90 100644 (file)
@@ -18,7 +18,6 @@
 // Authors:
 //   Markus Fasel <M.Fasel@gsi.de>
 //
-//
 #include <TBits.h>
 #include <TString.h>
 
@@ -37,6 +36,7 @@
 #include "AliMCParticle.h"
 #include "AliPIDResponse.h"
 #include "AliVEvent.h"
+#include "AliHFEpidTRD.h"
 #include "TTreeStream.h"
 
 #include "AliHFEdebugTreeTask.h"
@@ -47,6 +47,7 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask():
   AliAnalysisTaskSE(),
   fTrackCuts(NULL),
   fSignalCuts(NULL),
+  fTRDpid(NULL),
   fNclustersTPC(70),
   fNclustersTPCPID(0),
   fNclustersITS(2),
@@ -60,6 +61,7 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name):
   AliAnalysisTaskSE(name),
   fTrackCuts(NULL),
   fSignalCuts(NULL),
+  fTRDpid(NULL),
   fNclustersTPC(70),
   fNclustersTPCPID(0),
   fNclustersITS(2),
@@ -70,7 +72,8 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name):
 }
 
 AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){
-  if(fDebugTree) delete fDebugTree;
+    if(fDebugTree) delete fDebugTree;
+    if(fTRDpid) delete fTRDpid;
 }
 
 void AliHFEdebugTreeTask::UserCreateOutputObjects(){
@@ -93,6 +96,9 @@ void AliHFEdebugTreeTask::UserCreateOutputObjects(){
   fTrackCuts->SetUseMixedVertex(kTRUE);
   fTrackCuts->SetVertexRange(10.);
   fTrackCuts->Initialize();
+
+  fTRDpid = new AliHFEpidTRD("QAtrdPID");
+
 }
 
 void AliHFEdebugTreeTask::UserExec(Option_t *){
@@ -281,13 +287,16 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
     // TRD clusters and tracklets
     UChar_t nclustersTRD = track->GetTRDncls();
-    UChar_t ntracklets = track->GetTRDntrackletsPID();
+    UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
     // ITS and TRD acceptance maps
     UChar_t hasClusterITS[6], hasTrackletTRD[6];
     UChar_t itsPixel = track->GetITSClusterMap();
     for(Int_t icl = 0; icl < 6; icl++) hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0;
+    Double_t trddEdxSum[6];
+    for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;}
     for(Int_t itl = 0; itl < 6; itl++){
-      Int_t nSliceNonZero = 0;
+       Int_t nSliceNonZero = 0;
+        trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge
       for(Int_t islice = 0; islice < 8; islice++){
         if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++;
       }
@@ -298,6 +307,9 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
     track->GetTRDpid(pidprobs);
     Double_t likeEleTRD = pidprobs[0];
     Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]);
+    Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6);
+    Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6);
+
     // DCA
     Float_t b[2] = {0.,0.};
     Float_t bCov[3] = {0.,0.,0.};
@@ -326,7 +338,7 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
                   << "pt="                  << transversemomentum
                   << "eta="                 << eta
                   << "phi="                 << phi
-                  << "ntracklets="          << ntracklets
+                  << "ntracklets="          << ntrackletsTRDPID
                   << "nclustersTPC="        << nclustersTPC
                   << "nclustersTPCPID="     << nclustersTPCPID
                   << "nclustersTPCshared="  << nclustersTPCshared
@@ -348,11 +360,19 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
                   << "trd3="                << hasTrackletTRD[3]
                   << "trd4="                << hasTrackletTRD[4]
                   << "trd5="                << hasTrackletTRD[5]
+                 << "TRDdEdxl0="           << trddEdxSum[0]
+                 << "TRDdEdxl1="           << trddEdxSum[1]
+                 << "TRDdEdxl2="           << trddEdxSum[2]
+                 << "TRDdEdxl3="           << trddEdxSum[3]
+                 << "TRDdEdxl4="           << trddEdxSum[4]
+                 << "TRDdEdxl5="           << trddEdxSum[5]
                   << "TOFsigmaEl="          << nSigmaTOF
                   << "TPCsigmaEl="          << nSigmaTPC
                   << "TPCdEdx="             << tPCdEdx
                   << "TRDlikeEl="           << likeEleTRD
                   << "TRDlikeEln="          << likeEleTRDn
+                 << "trdtruncmean1="       << trdtruncmean1
+                  << "trdtruncmean2="       << trdtruncmean2
                   << "dcaR="                << b[0]
                   << "dcaZ="                << b[1]
                   << "dca="                 << dca
index 07f578c..47a7907 100644 (file)
@@ -25,6 +25,7 @@ class AliHFEcuts;
 class AliHFEsignalCuts;
 class TString;
 class TTreeSRedirector;
+class AliHFEpidTRD;
 
 class AliHFEdebugTreeTask : public AliAnalysisTaskSE{
   public:
@@ -48,6 +49,7 @@ class AliHFEdebugTreeTask : public AliAnalysisTaskSE{
     
     AliHFEcuts *fTrackCuts;           // Track
     AliHFEsignalCuts *fSignalCuts;    // Signal Cuts
+    AliHFEpidTRD *fTRDpid;            // TRD PID
     Int_t fNclustersTPC;              // Min Number of clusters in TPC
     Int_t fNclustersTPCPID;           // Min Number of clusters for TPC PID
     Int_t fNclustersITS;              // Min Number of clusters in ITS
index 46ea73a..bb34fee 100644 (file)
@@ -158,15 +158,15 @@ void AliHFEemcalPIDqa::Initialize(){
   const Double_t kTPCSigMax = 140.;
 
   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, Centrality, select)
-  Int_t nBins[12] = {AliPID::kSPECIES + 1, 500, 500,          400, 630,   200,   600,  300, 125, kCentralityBins, 6, 2};
+  Int_t nBins[12] = {AliPID::kSPECIES + 1, 500, 500,          400, 630,   200,   600,  300, 100, kCentralityBins, 6, 2};
   Double_t min[12] = {-1,               kMinP, kMinP,  kTPCSigMim,  0.,  -1.0,  -8.0,    0,   0,            0, -3.0, 0.};
-  Double_t max[12] = {AliPID::kSPECIES, kMaxP, kMaxP,  kTPCSigMax, 6.3,   1.0,   4.0,  3.0, 0.5,           11., 3.0, 2.};
+  Double_t max[12] = {AliPID::kSPECIES, kMaxP, kMaxP,  kTPCSigMax, 6.3,   1.0,   4.0,  3.0, 0.1,           11., 3.0, 2.};
   fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; phi ; eta ; nSig ; E/p ; Rmatch ;Centrality; PID Step; ", 12, nBins, min, max);
 
   //2nd histogram: EMCAL signal - E/p 
-  Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 125, 600, 2};
+  Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 100, 600, 2};
   Double_t min2[7] = {-1, kMinP, kMinP, 0,  0, -8., 0.};
-  Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5,  0.5, 4., 2.};
+  Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5,  0.1, 4., 2.};
   fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; TPCnsigma; PID Step", 7, nBins2, min2, max2);
     
 }
index 5ee7b15..00c1477 100644 (file)
@@ -82,19 +82,19 @@ class AliHFEextraCuts: public AliCFCutBase{
     Bool_t GetCheckITSstatus() const { return fCheck; };
     Int_t GetDebugLevel() const { return fDebugLevel; };
     void GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy); // temporary moved from protected to publich for IP QA 
+    Int_t GetITSstatus(AliVTrack *track, Int_t layer);
+    Bool_t CheckITSstatus(Int_t itsStatus) const;
     
   protected:
     virtual void AddQAHistograms(TList *qaList);
     Bool_t CheckRecCuts(AliVTrack *track);
     Bool_t CheckMCCuts(AliVParticle * /*track*/) const;
-    Bool_t CheckITSstatus(Int_t itsStatus) const;
     void FillQAhistosRec(AliVTrack *track, UInt_t when);
 //     void FillQAhistosMC(AliMCParticle *track, UInt_t when);
     void FillCutCorrelation(ULong64_t survivedCut);
     void PrintBitMap(Int_t bitmap);
     
     // Getter Functions for ESD/AOD compatible mode
-    Int_t GetITSstatus(AliVTrack *track, Int_t layer);
     UInt_t GetTPCncls(AliVTrack *track);
     UInt_t GetTPCnclusdEdx(AliVTrack *track);
     Bool_t GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track);
index 418eccc..6ddce29 100644 (file)
@@ -32,6 +32,7 @@
 #include <TH2F.h>
 #include <TList.h>
 #include <TParticle.h>
+#include "TTreeStream.h"
 
 #include <AliLog.h>
 #include <AliMCEvent.h>
@@ -56,6 +57,16 @@ ClassImp(AliHFEmcQA)
         ,fCentrality(0)
         ,fIsPbPb(kFALSE)
         ,fIsppMultiBin(kFALSE)
+        ,fContainerStep(0)
+        ,fIsDebugStreamerON(kFALSE)
+        ,fRecPt(-999)
+        ,fRecEta(-999)
+        ,fRecPhi(-999)
+        ,fLyrhit(0)
+        ,fLyrstat(0)
+        ,fHfeImpactR(-999)
+        ,fHfeImpactnsigmaR(-999)
+        ,fTreeStream(NULL)
 {
         // Default constructor
   for(Int_t mom = 0; mom < 9; mom++){
@@ -83,6 +94,16 @@ TObject(p)
         ,fCentrality(0)
         ,fIsPbPb(kFALSE)
         ,fIsppMultiBin(kFALSE)
+        ,fContainerStep(0)
+        ,fIsDebugStreamerON(kFALSE)
+        ,fRecPt(-999)
+        ,fRecEta(-999)
+        ,fRecPhi(-999)
+        ,fLyrhit(0)
+        ,fLyrstat(0)
+        ,fHfeImpactR(0)
+        ,fHfeImpactnsigmaR(0)
+        ,fTreeStream(NULL)
 {
         // Copy constructor
   for(Int_t mom = 0; mom < 9; mom++){
@@ -113,6 +134,7 @@ AliHFEmcQA::~AliHFEmcQA()
 {
         // Destructor
 
+        if(fTreeStream && fIsDebugStreamerON) delete fTreeStream;
         AliInfo("Analysis Done.");
 }
 
@@ -147,15 +169,9 @@ void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
   fQAhistos = qaList;
   fQAhistos->SetName("MCqa");
 
-  CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
-  CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kCharm);               // create histograms for charm
+  CreateHistograms(AliHFEmcQA::kBeauty);              // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers);              // create histograms for beauty
  
 // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles
   const Int_t nbspecies = 9;
@@ -229,10 +245,14 @@ void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
 
   fQAhistos->Add(fMCQACollection->GetList());
 
+  if(!fTreeStream && fIsDebugStreamerON){
+   fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName()));
+  }
+
 }
   
 //__________________________________________
-void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) 
+void AliHFEmcQA::CreateHistograms(const Int_t kquark) 
 {
   // create histograms
 
@@ -266,6 +286,11 @@ void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
     kqTypeLabel[kMisID-4]="miside";
   }
 
+  TString kqEtaRangeLabel[fgkEtaRanges];
+  kqEtaRangeLabel[0] = "mcqa_";
+  kqEtaRangeLabel[1] = "mcqa_barrel_";
+  kqEtaRangeLabel[2] = "mcqa_unitY_";
+
   const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning
   const Int_t nptbinning1 = 35;
   Int_t iBin[2];
@@ -285,64 +310,62 @@ void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt)
 
   TString hname; 
   if(kquark == kOthers){
+   for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
     for (Int_t iqType = 0; iqType < 4; iqType++ ){
-       hname = hnopt+"Pt_"+kqTypeLabel[iqType];
-       //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL
+       hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
        fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25);
-       hname = hnopt+"Y_"+kqTypeLabel[iqType];
+       hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
        fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
-       hname = hnopt+"Eta_"+kqTypeLabel[iqType];
+       hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
        fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
        // Fill List
        if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
     }
-    return;
+   }
+   return;
   }
-  for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
+  for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
+   for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){
      if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron
-     hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
+     hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType];
      fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
-     hname = hnopt+"Pt_"+kqTypeLabel[iqType];
+     hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType];
      fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[1],kPtbinning1); // new binning
-     //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); // old log binning
-     //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.,30.); // mj to compare with FONLL
-     //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL
-     hname = hnopt+"Y_"+kqTypeLabel[iqType];
+     hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType];
      fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
-     hname = hnopt+"Eta_"+kqTypeLabel[iqType];
+     hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType];
      fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5);
      // Fill List
      if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos);
+   }
   }
 
-  if (icut == 0){ 
-    hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
-    fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5);
-    hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
-    fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
+  for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){
+    hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); 
+    hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
+    hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
+    hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1);
+  
+    hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
+    hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
+    hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
+    hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark];
+    fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
+    if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
   }
-  hname = hnopt+"PtCorr_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
-  hname = hnopt+"PtCorrDp_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
-  hname = hnopt+"PtCorrD0_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
-  hname = hnopt+"PtCorrDrest_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning
-  //fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]);
+
+  hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark];
+  fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5);
+  hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark];
+  fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
   
-  hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1);
-  hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1);
-  hname = hnopt+"eDistance_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
-  hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark];
-  fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2);
-  if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos);
 }
 
 //__________________________________________
@@ -791,7 +814,7 @@ void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark)
 }
 
 //__________________________________________
-void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut) 
+void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed) 
 {
     // decay electron kinematics
     
@@ -807,14 +830,27 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
       return;
     }
 
+    Double_t eabsEta = TMath::Abs(mcpart->Eta());
+    Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
+
     if(kquark==kOthers){
       Int_t esource = -1;
       if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4;
       else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6
       if(esource==0|| esource==1 || esource==2 || esource==3){
-        fHist[iq][esource][icut].fPt->Fill(mcpart->Pt());
-        fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-        fHist[iq][esource][icut].fEta->Fill(mcpart->Eta());
+        fHist[iq][esource][0].fPt->Fill(mcpart->Pt());
+        fHist[iq][esource][0].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+        fHist[iq][esource][0].fEta->Fill(mcpart->Eta());
+        if(eabsEta<0.9){
+          fHist[iq][esource][1].fPt->Fill(mcpart->Pt());
+          fHist[iq][esource][1].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+          fHist[iq][esource][1].fEta->Fill(mcpart->Eta());
+        }
+        if(eabsY<0.5){
+          fHist[iq][esource][2].fPt->Fill(mcpart->Pt());
+          fHist[iq][esource][2].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+          fHist[iq][esource][2].fEta->Fill(mcpart->Eta());
+        }
         return; 
       }
       else {
@@ -886,22 +922,62 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
 
                   if (kquark == kCharm) return;
                   // fill electron kinematics
-                  fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-                  fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
-                  fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-                  fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
+                  fHist[iq][kElectron2nd][0].fPdgCode->Fill(mcpart->GetPdgCode());
+                  fHist[iq][kElectron2nd][0].fPt->Fill(mcpart->Pt());
+                  fHist[iq][kElectron2nd][0].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                  fHist[iq][kElectron2nd][0].fEta->Fill(mcpart->Eta());
 
                   // fill mother hadron kinematics
-                  fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); 
-                  fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
-                  fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
-                  fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
-                 
+                  fHist[iq][kDeHadron][0].fPdgCode->Fill(grandMaPDG); 
+                  fHist[iq][kDeHadron][0].fPt->Fill(grandMa->Pt());
+                  fHist[iq][kDeHadron][0].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+                  fHist[iq][kDeHadron][0].fEta->Fill(grandMa->Eta());
+
+                  if(eabsEta<0.9){
+                    fHist[iq][kElectron2nd][1].fPdgCode->Fill(mcpart->GetPdgCode());
+                    fHist[iq][kElectron2nd][1].fPt->Fill(mcpart->Pt());
+                    fHist[iq][kElectron2nd][1].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                    fHist[iq][kElectron2nd][1].fEta->Fill(mcpart->Eta());
+
+                    // fill mother hadron kinematics
+                    fHist[iq][kDeHadron][1].fPdgCode->Fill(grandMaPDG); 
+                    fHist[iq][kDeHadron][1].fPt->Fill(grandMa->Pt());
+                    fHist[iq][kDeHadron][1].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+                    fHist[iq][kDeHadron][1].fEta->Fill(grandMa->Eta());
+                  }
+
+                  if(eabsY<0.5){
+                    fHist[iq][kElectron2nd][2].fPdgCode->Fill(mcpart->GetPdgCode());
+                    fHist[iq][kElectron2nd][2].fPt->Fill(mcpart->Pt());
+                    fHist[iq][kElectron2nd][2].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                    fHist[iq][kElectron2nd][2].fEta->Fill(mcpart->Eta());
+
+                    // fill mother hadron kinematics
+                    fHist[iq][kDeHadron][2].fPdgCode->Fill(grandMaPDG); 
+                    fHist[iq][kDeHadron][2].fPt->Fill(grandMa->Pt());
+                    fHist[iq][kDeHadron][2].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+                    fHist[iq][kDeHadron][2].fEta->Fill(grandMa->Eta());
+                  }
+
                   // ratio between pT of electron and pT of mother B hadron 
-                  if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                  if(grandMa->Pt()) {
+                    fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                    if(eabsEta<0.9){
+                      fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                    }
+                    if(eabsY<0.5){
+                      fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
+                    }
+                  }
 
                   // distance between electron production point and primary vertex
-                  fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  if(eabsEta<0.9){
+                    fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  }
+                  if(eabsY<0.5){
+                    fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy);
+                  }
                   return;
                 }
              } 
@@ -913,51 +989,103 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde
          for (Int_t i=0; i<fNparents; i++){
             if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 
-//mj weighting to consider measured spectra!!!
-              Double_t mpt=partMotherCopy->Pt();
-              Double_t wfactor=2*(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225)); // 2* considering in pythia I used particle + antiparticle differently from the measurement
-              //Double_t wfactor=(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225));
-              // fill electron kinematics
-              if(iq==0){
-                fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode(),wfactor);
-                fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt(),wfactor);
-                fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart),wfactor);
-                fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta(),wfactor);  
-
-                fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-                fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
-                fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-                fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());  
-              } 
-              else{
-                fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-                fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
-                fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-                fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());  
-              }
-//--------------
+              fHist[iq][kElectron][0].fPdgCode->Fill(mcpart->GetPdgCode());
+              fHist[iq][kElectron][0].fPt->Fill(mcpart->Pt());
+              fHist[iq][kElectron][0].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+              fHist[iq][kElectron][0].fEta->Fill(mcpart->Eta());  
+
               // fill mother hadron kinematics
-              fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); 
-              fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
-              fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
-              fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
+              fHist[iq][keHadron][0].fPdgCode->Fill(maPdgcodeCopy); 
+              fHist[iq][keHadron][0].fPt->Fill(partMotherCopy->Pt());
+              fHist[iq][keHadron][0].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+              fHist[iq][keHadron][0].fEta->Fill(partMotherCopy->Eta());
+
+              if(eabsEta<0.9){
+                fHist[iq][kElectron][1].fPdgCode->Fill(mcpart->GetPdgCode());
+                fHist[iq][kElectron][1].fPt->Fill(mcpart->Pt());
+                fHist[iq][kElectron][1].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                fHist[iq][kElectron][1].fEta->Fill(mcpart->Eta());  
+
+                // fill mother hadron kinematics
+                fHist[iq][keHadron][1].fPdgCode->Fill(maPdgcodeCopy); 
+                fHist[iq][keHadron][1].fPt->Fill(partMotherCopy->Pt());
+                fHist[iq][keHadron][1].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+                fHist[iq][keHadron][1].fEta->Fill(partMotherCopy->Eta());
+              }
+
+              if(eabsY<0.5){
+                fHist[iq][kElectron][2].fPdgCode->Fill(mcpart->GetPdgCode());
+                fHist[iq][kElectron][2].fPt->Fill(mcpart->Pt());
+                fHist[iq][kElectron][2].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+                fHist[iq][kElectron][2].fEta->Fill(mcpart->Eta());  
+
+                // fill mother hadron kinematics
+                fHist[iq][keHadron][2].fPdgCode->Fill(maPdgcodeCopy); 
+                fHist[iq][keHadron][2].fPt->Fill(partMotherCopy->Pt());
+                fHist[iq][keHadron][2].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+                fHist[iq][keHadron][2].fEta->Fill(partMotherCopy->Eta());
+              }
 
               // ratio between pT of electron and pT of mother B or direct D hadron 
-              if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
-              fHistComm[iq][icut].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
-              if(TMath::Abs(partMotherCopy->GetPdgCode())==411) fHistComm[iq][icut].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
-              else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) fHistComm[iq][icut].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
-              else fHistComm[iq][icut].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              if(partMotherCopy->Pt()) {
+                 fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
+                 if(eabsEta<0.9){
+                   fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
+                 }
+                 if(eabsY<0.5){
+                   fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
+                 }
+              }
+              fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              if(eabsEta<0.9){
+                fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              }
+              if(eabsY<0.5){
+                fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt());
+              }
+              if(TMath::Abs(partMotherCopy->GetPdgCode())==411) {
+                fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                if(eabsEta<0.9){
+                  fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+                if(eabsY<0.5){
+                  fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+              }
+              else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) {
+                fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                if(eabsEta<0.9){
+                  fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+                if(eabsY<0.5){
+                  fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+              }
+              else {
+                fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                if(eabsEta<0.9){
+                  fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+                if(eabsY<0.5){
+                  fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt());
+                }
+              }
 
               // distance between electron production point and primary vertex
-              fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              if(eabsEta<0.9){
+                fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              }
+              if(eabsY<0.5){
+                fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
+              }
             }
          }
     } // end of if
 }
 
 //____________________________________________________________________
-void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut)
+void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed)
 {
   // decay electron kinematics
 
@@ -976,6 +1104,9 @@ void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I
 
   if ( abs(mcpart->GetPdgCode()) != kdecayed ) return;
 
+  Double_t eabsEta = TMath::Abs(mcpart->Eta());
+  Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart));
+
   // mother
   Int_t iLabel = mcpart->GetMother();
   if (iLabel<0){
@@ -1019,16 +1150,43 @@ void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I
 
             if (kquark == kCharm) return;
             // fill electron kinematics
-            fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-            fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt());
-            fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-            fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta());
+            fHist[iq][kElectron2nd][0].fPdgCode->Fill(mcpart->GetPdgCode());
+            fHist[iq][kElectron2nd][0].fPt->Fill(mcpart->Pt());
+            fHist[iq][kElectron2nd][0].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+            fHist[iq][kElectron2nd][0].fEta->Fill(mcpart->Eta());
 
             // fill mother hadron kinematics
-            fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG);
-            fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt());
-            fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa));
-            fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta());
+            fHist[iq][kDeHadron][0].fPdgCode->Fill(grandMaPDG);
+            fHist[iq][kDeHadron][0].fPt->Fill(grandMa->Pt());
+            fHist[iq][kDeHadron][0].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+            fHist[iq][kDeHadron][0].fEta->Fill(grandMa->Eta());
+
+            if(eabsEta<0.9){
+              // fill electron kinematics
+              fHist[iq][kElectron2nd][1].fPdgCode->Fill(mcpart->GetPdgCode());
+              fHist[iq][kElectron2nd][1].fPt->Fill(mcpart->Pt());
+              fHist[iq][kElectron2nd][1].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+              fHist[iq][kElectron2nd][1].fEta->Fill(mcpart->Eta());
+
+              // fill mother hadron kinematics
+              fHist[iq][kDeHadron][1].fPdgCode->Fill(grandMaPDG);
+              fHist[iq][kDeHadron][1].fPt->Fill(grandMa->Pt());
+              fHist[iq][kDeHadron][1].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+              fHist[iq][kDeHadron][1].fEta->Fill(grandMa->Eta());
+            }
+            if(eabsY<0.5){
+              // fill electron kinematics
+              fHist[iq][kElectron2nd][2].fPdgCode->Fill(mcpart->GetPdgCode());
+              fHist[iq][kElectron2nd][2].fPt->Fill(mcpart->Pt());
+              fHist[iq][kElectron2nd][2].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+              fHist[iq][kElectron2nd][2].fEta->Fill(mcpart->Eta());
+
+              // fill mother hadron kinematics
+              fHist[iq][kDeHadron][2].fPdgCode->Fill(grandMaPDG);
+              fHist[iq][kDeHadron][2].fPt->Fill(grandMa->Pt());
+              fHist[iq][kDeHadron][2].fY->Fill(AliHFEtools::GetRapidity(grandMa));
+              fHist[iq][kDeHadron][2].fEta->Fill(grandMa->Eta());
+            }
 
             return;
           }
@@ -1042,16 +1200,43 @@ void  AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I
        if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){
 
          // fill electron kinematics
-         fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode());
-         fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt());
-         fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart));
-         fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta());
+         fHist[iq][kElectron][0].fPdgCode->Fill(mcpart->GetPdgCode());
+         fHist[iq][kElectron][0].fPt->Fill(mcpart->Pt());
+         fHist[iq][kElectron][0].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+         fHist[iq][kElectron][0].fEta->Fill(mcpart->Eta());
 
          // fill mother hadron kinematics
-         fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy);
-         fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt());
-         fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
-         fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta());
+         fHist[iq][keHadron][0].fPdgCode->Fill(maPdgcodeCopy);
+         fHist[iq][keHadron][0].fPt->Fill(partMotherCopy->Pt());
+         fHist[iq][keHadron][0].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+         fHist[iq][keHadron][0].fEta->Fill(partMotherCopy->Eta());
+
+         if(eabsEta<0.9){
+           // fill electron kinematics
+           fHist[iq][kElectron][1].fPdgCode->Fill(mcpart->GetPdgCode());
+           fHist[iq][kElectron][1].fPt->Fill(mcpart->Pt());
+           fHist[iq][kElectron][1].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+           fHist[iq][kElectron][1].fEta->Fill(mcpart->Eta());
+
+           // fill mother hadron kinematics
+           fHist[iq][keHadron][1].fPdgCode->Fill(maPdgcodeCopy);
+           fHist[iq][keHadron][1].fPt->Fill(partMotherCopy->Pt());
+           fHist[iq][keHadron][1].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+           fHist[iq][keHadron][1].fEta->Fill(partMotherCopy->Eta());
+         }
+         if(eabsY<0.5){
+           // fill electron kinematics
+           fHist[iq][kElectron][2].fPdgCode->Fill(mcpart->GetPdgCode());
+           fHist[iq][kElectron][2].fPt->Fill(mcpart->Pt());
+           fHist[iq][kElectron][2].fY->Fill(AliHFEtools::GetRapidity(mcpart));
+           fHist[iq][kElectron][2].fEta->Fill(mcpart->Eta());
+
+           // fill mother hadron kinematics
+           fHist[iq][keHadron][2].fPdgCode->Fill(maPdgcodeCopy);
+           fHist[iq][keHadron][2].fPt->Fill(partMotherCopy->Pt());
+           fHist[iq][keHadron][2].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy));
+           fHist[iq][keHadron][2].fEta->Fill(partMotherCopy->Eta());
+         }
 
        }
     }
@@ -1486,6 +1671,21 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
   else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime
   else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5;         //rho
 
+  Double_t datamc[24]={-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999, -999, -999, -999, -999, -999, -999, -999, -999};
+  Double_t xr[3]={-999,-999,-999};
+  datamc[0] = mesonID;
+  datamc[17] = mctrack->Pt(); //electron pt
+  datamc[18] = mctrack->Eta(); //electron eta
+
+  mctrack->XvYvZv(xr);
+  datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+  datamc[10] = xr[2];
+
+  TParticle *mcpart = mctrack->Particle();
+  if(mcpart){
+    datamc[14] = mcpart->GetUniqueID();
+  }
+
   if(!(mArr<0)){
      if(mesonID>=kGammaPi0) {  // conversion electron, be careful with the enum odering 
         Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label
@@ -1494,6 +1694,37 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
           if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
             mesonPt = mctrackmother->Pt(); //meson pt
             bgcategory = 1.;
+            datamc[1] = bgcategory;
+            datamc[2] = mesonPt;
+            mctrackmother->XvYvZv(xr);
+            datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+            datamc[12] = xr[2];
+
+            mcpart = mctrackmother->Particle();
+            if(mcpart){
+              datamc[15] = mcpart->GetUniqueID();
+            }
+            if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+              bgcategory = 2.;
+              datamc[1] = bgcategory;
+              //printf("I should be gamma meson = %d  mesonlabel= %d  NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries()); 
+              glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+              if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                datamc[3]=mctrackmother->PdgCode();
+                datamc[4]=mctrackmother->Pt();
+                if(TMath::Abs(mctrackmother->PdgCode())==310){
+                  bgcategory = 3.;
+                  glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+                  if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                    datamc[5]=mctrackmother->PdgCode();
+                    glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+                    if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                       datamc[6]=mctrackmother->PdgCode();
+                    }
+                  }
+                }
+              }
+            } 
           } 
         }
      }
@@ -1502,9 +1733,40 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
           mesonPt=mctrackmother->Pt(); //meson pt
           bgcategory = -1.;
+          datamc[1] = bgcategory;
+          datamc[2] = mesonPt;
+          mctrackmother->XvYvZv(xr);
+          datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+          datamc[12] = xr[2];
+
+          mcpart = mctrackmother->Particle();
+          if(mcpart){
+            datamc[15] = mcpart->GetUniqueID();
+          }
+          if(glabel>fMCEvent->GetNumberOfPrimaries()) {
+            bgcategory = -2.;
+            datamc[1] = bgcategory;
+            glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+            if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+              datamc[3]=mctrackmother->PdgCode();
+              datamc[4]=mctrackmother->Pt();
+              if(TMath::Abs(mctrackmother->PdgCode())==310){
+               bgcategory = -3.;
+               glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother
+               if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                 datamc[5]=mctrackmother->PdgCode();
+                 glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother
+                 if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
+                   datamc[6]=mctrackmother->PdgCode();
+                 }
+               }
+              }  
+            }
+          }
         }
      }
 
+
      if(fIsPbPb){
        if(fCentrality < 0)return 0.;
        weightElecBg=fElecBackgroundFactor[iBgLevel][fCentrality][mArr][kBgPtBins-1];                        
@@ -1525,7 +1787,54 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve
        }
      }    
   }
-  return bgcategory*weightElecBg;
+
+  datamc[13] = weightElecBg;
+  datamc[16] = Double_t(fContainerStep);
+
+  datamc[7] = fHfeImpactR;
+  datamc[8] = fHfeImpactnsigmaR;
+
+  datamc[19] = fRecPt;
+  datamc[20] = fRecEta;
+  datamc[21] = fRecPhi;
+  datamc[22] = fLyrhit;
+  datamc[23] = fLyrstat;
+
+  if(fIsDebugStreamerON && fTreeStream){
+   if(!iBgLevel){
+    (*fTreeStream)<<"nonhfeQA"<<
+        "mesonID="<<datamc[0]<<
+        "bgcategory="<<datamc[1]<<
+        "mesonPt="<<datamc[2]<<
+        "mesonMomPdg="<<datamc[3]<<
+        "mesonMomPt="<<datamc[4]<<
+        "mesonGMomPdg="<<datamc[5]<<
+        "mesonGGMomPdg="<<datamc[6]<<
+        "eIPAbs="<<datamc[7]<<
+        "eIPSig="<<datamc[8]<<
+        "eR="<<datamc[9]<<
+        "eZ="<<datamc[10]<<
+        "mesonR="<<datamc[11]<<
+        "mesonZ="<<datamc[12]<<
+        "weightElecBg="<<datamc[13]<< 
+        "eUniqID="<<datamc[14]<<
+        "mesonUniqID="<<datamc[15]<<
+        "containerStep="<<datamc[16]<<
+        "emcpt="<<datamc[17]<<
+        "emceta="<<datamc[18]<<
+        "erecpt="<<datamc[19]<<
+        "ereceta="<<datamc[20]<<
+        "erecphi="<<datamc[21]<< 
+        "itshit="<<datamc[22]<<
+        "itsstat="<<datamc[23]
+        << "\n";
+   }
+  }
+
+  Double_t returnval = bgcategory*weightElecBg;
+  if(TMath::Abs(bgcategory)>1) returnval = bgcategory/2.;
+
+  return returnval;
 }
 
 //__________________________________________
index b8cf76d..55d2c7a 100644 (file)
@@ -40,6 +40,7 @@ class AliGenEventHeader;
 class AliMCParticle;
 class AliAODMCParticle;
 class AliHFEcollection;
+class TTreeSRedirector;
 
 //________________________________________________________________
 class AliHFEmcQA: public TObject {
@@ -68,7 +69,7 @@ class AliHFEmcQA: public TObject {
     TList *GetList() const { return fQAhistos; };
     void PostAnalyze() const;
     void CreatDefaultHistograms(TList * const qaList); // create default histograms  
-    void CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt=""); // create histograms for mc qa analysis
+    void CreateHistograms(const Int_t kquark); // create histograms for mc qa analysis
     void SetMCEvent(AliMCEvent* const mcEvent){fMCEvent = mcEvent;} 
     void SetGenEventHeader(AliGenEventHeader* const mcHeader){fMCHeader=mcHeader;} // set stack pointer
     void SetMCArray(TClonesArray* const mcarry){fMCArray=mcarry;} // set mcarray pointer
@@ -76,16 +77,21 @@ class AliHFEmcQA: public TObject {
 
     void GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark); // get heavy quark kinematics distribution
     void GetHadronKine(TParticle *part, const Int_t kquark); // get heavy hadron kinematics distribution
-    void GetDecayedKine(TParticle *part, const Int_t kquark, const Int_t kdecayed, Int_t icut); // get decay electron kinematics distribution
-    void GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut); // get decay electron kinematics for AOD 
+    void GetDecayedKine(TParticle *part, const Int_t kquark, const Int_t kdecayed); // get decay electron kinematics distribution
+    void GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed); // get decay electron kinematics for AOD 
     void GetMesonKine(); // get meson and its decay electron pt spectra
     void EndOfEventAna(const Int_t kquark); // run analysis which should be done at the end of the event loop
     Int_t GetSource(const TParticle * const mcpart); // return source id 
     Int_t GetElecSource(TParticle * const mcpart); // return electron source id 
     Int_t GetSource(const AliAODMCParticle * const mcpart); // return electron source id for AOD
     Double_t GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel); // return best/lower/upper weighting factor for electron's mother meson
+    void EnableDebugStreamer() { fIsDebugStreamerON = kTRUE;};
 
     void SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit);
+    void SetContainerStep(Int_t containerStep) { fContainerStep = containerStep;};
+    void SetHFEImpactParameters(Double_t hfeimpactR, Double_t hfeimpactnsigmaR) {fHfeImpactR = hfeimpactR; fHfeImpactnsigmaR = hfeimpactnsigmaR; };
+    void SetTrkKine(Double_t pt, Double_t eta, Double_t phi) {fRecPt = pt; fRecEta = eta; fRecPhi = phi;};
+    void SetITSInfo(Double_t ilyrhit, Double_t ilyrstat) { fLyrhit = ilyrhit; fLyrstat = ilyrstat;};
 
     void SetCentrality(Int_t centrality) { fCentrality = centrality; };
     void SetPbPb() { fIsPbPb = kTRUE; };
@@ -110,6 +116,7 @@ class AliHFEmcQA: public TObject {
     static const Int_t fgkMaxGener=10; // ancester level wanted to be checked 
     static const Int_t fgkMaxIter=100; // number of iteration to find out matching particle 
     static const Int_t fgkqType=7; // number of particle type to be checked
+    static const Int_t fgkEtaRanges=3; // cuts for different eta ranges
 
     struct AliHists{
       TH1F *fPdgCode; // histogram to store particle pdg code
@@ -206,6 +213,19 @@ private:
     Int_t              fCentrality;  // Centrality
     Bool_t             fIsPbPb;        // Analysis Type: pp or PbPb
     Bool_t             fIsppMultiBin;  // pp multiplicity bin analysis
+    Int_t              fContainerStep; // step the weighting factor called
+    Bool_t             fIsDebugStreamerON; // check if the debugstreamer is on
+
+    Double_t           fRecPt;  //reconstructed pt
+    Double_t           fRecEta; //reconstructed eta
+    Double_t           fRecPhi; //reconstructed phi
+    Double_t           fLyrhit; //its layer hit
+    Double_t           fLyrstat; // its layer status
+
+    Double_t fHfeImpactR;              // absolute impact parameter R
+    Double_t fHfeImpactnsigmaR;        // absolute impact parameter sigma R 
+    TTreeSRedirector *fTreeStream;        //! TreeStream
 
   ClassDef(AliHFEmcQA,1);
 };
index 8453be9..027c2bd 100644 (file)
@@ -159,8 +159,8 @@ Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanage
         feopMaxCut = CalEopCutMax(track->GetRecTrack(),1); 
        }
 
-  //printf("eop cuts org; %g; %g \n",feopMim,feopMax);
-  //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);
+ //printf("eop cuts org; %g; %g \n",feopMim,feopMax);
+ //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut);
 
   Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
   AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
@@ -198,20 +198,23 @@ Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t fla
 
   if(flageop<0.5)
     { //<--- new para for updated non-liniarity
-     double meanP[3] = {1.02902,1.22293,1.67287e-01}; 
-     double sigP[3] = {5.36355e-02,-3.25999e-02,4.42442e-01}; 
+     double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02}; 
+     double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04}; 
      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
      maxCut = mean+3.0*sig; 
     }
   else
     {
-     double meanP[3] = {0.99,1.299,0.35}; 
-     double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; 
+     double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01}; 
+     double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01}; 
      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
      maxCut = mean+3.0*sig; 
     }
+
+  //printf("eop cuts org; %g; %g ; %g\n",feopMim,feopMax,pt);
+
      return maxCut;
 }
 
@@ -226,22 +229,20 @@ Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t fla
   if(flageop<0.5) // real
      { 
       //printf("real"); 
-     //double meanP[3] = {1.056,1.169,0.1339}; 
-     //double sigP[3] = {3.37e-05,1.298e-04,1.403e-04}; 
-     double meanP[3] = {1.02902,1.22293,1.67287e-01}; 
-     double sigP[3] = {5.36355e-02,-3.25999e-02,4.42442e-01}; 
+     double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02}; 
+     double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04}; 
      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
-     mimCut = mean+0.0*sig;
+     mimCut = mean-3.0*sig;
     }
    else // MC
     { 
       //printf("MC"); 
-     double meanP[3] = {0.99,1.299,0.35}; 
-     double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; 
+     double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01}; 
+     double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01}; 
      double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); 
      double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); 
-     mimCut = mean+0.0*sig;
+     mimCut = mean-3.0*sig;
     }
      return mimCut;
 }
index 9a96f69..f1f8234 100644 (file)
@@ -1309,7 +1309,7 @@ Int_t AliHFEpidQA::GetCentrality(AliVEvent*  const fInputEvent){
 }
 
 //___________________________________________________
-Int_t AliHFEpidQA::GetMultiplicityITS(AliVEvent*  const fInputEvent)
+Int_t AliHFEpidQA::GetMultiplicityITS(AliVEvent*  const fInputEvent) const
 {
   //
   // Definition of the Multiplicity according to the JPSI group (F. Kramer)
index 721b5d9..fcd242f 100644 (file)
@@ -80,7 +80,7 @@ class AliHFEpidQA : public TObject{
     AliHFEtrdPIDqa *GetTRDQA() const { return fTRDpidQA; }
 
     Int_t GetCentrality(AliVEvent*  const fInputEvent);
-    Int_t GetMultiplicityITS(AliVEvent*  const fInputEvent);
+    Int_t GetMultiplicityITS(AliVEvent*  const fInputEvent) const;
     void SetPbPb() { fIsPbPb = kTRUE; };
     void SetPP() { fIsPbPb = kFALSE; };
     void SetPPMultiBin() { fIsppMultiBin = kFALSE; };
index bddaf85..02990c7 100644 (file)
@@ -320,6 +320,7 @@ Int_t AliHFEsignalCuts::GetElecSource(const AliVParticle * const track) const {
   TClass *tracktype;
   const AliVParticle *mctrack = NULL;
   TParticle *mcpart = NULL;
+  //AliMCParticle *esdmcmother = NULL;
   if((tracktype = track->IsA()) == AliESDtrack::Class() || tracktype == AliAODTrack::Class()){
     mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel()));
   } else  mctrack = track;
@@ -331,6 +332,23 @@ Int_t AliHFEsignalCuts::GetElecSource(const AliVParticle * const track) const {
     if(esdmc){
       mcpart = esdmc->Particle();
       eSource=fMCQA->GetElecSource(mcpart);
+/* // considering secondary pions
+      if(eSource>=AliHFEmcQA::kGammaPi0) {  // conversion electron, be careful with the enum odering 
+        Int_t glabel=TMath::Abs(esdmc->GetMother()); // gamma label
+        if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
+          glabel=TMath::Abs(esdmcmother->GetMother()); // gamma's mother's label
+          if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
+            if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kElse;
+          }
+        }
+      }
+      else if(eSource==AliHFEmcQA::kPi0 || (eSource>=AliHFEmcQA::kEta && eSource<=AliHFEmcQA::kRho0) ){ // nonHFE except for the conversion electron
+        Int_t glabel=TMath::Abs(esdmc->GetMother());
+        if((esdmcmother= dynamic_cast<AliMCParticle *>(fMC->GetTrack(glabel)))){
+          if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kElse;
+        }
+      }
+*/
     }
   } else {
     return -1;
index 2c59c9d..61af42b 100644 (file)
@@ -1,3 +1,4 @@
+
 /**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
@@ -76,11 +77,11 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
   fIPanaCharmBgSubtract(kFALSE),
   fIPanaConversionBgSubtract(kFALSE),
   fIPanaNonHFEBgSubtract(kFALSE),
-  fNonHFEbgMethod2(kFALSE),
+  fIPParameterizedEff(kFALSE),
   fNonHFEsyst(0),
+  fBeauty2ndMethod(kFALSE),
   fNbDimensions(1),
   fNMCEvents(0),
-  fNMCbgEvents(0),
   fStepMC(-1),
   fStepTrue(-1),
   fStepData(-1),
@@ -91,8 +92,8 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
   fChargeChoosen(kAllCharge),
   fNCentralityBinAtTheEnd(0),
   fHadronEffbyIPcut(NULL),
-  fConversionEff(NULL),
-  fNonHFEEff(NULL),
+  fBSpectrum2ndMethod(NULL),
+  fkBeauty2ndMethodfilename(""),
   fBeamType(0),
   fDebugLevel(0),
   fWriteToFile(kFALSE)
@@ -102,9 +103,25 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
   //
 
   for(Int_t k = 0; k < 20; k++){
-    fNEvents[k] = 0;
-    fLowBoundaryCentralityBinAtTheEnd[k] = 0;
-    fHighBoundaryCentralityBinAtTheEnd[k] = 0;
+      fNEvents[k] = 0;
+      fNMCbgEvents[k] = 0;
+      fLowBoundaryCentralityBinAtTheEnd[k] = 0;
+      fHighBoundaryCentralityBinAtTheEnd[k] = 0;
+      if(k<kCentrality)
+      {
+          fEfficiencyTOFPIDD[k] = 0;
+          fEfficiencyIPCharmD[k] = 0;     
+          fEfficiencyIPBeautyD[k] = 0;    
+          fEfficiencyIPConversionD[k] = 0;
+          fEfficiencyIPNonhfeD[k] = 0;   
+
+         fConversionEff[k] = 0;
+         fNonHFEEff[k] = 0;
+         fCharmEff[k] = 0;
+         fBeautyEff[k] = 0;
+         fEfficiencyCharmSigD[k] = 0;
+          fEfficiencyBeautySigD[k] = 0;
+      }
   }
   memset(fEtaRange, 0, sizeof(Double_t) * 2);
   memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
@@ -136,9 +153,20 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
   // If no inclusive spectrum, then take only efficiency map for beauty electron
   // and the appropriate correlation matrix
   //
+
+
+  if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
+
   Int_t kNdim = 3;
+  Int_t kNcentr =1;
+  Int_t ptpr =0;
   if(fBeamType==0) kNdim=3;
-  if(fBeamType==1) kNdim=4;
+  if(fBeamType==1)
+  {
+      kNdim=4;
+      kNcentr=11;
+      ptpr=1;
+  }
 
   Int_t dims[kNdim];
   // Get the requested format
@@ -183,7 +211,8 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
     datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
   }
   else{
-    datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
+
+      datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
   }
   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
@@ -219,10 +248,21 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
       for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
        for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
           nonHFEtempContainer =  bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
-         fNonHFESourceContainer[iSource][iLevel] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
-          convtempContainer =  bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
-         fConvSourceContainer[iSource][iLevel] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
-          if((!fConvSourceContainer[iSource][iLevel])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;   
+         convtempContainer =  bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
+         for(Int_t icentr=0;icentr<kNcentr;icentr++)
+         {
+             if(fBeamType==0)
+             {
+                 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
+                 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
+             }
+             if(fBeamType==1)
+             {
+                 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr);
+                 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr);
+             }
+             if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
+         }
        }
       }
     }
@@ -248,17 +288,7 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
    mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
    SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
 
-   source = 2; //conversion
-   mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
-   AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
-   efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
-   fConversionEff = (TH1D*)efficiencyConv->Project(0);
-
-   source = 3; //non HFE except for the conversion
-   mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
-   AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
-   efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step #
-   fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0);
+   SetParameterizedEff(mccontainermc, mccontainermcbg, dims);
 
    if(!fNonHFEsyst){
      AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
@@ -347,6 +377,17 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
 
 
 //____________________________________________________________
+void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){
+  //
+  // get spectrum for beauty 2nd method
+  //
+  //
+    TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
+    fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
+}
+
+
+//____________________________________________________________
 Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
   //
   // Correct the spectrum for efficiency and unfolding
@@ -708,12 +749,15 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
   AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
   TGraphErrors *gNormalizedRawSpectrum = 0x0;
   if(subtractcontamination) {
-    dataspectrumaftersubstraction = SubtractBackground(kTRUE);
-    unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
-    dataGridAfterFirstSteps = dataspectrumaftersubstraction;
-    gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
+      if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
+      else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
+      unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
+      dataGridAfterFirstSteps = dataspectrumaftersubstraction;
+      gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
   }
 
+  printf("after normalize getting IP \n");
+
   /////////////////////////////////////////////////////////////////////////////////////////
   // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
   /////////////////////////////////////////////////////////////////////////////////////////
@@ -724,7 +768,7 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
 
   if(fEfficiencyFunction){
     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
-    dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;  
+    dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
   }
   else if(dataContainerV0){
     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
@@ -732,6 +776,7 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
   }  
  
 
+
   ///////////////
   // Unfold
   //////////////
@@ -754,13 +799,19 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
   /////////////////////
   // Simply correct
   ////////////////////
+
   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
   
 
   //////////
   // Plot
   //////////
+
   if(fDebugLevel > 0.0) {
+
+    Int_t ptpr = 0;
+    if(fBeamType==0) ptpr=0;
+    if(fBeamType==1) ptpr=1;
   
     TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
     ccorrected->Divide(2,1);
@@ -787,8 +838,8 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
     legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
     legcorrected->Draw("same");
     ccorrected->cd(2);
-    TH1D *correctedTH1D = correctedspectrum->Projection(0);
-    TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
+    TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
+    TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
     TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
     ratiocorrected->SetName("ratiocorrected");
     ratiocorrected->SetTitle("");
@@ -799,15 +850,17 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
     ratiocorrected->Draw();
     if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps");
 
-    if(fNonHFEsyst){
-      CalculateNonHFEsyst();
+    if(fBeamType == 0){
+       if(fNonHFEsyst){
+           CalculateNonHFEsyst(0);
+       }
     }
 
-
     // Dump to file if needed
 
     if(fDumpToFile) {
-      
+        // to do centrality dependent
+
       TFile *out;
       out = new TFile("finalSpectrum.root","recreate");
       out->cd();
@@ -824,16 +877,73 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
       alltogetherCorrection->Write();
       //
-      if(unnormalizedRawSpectrum) {
-       unnormalizedRawSpectrum->SetName("beautyAfterIP");
-       unnormalizedRawSpectrum->Write();
-      }
+      unnormalizedRawSpectrum->SetName("beautyAfterIP");
+      unnormalizedRawSpectrum->Write();
       
       if(gNormalizedRawSpectrum){
         gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
         gNormalizedRawSpectrum->Write();
       }
 
+      if(fBeamType==1) {
+
+         TGraphErrors* correctedspectrumDc[kCentrality];
+          TGraphErrors* alltogetherspectrumDc[kCentrality];
+         for(Int_t i=0;i<kCentrality-1;i++)
+         {
+             correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
+             correctedspectrumDc[i] = Normalize(correctedspectrum,i);
+             correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
+             correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
+             alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
+             alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
+             alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
+             alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
+             correctedspectrumDc[i]->Write();
+             alltogetherspectrumDc[i]->Write();
+             
+
+             TH1D *centrcrosscheck = correctedspectrum->Projection(0);
+             centrcrosscheck->SetTitle(Form("centrality_%i",i));
+             centrcrosscheck->SetName(Form("centrality_%i",i));
+             centrcrosscheck->Write();
+
+             TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr);
+             TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr);
+
+             TH1D *centrcrosscheck2 =  (TH1D *) alltogetherCorrection->Project(0);
+             centrcrosscheck2->SetTitle(Form("centrality2_%i",i));
+             centrcrosscheck2->SetName(Form("centrality2_%i",i));
+             centrcrosscheck2->Write();
+
+             TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone();
+             ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1);
+             ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
+             ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
+              ratiocorrectedc->Write();
+
+              fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
+             fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
+              fEfficiencyCharmSigD[i]->Write();
+             fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
+             fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
+              fEfficiencyBeautySigD[i]->Write();
+             fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
+             fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
+              fCharmEff[i]->Write();
+             fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i));
+             fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i));
+              fBeautyEff[i]->Write();
+             fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i));
+             fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i));
+              fConversionEff[i]->Write();
+             fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i));
+             fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i));
+              fNonHFEEff[i]->Write();
+         }
+
+      }
+
       out->Close(); 
       delete out;
     }
@@ -847,7 +957,20 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
   //
   // Apply background subtraction
   //
-  
+
+    Int_t ptpr = 0;
+    Int_t nbins=1;
+    if(fBeamType==0)
+    {
+       ptpr=0;
+        nbins=1;
+    }
+    if(fBeamType==1)
+    {
+       ptpr=1;
+       nbins=2;
+    }
+
   // Raw spectrum
   AliCFContainer *dataContainer = GetContainer(kDataContainer);
   if(!dataContainer){
@@ -872,27 +995,39 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
 
   if(!fInclusiveSpectrum){
     //Background subtraction for IP analysis
-    TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(0);
+    TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
     CorrectFromTheWidth(measuredTH1Draw);
     TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500);
     rawspectra->cd();
-    rawspectra->SetLogx();
     rawspectra->SetLogy();
     TLegend *lRaw = new TLegend(0.65,0.65,0.95,0.95);
     measuredTH1Draw->SetMarkerStyle(20);
     measuredTH1Draw->Draw();
     lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
+    TH1D* htemp;
+    Int_t* bins=new Int_t[1];
     if(fIPanaHadronBgSubtract){
       // Hadron background
-      printf("Hadron background for IP analysis subtracted!\n");
-      Int_t* bins=new Int_t[1];
-      TH1D* htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
-      bins[0]=htemp->GetNbinsX();
-      AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins);
+       printf("Hadron background for IP analysis subtracted!\n");
+       if(fBeamType==0)
+       {
+
+           htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
+           bins[0]=htemp->GetNbinsX();
+       }
+       if(fBeamType==1)
+       {
+           bins=new Int_t[2];
+           htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
+           bins[0]=htemp->GetNbinsX();
+           htemp  = (TH1D *) fHadronEffbyIPcut->Projection(1);
+           bins[1]=htemp->GetNbinsX();
+       }
+      AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
       hbgContainer->SetGrid(fHadronEffbyIPcut);
       backgroundGrid->Multiply(hbgContainer,1);
       // draw raw hadron bg spectra
-      TH1D *hadronbg= (TH1D *) backgroundGrid->Project(0);
+      TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr);
       CorrectFromTheWidth(hadronbg);
       hadronbg->SetMarkerColor(7);
       hadronbg->SetMarkerStyle(20);
@@ -907,7 +1042,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
       printf("Charm background for IP analysis subtracted!\n");
       AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
       // draw charm bg spectra
-      TH1D *charmbg= (TH1D *) charmbgContainer->Project(0);
+      TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr);
       CorrectFromTheWidth(charmbg);
       charmbg->SetMarkerColor(3);
       charmbg->SetMarkerStyle(20);
@@ -919,9 +1054,9 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
     }
     if(fIPanaConversionBgSubtract){
       // Conversion background
-      AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
+      AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground(); 
       // draw conversion bg spectra
-      TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(0);
+      TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr);
       CorrectFromTheWidth(conversionbg);
       conversionbg->SetMarkerColor(4);
       conversionbg->SetMarkerStyle(20);
@@ -936,7 +1071,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
       // NonHFE background
       AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
       // draw Dalitz/dielectron bg spectra
-      TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(0);
+      TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr);
       CorrectFromTheWidth(nonhfebg);
       nonhfebg->SetMarkerColor(6);
       nonhfebg->SetMarkerStyle(20);
@@ -947,7 +1082,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
       spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
       printf("Non HFE background subtraction is preliminary!\n");
     }
-    TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(0);
+    TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
     CorrectFromTheWidth(rawbgsubtracted);
     rawbgsubtracted->SetMarkerStyle(24);
     rawspectra->cd();
@@ -968,16 +1103,16 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
 
   if(fDebugLevel > 0) {
 
-    Int_t ptpr;
-    if(fBeamType==0) ptpr=0;
-    if(fBeamType==1) ptpr=1;
+    Int_t ptprd;
+    if(fBeamType==0) ptprd=0;
+    if(fBeamType==1) ptprd=1;
     
     TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
     cbackgroundsubtraction->Divide(3,1);
     cbackgroundsubtraction->cd(1);
     gPad->SetLogy();
-    TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr);
-    TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
+    TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd);
+    TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd);
     CorrectFromTheWidth(measuredTH1Daftersubstraction);
     CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
     measuredTH1Daftersubstraction->SetStats(0);
@@ -1019,7 +1154,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
     }
     ratiomeasuredcontamination->Draw("P");
     cbackgroundsubtraction->cd(3);
-    TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr);
+    TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd);
     CorrectFromTheWidth(measuredTH1background);
     measuredTH1background->SetStats(0);
     measuredTH1background->SetTitle("");
@@ -1133,9 +1268,20 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
   //
   // calculate charm background
   //
+    Int_t ptpr = 0;
+    Int_t nDim = 1;
+    if(fBeamType==0)
+    {
+       ptpr=0;
+    }
+    if(fBeamType==1)
+    {
+       ptpr=1;
+       nDim=2;
+    }
 
   Double_t evtnorm=0;
-  if(fNMCbgEvents) evtnorm= double(fNEvents[0])/double(fNMCbgEvents);
+  if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]);
 
   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
   if(!mcContainer){
@@ -1150,29 +1296,58 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
 
   AliCFDataGrid *charmBackgroundGrid= 0x0;
   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
-  TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
 
+  TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
   Int_t* bins=new Int_t[1];
   bins[0]=charmbgaftertofpid->GetNbinsX();
-  AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",1,bins);
-  ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
-  TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
 
-  charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
-  TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
+  if(fBeamType==1)
+  {
+      charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
+      bins=new Int_t[2];
+      bins[0]=charmbgaftertofpid->GetNbinsX();
+      charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
+      bins[1]=charmbgaftertofpid->GetNbinsX();
 
-  AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins);
+  }
+
+  AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
+  ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
+  TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
+
+  if(fBeamType==0)charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
+  Double_t contents[2];
+  AliCFDataGrid *eventTemp = new AliCFDataGrid("eventTemp","eventTemp",nDim,bins);
+  if(fBeamType==1)
+  {
+      for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
+      {
+         for(Int_t kpt=0;kpt<bins[1];kpt++)
+         {
+             Double_t evtnormPbPb=0;
+             if(fNMCbgEvents[kCentr]) evtnormPbPb= double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
+             contents[0]=kCentr;
+             contents[1]=kpt;
+              eventTemp->Fill(contents,evtnormPbPb);
+         }
+      }
+   charmBackgroundGrid->Multiply(eventTemp,1);
+  }
+  TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
+
+  AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
   weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
-  TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
+  TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr);
   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
-  TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
+  TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr);
 
   // Efficiency (set efficiency to 1 for only folding) 
   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
   efficiencyD->CalculateEfficiency(0,0);
 
-  // Folding 
-  Int_t nDim = 1;
+  // Folding
+  if(fBeamType==0)nDim = 1;
+  if(fBeamType==1)nDim = 2;
   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
   folding.SetMaxNumberOfIterations(1);
   folding.Unfold();
@@ -1181,7 +1356,7 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
   THnSparse* result=(THnSparse*)result1->Clone();
   charmBackgroundGrid->SetGrid(result);
-  TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
+  TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr);
 
   //Charm background evaluation plots
 
@@ -1189,7 +1364,18 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
   cCharmBgEval->Divide(3,1);
 
   cCharmBgEval->cd(1);
-  charmbgaftertofpid->Scale(evtnorm);
+
+  if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
+  if(fBeamType==1)
+  {
+      Double_t evtnormPbPb=0;
+      for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
+      {
+         if(fNMCbgEvents[kCentr]) evtnormPbPb= evtnormPbPb+double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
+      }
+      charmbgaftertofpid->Scale(evtnormPbPb);
+  }
+
   CorrectFromTheWidth(charmbgaftertofpid);
   charmbgaftertofpid->SetMarkerStyle(25);
   charmbgaftertofpid->Draw("p");
@@ -1248,15 +1434,15 @@ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
   //
   
   Double_t evtnorm[1] = {0.0};
-  if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
+  if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
   printf("check event!!! %lf \n",evtnorm[0]);
   
   AliCFContainer *backgroundContainer = 0x0;
   
   if(fNonHFEsyst){
-    backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0]->Clone();
+    backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone();
     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
-      backgroundContainer->Add(fConvSourceContainer[iSource][0]);
+      backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent
     }  
   }
   else{    
@@ -1269,22 +1455,46 @@ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
   }
   
   Int_t stepbackground = 3;
+  if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4;
+
   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
-  //backgroundGrid->Scale(evtnorm);
-  //workaround for statistical errors: "Scale" method does not work for our errors, since Sumw2 is not defined for AliCFContainer
-  Int_t nBin[1];
-  
-  Int_t bins[1];
-  bins[0]=fConversionEff->GetNbinsX();
-  
-  for(Long_t iBin=1; iBin<= bins[0];iBin++){
-    nBin[0]=iBin;
-    backgroundGrid->SetElementError(nBin, backgroundGrid->GetElementError(nBin)*evtnorm[0]);
-    backgroundGrid->SetElement(nBin,backgroundGrid->GetElement(nBin)*evtnorm[0]);
+  Int_t *nBinpp=new Int_t[1];
+  Int_t* binspp=new Int_t[1];
+  binspp[0]=fConversionEff[0]->GetNbinsX();  // number of pt bins
+
+  Int_t *nBinPbPb=new Int_t[2];
+  Int_t* binsPbPb=new Int_t[2];
+  binsPbPb[1]=fConversionEff[0]->GetNbinsX();  // number of pt bins
+  binsPbPb[0]=12;
+
+  Int_t looppt=binspp[0];
+  if(fBeamType==1) looppt=binsPbPb[1];
+
+  for(Long_t iBin=1; iBin<= looppt;iBin++){
+      if(fBeamType==0)
+      {
+         nBinpp[0]=iBin;
+         backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
+         backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
+      }
+      if(fBeamType==1)
+      {
+          // loop over centrality
+         for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
+             nBinPbPb[0]=iiBin;
+             nBinPbPb[1]=iBin;
+              Double_t evtnormPbPb=0;
+              if(fNMCbgEvents[iiBin]) evtnormPbPb= double(fNEvents[iiBin])/double(fNMCbgEvents[iiBin]);
+             backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
+             backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
+         }
+      }
   }
   //end of workaround for statistical errors
-  
-  AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,&bins[0]);
+
+  AliCFDataGrid *weightedConversionContainer;
+  if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp);
+  else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
   weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
   backgroundGrid->Multiply(weightedConversionContainer,1.0);
   
@@ -1299,14 +1509,14 @@ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
   //
 
   Double_t evtnorm[1] = {0.0};
-  if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
+  if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
   printf("check event!!! %lf \n",evtnorm[0]);
   
   AliCFContainer *backgroundContainer = 0x0;
   if(fNonHFEsyst){
-    backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0]->Clone();
+    backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone();
     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
-      backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);
+      backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
     } 
   }
   else{    
@@ -1320,22 +1530,45 @@ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
   
   
   Int_t stepbackground = 3;
+  if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4;
+
   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
-  //backgroundGrid->Scale(evtnorm);
-  //workaround for statistical errors: "Scale" method does not work for our errors, since Sumw2 is not defined for AliCFContainer
-  Int_t nBin[1];
-  
-  Int_t bins[1];
-  bins[0]=fConversionEff->GetNbinsX();
-  
-  for(Long_t iBin=1; iBin<= bins[0];iBin++){
-    nBin[0]=iBin;
-    backgroundGrid->SetElementError(nBin, backgroundGrid->GetElementError(nBin)*evtnorm[0]);
-    backgroundGrid->SetElement(nBin,backgroundGrid->GetElement(nBin)*evtnorm[0]);
+  Int_t *nBinpp=new Int_t[1];
+  Int_t* binspp=new Int_t[1];
+  binspp[0]=fConversionEff[0]->GetNbinsX();  // number of pt bins
+
+  Int_t *nBinPbPb=new Int_t[2];
+  Int_t* binsPbPb=new Int_t[2];
+  binsPbPb[1]=fConversionEff[0]->GetNbinsX();  // number of pt bins
+  binsPbPb[0]=12;
+
+  Int_t looppt=binspp[0];
+  if(fBeamType==1) looppt=binsPbPb[1];
+
+
+  for(Long_t iBin=1; iBin<= looppt;iBin++){
+      if(fBeamType==0)
+      {
+         nBinpp[0]=iBin;
+         backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
+         backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
+      }
+      if(fBeamType==1)
+      {
+         for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
+             nBinPbPb[0]=iiBin;
+             nBinPbPb[1]=iBin;
+             Double_t evtnormPbPb=0;
+              if(fNMCbgEvents[iiBin]) evtnormPbPb= double(fNEvents[iiBin])/double(fNMCbgEvents[iiBin]);
+             backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
+             backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
+         }
+      }
   }
   //end of workaround for statistical errors
-  
-  AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,&bins[0]);
+  AliCFDataGrid *weightedNonHFEContainer;
+  if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
+  else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
   weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
   backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
 
@@ -1361,16 +1594,13 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons
       AliError("Data Container not available");
       return NULL;
     }
-
     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
   } 
-
   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
   result->SetName("ParametrizedEfficiencyBefore");
   THnSparse *h = result->GetGrid();
   Int_t nbdimensions = h->GetNdimensions();
   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
-
   AliCFContainer *dataContainer = GetContainer(kDataContainer);
   if(!dataContainer){
     AliError("Data Container not available");
@@ -1384,7 +1614,6 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons
   memset(coord, 0, sizeof(Int_t) * nbdimensions);
   Double_t* points = new Double_t[nbdimensions];
 
-
   ULong64_t nEntries = h->GetNbins();
   for (ULong64_t i = 0; i < nEntries; ++i) {
     
@@ -1421,7 +1650,7 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons
   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
 
   if(fDebugLevel > 0) {
-    
+
     TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
     cEfficiencyParametrized->Divide(2,1);
     cEfficiencyParametrized->cd(1);
@@ -1458,9 +1687,9 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons
     //ratioefficiency->Draw();
 
     if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
+
   }
 
-  
   return resultt;
 
 }
@@ -1673,15 +1902,18 @@ TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
   // Efficiency
   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
-  
-  // Consider parameterized IP cut efficiency
-  if(!fInclusiveSpectrum){
-    Int_t* bins=new Int_t[1];
-    bins[0]=35;
 
-    AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
-    beffContainer->SetGrid(GetBeautyIPEff());
-    efficiencyD->Multiply(beffContainer,1);  
+  if(!fBeauty2ndMethod)
+  {
+      // Consider parameterized IP cut efficiency
+      if(!fInclusiveSpectrum){
+         Int_t* bins=new Int_t[1];
+         bins[0]=35;
+
+         AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
+         beffContainer->SetGrid(GetBeautyIPEff());
+         efficiencyD->Multiply(beffContainer,1);
+      }
   }
   
 
@@ -1774,14 +2006,18 @@ AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpe
   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 
-  // Consider parameterized IP cut efficiency
-  if(!fInclusiveSpectrum){
-    Int_t* bins=new Int_t[1];
-    bins[0]=35;
-  
-    AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
-    beffContainer->SetGrid(GetBeautyIPEff());
-    efficiencyD->Multiply(beffContainer,1);
+
+  if(!fBeauty2ndMethod)
+  {
+      // Consider parameterized IP cut efficiency
+      if(!fInclusiveSpectrum){
+         Int_t* bins=new Int_t[1];
+         bins[0]=35;
+
+         AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
+         beffContainer->SetGrid(GetBeautyIPEff());
+         efficiencyD->Multiply(beffContainer,1);
+      }
   }
 
   // Data in the right format
@@ -1958,6 +2194,7 @@ TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i)
     TH1D* projection = (TH1D *) spectrum->Project(ptpr);
     CorrectFromTheWidth(projection);
     TGraphErrors *graphError = NormalizeTH1(projection,i);
+
     return graphError;
     
   }
@@ -1990,7 +2227,6 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
       n = input->GetBinContent(ibin);
       dN = input->GetBinError(ibin);
-
       // New point
       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
       errdN = 1./(2. * TMath::Pi() * p);
@@ -2005,7 +2241,7 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
     spectrumNormalized->SetMarkerStyle(22);
     spectrumNormalized->SetMarkerColor(kBlue);
     spectrumNormalized->SetLineColor(kBlue);
-    
+
     return spectrumNormalized;
     
   }
@@ -2038,7 +2274,6 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) cons
       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
       n = input->GetBinContent(ibin);
       dN = input->GetBinError(ibin);
-
       // New point
       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
       errdN = 1./(2. * TMath::Pi() * p);
@@ -2080,13 +2315,14 @@ AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type)
   return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
 }
 //____________________________________________________________
-AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge) {
+AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centrality) {
   //
   // Slice bin for a given source of electron
   // nDim is the number of dimension the corrections are done
   // dimensions are the definition of the dimensions
   // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
   // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
+  // centrality (-1 means we do not cut on centrality)
   //
   
   Double_t *varMin = new Double_t[container->GetNVar()],
@@ -2124,7 +2360,12 @@ AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, In
       fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
       AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
     }
-    
+    if(ivar == 5){
+       if((centrality>= 0) && (centrality<container->GetNBins(ivar))) {
+           varMin[ivar] = binLimits[centrality];
+           varMax[ivar] = binLimits[centrality];
+       }
+    }
 
     delete[] binLimits;
   }
@@ -2255,86 +2496,353 @@ THnSparse* AliHFEspectrum::GetCharmWeights(){
   // Measured D->e based weighting factors
   //
 
-  const Int_t nDim=1;
-  Int_t nBin[nDim] = {35};
-
+  const Int_t nDimpp=1;
+  Int_t nBinpp[nDimpp] = {35};
   Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
+  const Int_t nDimPbPb=2;
+  Int_t nBinPbPb[nDimPbPb] = {11,35};
+  Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+  Int_t loopcentr=1;
+  Int_t looppt=nBinpp[0];
+  if(fBeamType==0)
+  {
+      fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
+      fWeightCharm->SetBinEdges(0, ptbinning1);
+  }
+  if(fBeamType==1)
+  {
+      fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
+      fWeightCharm->SetBinEdges(1, ptbinning1);
+      fWeightCharm->SetBinEdges(0, kCentralityRange);
+      loopcentr=nBinPbPb[0];
+  }
+
 
-  fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin);
-  for(Int_t idim = 0; idim < nDim; idim++)
-     fWeightCharm->SetBinEdges(idim, ptbinning1);
      Double_t weight[35]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385};
 //{1.11865,1.11444,1.13395,1.15751,1.13412,1.14446,1.15372,1.09458,1.13446,1.11707,1.1352,1.10979,1.08887,1.11164,1.10772,1.10727,1.07027,1.07016,1.04826,1.03248,1.0093,0.973534,0.932574,0.869838,0.799316,0.734015,0.643001,0.584319,0.527147,0.495661,0.470002,0.441129,0.431156,0.417242,0.39246,0.379611,0.390845,0.368857,0.34959,0.387186,0.376364,0.384437,0.349585,0.383225};
   //points
   Double_t pt[1];
-  for(int i=0; i<nBin[0]; i++){
-    pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
-    fWeightCharm->Fill(pt,weight[i]);
+  Double_t contents[2];
+
+  for(int icentr=0; icentr<loopcentr; icentr++)
+  {
+      for(int i=0; i<looppt; i++){
+         pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
+         if(fBeamType==1)
+         {
+             contents[0]=icentr;
+             contents[1]=pt[0];
+         }
+         if(fBeamType==0)
+         {
+             contents[0]=pt[0];
+         }
+
+         fWeightCharm->Fill(contents,weight[i]);
+      }
+  }
+
+  Int_t nDimSparse = fWeightCharm->GetNdimensions();
+  Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+  Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+  for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
+      nBins *= binsvar[iVar];
+  }
+
+  Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+  // loop that sets 0 error in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bintmp = iBin ;
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+      bintmp /= binsvar[iVar] ;
+    }
+    fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
   }
-  Int_t* ibins[nDim];
-  for(Int_t ivar = 0; ivar < nDim; ivar++)
-    ibins[ivar] = new Int_t[nBin[ivar] + 1];
-  fWeightCharm->SetBinError(ibins[0],0);
 
   return fWeightCharm;
 }
 
 //____________________________________________________________________________
-THnSparse* AliHFEspectrum::GetBeautyIPEff(){
-  //
-  // Return beauty electron IP cut efficiency
-  //
+void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, Int_t *dimensions){
+
+   // TOF PID efficiencies
+   Int_t ptpr;
+   if(fBeamType==0) ptpr=0;
+   if(fBeamType==1) ptpr=1;
+
+   Int_t loopcentr=1;
+   const Int_t nCentralitybinning=11; //number of centrality bins
+   if(fBeamType==1)
+   {
+     loopcentr=nCentralitybinning;
+   }
 
-  const Int_t nDim=1;
-  Int_t nBin[nDim] = {35};
+   TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,20.);
+   TF1 *fipfit2 = new TF1("fipfit2","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,10.0);
+   TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",1.0,7.);
+
+   TCanvas * cefficiencyParam = new TCanvas("efficiencyParam","efficiencyParam",1200,600);
+   cefficiencyParam->Divide(2,1);
+   cefficiencyParam->cd(1);
+
+   AliCFContainer *mccontainermcD = 0x0;
+   TH1D* efficiencysigTOFPIDD[nCentralitybinning];
+   TH1D* efficiencyTOFPIDD[nCentralitybinning];
+   Int_t source = -1; //get parameterized TOF PID efficiencies
+
+   for(int icentr=0; icentr<loopcentr; icentr++) {
+      // signal sample
+      if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
+      if(fBeamType==1) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+      AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
+      efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+      // mb sample for double check
+      if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+      if(fBeamType==1)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+      AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
+      efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+      efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
+      efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
+      efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
+      efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
+
+      fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+      fittofpid->SetLineColor(3);
+      efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
+      efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
+      fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
+   }
 
-  Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
+   // draw (for PbPb, only 1st bin)
+   efficiencysigTOFPIDD[0]->SetTitle("");
+   efficiencysigTOFPIDD[0]->SetStats(0);
+   efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
+   efficiencysigTOFPIDD[0]->SetMarkerColor(2);
+   efficiencysigTOFPIDD[0]->SetLineColor(2);
+   efficiencysigTOFPIDD[0]->Draw();
 
-  THnSparseF *ipcut = new THnSparseF("beff", "b eff; pt[GeV/c]", nDim, nBin);
-  for(Int_t idim = 0; idim < nDim; idim++)
-     ipcut->SetBinEdges(idim, ptbinning1);
-  Double_t pt[1];
-  Double_t weight;
-  for(int i=0; i<nBin[0]; i++){
-    pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
-    weight=0.612*(0.385+0.111*log(pt[0])-0.00435*log(pt[0])*log(pt[0]))*tanh(5.42*pt[0]-1.22);  // for 3 sigma cut   
-    //weight=0.786*(0.493+0.0283*log(pt[0])-0.0190*log(pt[0])*log(pt[0]))*tanh(1.84*pt[0]-0.00368);  // for 2 sigma cut   
-    ipcut->Fill(pt,weight);
-  }
-  Int_t* ibins[nDim];
-  for(Int_t ivar = 0; ivar < nDim; ivar++)
-    ibins[ivar] = new Int_t[nBin[ivar] + 1];
-  ipcut->SetBinError(ibins[0],0);
+   fEfficiencyTOFPIDD[0]->Draw("same");
 
-  return ipcut;
+   efficiencyTOFPIDD[0]->SetTitle("");
+   efficiencyTOFPIDD[0]->SetStats(0);
+   efficiencyTOFPIDD[0]->SetMarkerStyle(24);
+   efficiencyTOFPIDD[0]->Draw("same");
+
+
+   cefficiencyParam->cd(2);
+
+   // IP cut efficiencies
+  
+   for(int icentr=0; icentr<loopcentr; icentr++)  {
+     source = 0; //charm
+     if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
+     if(fBeamType==1) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
+     efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+
+     source = 1; //beauty
+     if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
+     if(fBeamType==1) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
+     efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+
+     source = 0; //charm
+     if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
+     efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+
+     source = 1; //beauty
+     if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
+     efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+
+     source = 2; //conversion
+     if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
+     efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+
+     source = 3; //non HFE except for the conversion
+     if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
+     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
+     efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+     fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr);
+     fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr);
+     fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr);
+     fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr);
+     fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr);
+     fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr);
+   }
+
+   for(int icentr=0; icentr<loopcentr; icentr++)  {
+     fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+     fipfit->SetLineColor(2);
+     fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
+     fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
+     if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit");
+     else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
+
+     if(fIPParameterizedEff){
+       fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+       fipfit->SetLineColor(1);
+       fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
+       fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
+       fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
+
+       fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+       fipfit2->SetLineColor(3);
+       fConversionEff[icentr]->Fit(fipfit2,"R");
+       fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
+       fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfit2");
+
+       fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+       fipfit2->SetLineColor(4);
+       fNonHFEEff[icentr]->Fit(fipfit2,"R");
+       fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
+       fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfit2");
+     }
+   }
+
+   // draw (for PbPb, only 1st bin)
+
+   fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
+   fEfficiencyCharmSigD[0]->SetMarkerColor(1);
+   fEfficiencyCharmSigD[0]->SetLineColor(1);
+   fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
+   fEfficiencyBeautySigD[0]->SetMarkerColor(2);
+   fEfficiencyBeautySigD[0]->SetLineColor(2);
+   fCharmEff[0]->SetMarkerStyle(24);
+   fCharmEff[0]->SetMarkerColor(1);
+   fCharmEff[0]->SetLineColor(1);
+   fBeautyEff[0]->SetMarkerStyle(24);
+   fBeautyEff[0]->SetMarkerColor(2);
+   fBeautyEff[0]->SetLineColor(2);
+   fConversionEff[0]->SetMarkerStyle(24);
+   fConversionEff[0]->SetMarkerColor(3);
+   fConversionEff[0]->SetLineColor(3);
+   fNonHFEEff[0]->SetMarkerStyle(24);
+   fNonHFEEff[0]->SetMarkerColor(4);
+   fNonHFEEff[0]->SetLineColor(4);
+
+   fEfficiencyCharmSigD[0]->Draw();
+   fEfficiencyBeautySigD[0]->Draw("same");
+   fCharmEff[0]->Draw("same");
+   fBeautyEff[0]->Draw("same");
+   fConversionEff[0]->Draw("same");
+   fNonHFEEff[0]->Draw("same");
+
+   fEfficiencyIPBeautyD[0]->Draw("same");
+   if(fIPParameterizedEff){
+     fEfficiencyIPCharmD[0]->Draw("same");
+     fEfficiencyIPConversionD[0]->Draw("same");
+     fEfficiencyIPNonhfeD[0]->Draw("same");
+   }
 }
 
 //____________________________________________________________________________
-THnSparse* AliHFEspectrum::GetCharmEff(){
+THnSparse* AliHFEspectrum::GetBeautyIPEff(){
   //
-  // Return charm electron IP cut efficiency
+  // Return beauty electron IP cut efficiency
   //
 
-  const Int_t nDim=1;
-  Int_t nBin[nDim] = {35};
+  const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
+  const Int_t nCentralitybinning=11;//number of centrality bins
+  Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
+  Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+  Int_t ptpr = 0;
+  Int_t nDim=1;  //dimensions of the efficiency weighting grid
+  if(fBeamType==1)
+  {
+    ptpr=1;
+    nDim=2; //dimensions of the efficiency weighting grid
+  }
+  Int_t nBin[1] = {nPtbinning1};
+  Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
 
-  Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
 
-  THnSparseF *ipcut = new THnSparseF("ceff", "c eff; pt[GeV/c]", nDim, nBin);
+  THnSparseF *ipcut;
+  if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
+  else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
   for(Int_t idim = 0; idim < nDim; idim++)
-     ipcut->SetBinEdges(idim, ptbinning1);
+  {
+    if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
+    if(nDim==2)
+      {
+        ipcut->SetBinEdges(0, kCentralityRange);
+        ipcut->SetBinEdges(1, kPtRange);
+      }
+  }
   Double_t pt[1];
-  Double_t weight;
-  for(int i=0; i<nBin[0]; i++){
-    pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
-    weight=0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
-    //weight=0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
-    ipcut->Fill(pt,weight);
-  }
-  Int_t* ibins[nDim];
-  for(Int_t ivar = 0; ivar < nDim; ivar++)
-    ibins[ivar] = new Int_t[nBin[ivar] + 1];
-  ipcut->SetBinError(ibins[0],0);
+  Double_t weight[nCentralitybinning];
+  Double_t contents[2];
+
+  for(Int_t a=0;a<11;a++)
+  {
+      weight[a] = 1.0;
+  }
+
+
+  Int_t looppt=nBin[0];
+  Int_t loopcentr=1;
+  if(fBeamType==1)
+  {
+      loopcentr=nBinPbPb[0];
+  }
+
+
+  for(int icentr=0; icentr<loopcentr; icentr++)
+  {
+      for(int i=0; i<looppt; i++)
+      {
+         pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+          if(fEfficiencyIPBeautyD[icentr])
+            weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
+          else{
+            printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
+            weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1); 
+          }
+
+          if(fBeamType==1){
+              contents[0]=icentr;
+              contents[1]=pt[0];
+          }
+          if(fBeamType==0){
+              contents[0]=pt[0];
+          }
+          ipcut->Fill(contents,weight[icentr]);
+      }
+  } 
+
+
+  Int_t nDimSparse = ipcut->GetNdimensions();
+  Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+  Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+  for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
+      nBins *= binsvar[iVar];
+  }
+
+  Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+  // loop that sets 0 error in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bintmp = iBin ;
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+      bintmp /= binsvar[iVar] ;
+    }
+    ipcut->SetBinError(binfill,0.); // put 0 everywhere
+  }
 
   return ipcut;
 }
@@ -2344,54 +2852,255 @@ THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
   //
   // Return PID x IP cut efficiency
   //
+    const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
+    const Int_t nCentralitybinning=11;//number of centrality bins
+    Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
+    Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+    Int_t ptpr = 0;
+    Int_t nDim=1;  //dimensions of the efficiency weighting grid
+    if(fBeamType==1)
+    {
+       ptpr=1;
+       nDim=2; //dimensions of the efficiency weighting grid
+    }
+    Int_t nBin[1] = {nPtbinning1};
+    Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
+
+    THnSparseF *pideff;
+    if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
+    else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
+    for(Int_t idim = 0; idim < nDim; idim++)
+    {
+
+       if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
+       if(nDim==2)
+       {
+           pideff->SetBinEdges(0, kCentralityRange);
+           pideff->SetBinEdges(1, kPtRange);
+       }
+    }
 
-  const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
-  const Int_t nDim=1;//dimensions of the efficiency weighting grid
-  Int_t nBin[nDim] = {nPtbinning1};
+  Double_t pt[1];
+  Double_t weight[nCentralitybinning];
+  Double_t contents[2];
+
+  for(Int_t a=0;a<11;a++)
+  {
+      weight[a] = 1.0;
+  }
+
+  Int_t looppt=nBin[0];
+  Int_t loopcentr=1;
+  if(fBeamType==1)
+  {
+      loopcentr=nBinPbPb[0];
+  }
+
+  for(int icentr=0; icentr<loopcentr; icentr++)
+  {
+      Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
+      for(int i=0; i<looppt; i++)
+      {
+         pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+
+          Double_t tofpideff = 0.;
+          if(fEfficiencyTOFPIDD[icentr])
+            tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
+          else{
+            printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
+          }  
+
+          //tof pid eff x tpc pid eff x ip cut eff
+          if(fIPParameterizedEff){
+            if(source==0) {
+              if(fEfficiencyIPCharmD[icentr])
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
+              else{
+                printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
+              }
+            } 
+           else if(source==2) {
+              if(fEfficiencyIPConversionD[icentr])
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]); 
+              else{
+                printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
+              }
+            }
+           else if(source==3) {
+              if(fEfficiencyIPNonhfeD[icentr])
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]); 
+              else{
+                printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
+              }  
+            }
+          }
+          else{
+            if(source==0) weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1); //charm
+           else if(source==2) weight[icentr] = tofpideff*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
+           else if(source==3) weight[icentr] = tofpideff*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
+          }
+
+         if(fBeamType==1){
+             contents[0]=icentr;
+             contents[1]=pt[0];
+         }
+         if(fBeamType==0){
+             contents[0]=pt[0];
+         }
+
+         pideff->Fill(contents,weight[icentr]);
+      }
+  }
+
+  Int_t nDimSparse = pideff->GetNdimensions();
+  Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+  Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+  for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
+      nBins *= binsvar[iVar];
+  }
+
+  Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+  // loop that sets 0 error in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bintmp = iBin ;
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+      bintmp /= binsvar[iVar] ;
+    }
+    pideff->SetBinError(binfill,0.); // put 0 everywhere
+  }
 
-  Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
 
-  THnSparseF *pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
-  for(Int_t idim = 0; idim < nDim; idim++)
-     pideff->SetBinEdges(idim, kPtRange);
 
-  Double_t pt[1];
-  Double_t weight = 1.0;
-  Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
-  for(int i=0; i<nBin[0]; i++){
-    pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
-    Double_t tofeff = 0.9885*(0.8551+0.0070*log(pt[0])-0.0014*log(pt[0])*log(pt[0]))*tanh(0.8884*pt[0]+0.6061); // parameterized TOF PID eff
-
-    if(source==0) weight = tofeff*trdtpcPidEfficiency*0.299*(0.307+0.176*log(pt[0])+0.0176*log(pt[0])*log(pt[0]))*tanh(96.3*pt[0]-19.7); // for 3 sigma cut
-    //if(source==0) weight = tofeff*trdtpcPidEfficiency*0.477*(0.3757+0.184*log(pt[0])-0.0013*log(pt[0])*log(pt[0]))*tanh(436.013*pt[0]-96.2137); // for 2 sigma cut
-    //if(source==2) weight = tofeff*trdtpcPidEfficiency*0.5575*(0.3594-0.3051*log(pt[0])+0.1597*log(pt[0])*log(pt[0]))*tanh(8.2436*pt[0]-0.8125);
-    //if(source==3) weight = tofeff*trdtpcPidEfficiency*0.0937*(0.4449-0.3769*log(pt[0])+0.5031*log(pt[0])*log(pt[0]))*tanh(6.4063*pt[0]+3.1425); // multiply ip cut eff for Dalitz/dielectrons
-    if(source==2) weight = tofeff*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // multiply ip cut eff for conversion electrons
-    if(source==3) weight = tofeff*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // multiply ip cut eff for Dalitz/dielectrons
-    //printf("tofeff= %lf trdtpcPidEfficiency= %lf conversion= %lf nonhfe= %lf\n",tofeff,trdtpcPidEfficiency,fConversionEff->GetBinContent(i+1),fNonHFEEff->GetBinContent(i+1));
-
-    pideff->Fill(pt,weight);
-  }
-  Int_t* ibins[nDim];
-  for(Int_t ivar = 0; ivar < nDim; ivar++)
-    ibins[ivar] = new Int_t[nBin[ivar] + 1];
-  pideff->SetBinError(ibins[0],0);
 
   return pideff;
 }
+
+//__________________________________________________________________________
+AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
+ //
+ // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
+    //
+    Int_t ptpr = 0;
+    Int_t nDim = 1;
+    if(fBeamType==0)
+    {
+       ptpr=0;
+    }
+    if(fBeamType==1)
+    {
+       ptpr=1;
+       nDim=2;
+    }
+
+    const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
+    const Int_t nCentralitybinning=11;//number of centrality bins
+    Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
+    Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
+    Int_t nBin[1] = {nPtbinning1};
+    Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
+    //fBSpectrum2ndMethod->GetNbinsX()
+
+    AliCFDataGrid *rawBeautyContainer;
+    if(fBeamType==0)  rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
+    else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb);
+    //  printf("number of bins= %d\n",bins[0]);
+
+
+    
+    
+    THnSparseF *brawspectra;
+    if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
+    else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb);
+    if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange);
+    if(fBeamType==1)
+      {
+       //      brawspectra->SetBinEdges(0, centralityBins);
+       brawspectra->SetBinEdges(0, kCentralityRange);
+       brawspectra->SetBinEdges(1, kPtRange);
+      }
+    
+    Double_t pt[1];
+    Double_t yields= 0.;
+    Double_t valuesb[2];
+    
+    //Int_t looppt=nBin[0];
+    Int_t loopcentr=1;
+    if(fBeamType==1)
+      {
+       loopcentr=nBinPbPb[0];
+      }
+    
+    for(int icentr=0; icentr<loopcentr; icentr++)
+      {
+       
+       for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
+         pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+         
+         yields = fBSpectrum2ndMethod->GetBinContent(i+1);
+         
+         if(fBeamType==1)
+           {
+             valuesb[0]=icentr;
+             valuesb[1]=pt[0];
+           }
+         if(fBeamType==0) valuesb[0]=pt[0];
+         brawspectra->Fill(valuesb,yields);
+       }
+      }
+    
+    
+    
+    Int_t nDimSparse = brawspectra->GetNdimensions();
+    Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
+    Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+    
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
+      nBins *= binsvar[iVar];
+    }
+    
+    Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+    // loop that sets 0 error in each bin
+    for (Long_t iBin=0; iBin<nBins; iBin++) {
+      Long_t bintmp = iBin ;
+      for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+       bintmp /= binsvar[iVar] ;
+      }
+      brawspectra->SetBinError(binfill,0.); // put 0 everywhere
+    }
+    
+    
+    rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
+    TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr);
+    
+    new TCanvas;
+    fBSpectrum2ndMethod->SetMarkerStyle(24);
+    fBSpectrum2ndMethod->Draw("p");
+    hRawBeautySpectra->SetMarkerStyle(25);
+    hRawBeautySpectra->Draw("samep");
+    
+    return rawBeautyContainer;
+}
+
 //__________________________________________________________________________
-void AliHFEspectrum::CalculateNonHFEsyst(){
+void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
+  //
+  // Calculate non HFE sys
   //
-  //Calculate systematic errors for non-HFE and conversion background,
-  //based on correlated errors from pi0 input and uncorrelated errors 
-  //from m_t scaling
   //
+
   if(!fNonHFEsyst)
     return;
 
   Double_t evtnorm[1] = {0.0};
-  if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents);
+  if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
   
   AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
   AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
@@ -2402,22 +3111,23 @@ void AliHFEspectrum::CalculateNonHFEsyst(){
 
   Int_t stepbackground = 3;
   Int_t* bins=new Int_t[1];
-  bins[0]=fConversionEff->GetNbinsX();
+  bins[0]=fConversionEff[centrality]->GetNbinsX();
+   
   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
 
-  for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
-   
+  for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){   
     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
-      convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground); 
+      convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),*fConvSourceContainer[iSource][iLevel][centrality],stepbackground);
       weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
       convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
-
-      nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
+      
+      nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),*fNonHFESourceContainer[iSource][iLevel][centrality],stepbackground);
       weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
       nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
     }
-
+    
     bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
       bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
@@ -2427,11 +3137,12 @@ void AliHFEspectrum::CalculateNonHFEsyst(){
     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){//add other sources to get overall background from all meson decays
       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
     }
-        
+    
     bgLevelGrid[iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
     bgLevelGrid[iLevel]->Add(bgNonHFEGrid[iLevel]);
   }
   
+  
   //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated)
   AliCFDataGrid *bgErrorGrid[2];
   bgErrorGrid[0] = (AliCFDataGrid*)bgLevelGrid[1]->Clone();
@@ -2588,5 +3299,7 @@ void AliHFEspectrum::CalculateNonHFEsyst(){
   gOverallSystErrLow->Write(); 
 
   output->Close(); 
-  delete output;  
+  delete output;
+  
 }
+
index 8a1b70c..15f2c96 100644 (file)
@@ -52,7 +52,8 @@ class AliHFEspectrum : public TNamed{
     enum{
       kElecBgSources = 6,
       kBgLevels = 3,
-      kBgPtBins = 44
+      kBgPtBins = 44,
+      kCentrality = 12
        };
 
    enum Chargetype_t{
@@ -87,10 +88,12 @@ class AliHFEspectrum : public TNamed{
     void SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type);
     void SetEfficiencyFunction(TF1 *efficiencyFunction) { fEfficiencyFunction = efficiencyFunction; };
     void SetPbPbAnalysis(Bool_t isPbPb = kFALSE) { fBeamType=(Char_t) isPbPb; };
+
+    void SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, Int_t *dimensions);
     
     void SetNumberOfEvents(Int_t nEvents,Int_t i = 0) { fNEvents[i] = nEvents; };
     void SetNumberOfMCEvents(Int_t nEvents) { fNMCEvents = nEvents; };
-    void SetNumberOfMC2Events(Int_t nEvents) { fNMCbgEvents = nEvents; }; 
+    void SetNumberOfMC2Events(Int_t nEvents,Int_t i = 0) { fNMCbgEvents[i] = nEvents; };
     void SetMCEffStep(Int_t step) { fStepMC = step; };
     void SetMCTruthStep(Int_t step) { fStepTrue = step; };
     void SetStepToCorrect(Int_t step) { fStepData = step; };
@@ -106,8 +109,10 @@ class AliHFEspectrum : public TNamed{
     void SetLowHighBoundaryCentralityBinAtTheEnd(Int_t low, Int_t high, Int_t i) { fLowBoundaryCentralityBinAtTheEnd[i] = low; fHighBoundaryCentralityBinAtTheEnd[i] = high;};
 
     void SetBeautyAnalysis() { fInclusiveSpectrum = kFALSE; };
+    void CallInputFileForBeauty2ndMethod();
+    void SetInputFileForBeauty2ndMethod(const char *filenameb = "BSpectrum2ndmethod.root"){fkBeauty2ndMethodfilename = filenameb; };
+    void SetBeautyAnalysis2ndMethod(Bool_t beauty2ndmethod) { fBeauty2ndMethod = beauty2ndmethod; }
     void SetHadronEffbyIPcut(THnSparseF* hsHadronEffbyIPcut) { fHadronEffbyIPcut = hsHadronEffbyIPcut;};
-    void SetNonHFEBackground2ndMethod() { fNonHFEbgMethod2 = kTRUE; };
     void SetNonHFEsyst(Bool_t syst){ fNonHFEsyst = syst; };
 
     void SetStepGuessedUnfolding(Int_t stepGuessedUnfolding) { fStepGuessedUnfolding = stepGuessedUnfolding; };
@@ -117,24 +122,26 @@ class AliHFEspectrum : public TNamed{
   
     void SetDebugLevel(Int_t debugLevel, Bool_t writeToFile = kFALSE) { fDebugLevel = debugLevel; fWriteToFile = writeToFile; };
 
+
+    AliCFDataGrid* GetRawBspectra2ndMethod();
     AliCFDataGrid* GetCharmBackground();
     AliCFDataGrid* GetConversionBackground();
     AliCFDataGrid* GetNonHFEBackground();
     THnSparse* GetCharmWeights();
     THnSparse* GetBeautyIPEff();
-    THnSparse* GetCharmEff();
     THnSparse* GetPIDxIPEff(Int_t source);
-    void CalculateNonHFEsyst();
+    void CalculateNonHFEsyst(Int_t centrality = 0);
 
     void EnableIPanaHadronBgSubtract() { fIPanaHadronBgSubtract = kTRUE; };
     void EnableIPanaCharmBgSubtract() { fIPanaCharmBgSubtract = kTRUE; };
     void EnableIPanaConversionBgSubtract() { fIPanaConversionBgSubtract = kTRUE; };
     void EnableIPanaNonHFEBgSubtract() { fIPanaNonHFEBgSubtract = kTRUE; };
+    void EnableIPParameterizedEff() { fIPParameterizedEff = kTRUE; };
 
   protected:
        
     AliCFContainer *GetContainer(AliHFEspectrum::CFContainer_t contt);
-    AliCFContainer *GetSlicedContainer(AliCFContainer *cont, Int_t ndim, Int_t *dimensions,Int_t source=-1,Chargetype_t charge=kAllCharge);
+    AliCFContainer *GetSlicedContainer(AliCFContainer *cont, Int_t ndim, Int_t *dimensions,Int_t source=-1,Chargetype_t charge=kAllCharge,Int_t centrality=-1);
     THnSparseF *GetSlicedCorrelation(THnSparseF *correlationmatrix,Int_t nDim, Int_t *dimensions) const;
     TObject* GetSpectrum(const AliCFContainer * const c, Int_t step);
     TObject* GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0);
@@ -154,11 +161,16 @@ class AliHFEspectrum : public TNamed{
     THnSparseF *fCorrelation;     // Correlation Matrices
     AliCFDataGrid *fBackground;   // Background Grid
     TF1 *fEfficiencyFunction;     // Efficiency Function
+    TF1 *fEfficiencyTOFPIDD[kCentrality];       // TOF PID efficiency parameterized
+    TF1 *fEfficiencyIPCharmD[kCentrality];      // IP efficiency parameterized for charm
+    TF1 *fEfficiencyIPBeautyD[kCentrality];     // IP efficiency parameterized for beauty 
+    TF1 *fEfficiencyIPConversionD[kCentrality]; // IP efficiency parameterized for conversion
+    TF1 *fEfficiencyIPNonhfeD[kCentrality];     // IP efficiency parameterized for nonhfe 
 
     THnSparseF *fWeightCharm;     // Weight for charm bg
 
-    AliCFContainer *fConvSourceContainer[kElecBgSources][kBgLevels]; //container for conversion electrons, divided into different photon sources
-    AliCFContainer *fNonHFESourceContainer[kElecBgSources][kBgLevels];     //container for non-HF electrons, divided into different sources
+    AliCFContainer *fConvSourceContainer[kElecBgSources][kBgLevels][kCentrality]; //container for conversion electrons, divided into different photon sources
+    AliCFContainer *fNonHFESourceContainer[kElecBgSources][kBgLevels][kCentrality]; //container for non-HF electrons, divided into different sources
 
     Bool_t fInclusiveSpectrum;     // Inclusive Spectrum
     Bool_t fDumpToFile;           // Write Result in a file
@@ -171,13 +183,14 @@ class AliHFEspectrum : public TNamed{
     Bool_t fIPanaCharmBgSubtract;      // Charm background subtraction 
     Bool_t fIPanaConversionBgSubtract; // Conversion background subtraction
     Bool_t fIPanaNonHFEBgSubtract;     // nonHFE except for conversion background subtraction
-    Bool_t fNonHFEbgMethod2;           // switch for 2nd method to subtract non HFE background
+    Bool_t fIPParameterizedEff;        // switch to use parameterized efficiency for ip analysis
     Bool_t fNonHFEsyst;            // choose NonHFE background level (upper, lower, central)
+    Bool_t fBeauty2ndMethod;      // 2nd method to get beauty spectrum
 
     Int_t fNbDimensions;          // Number of dimensions for the correction
     Int_t fNEvents[20];           // Number of Events
     Int_t fNMCEvents;             // Number of MC Events
-    Int_t fNMCbgEvents;           // Number of BG MC Events
+    Int_t fNMCbgEvents[20];       // Number of BG MC Events
     Int_t fStepMC;                // MC step (for unfolding)
     Int_t fStepTrue;              // MC step of the final spectrum
     Int_t fStepData;              // Data Step (various applications)
@@ -195,8 +208,14 @@ class AliHFEspectrum : public TNamed{
     Int_t fHighBoundaryCentralityBinAtTheEnd[20];  // Boundary of the bins
 
     THnSparseF *fHadronEffbyIPcut;// container for hadron efficiency by IP cut
-    TH1D *fConversionEff;         // conversion IP cut eff
-    TH1D *fNonHFEEff;             // nonhfe IP cut eff
+    TH1D *fEfficiencyCharmSigD[kCentrality]; // charm IP cut eff from signal enhanced MC
+    TH1D *fEfficiencyBeautySigD[kCentrality]; // beauty IP cut eff from signal enhanced MC
+    TH1D *fConversionEff[kCentrality];     // conversion IP cut eff
+    TH1D *fNonHFEEff[kCentrality];         // nonhfe IP cut eff
+    TH1D *fCharmEff[kCentrality];              // charm IP cut eff
+    TH1D *fBeautyEff[kCentrality];             // beauty IP cut eff
+    TH1D *fBSpectrum2ndMethod;             // beauty spectrum for 2nd method
+    const char *fkBeauty2ndMethodfilename;      // name of file, which contains beauty spectrum for 2ndmethod
     Char_t fBeamType;             // beamtype; default -1; pp =0; PbPb=1
 
 
index c2a3d09..ec20ab3 100644 (file)
 #include <TMath.h>
 #include <TParticle.h>
 #include "TH1.h"
+#include "TH1D.h"
 #include "TH2.h"
+#include "TGraph.h"
+#include "TGraphErrors.h"
 #include "THnSparse.h"
 #include "TAxis.h"
+#include "TMath.h"
+#include "TString.h"
 
 #include "AliAODMCParticle.h"
 #include "AliAODpidUtil.h"
@@ -242,3 +247,159 @@ Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
   };
   return pid;
 }
+
+//__________________________________________
+TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
+{
+    // Creates a TH1D from TGraph g. The binwidth of the first bin has to
+    // specified. The others bins are calculated automatically. Supports also Graphs
+    // with non constant x steps. The axis of the Graph can be exchanged if
+    // exchange=kTRUE (modifies the Graph).
+
+    TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize);
+    if( result == 0) return result;
+
+    //------------------------------------------
+    // setup the final hist
+    Int_t nBinX = g->GetN();
+    Double_t* err = g->GetEY();           // error y is still ok even if exchanged
+    for(Int_t i = 0; i < nBinX; i ++){
+       result->SetBinError(i + 1, err[i]);
+    }
+    if(exchange){
+        AliHFEtools::ExchangeXYGraph(g);        // undo  what has been done in GraphToHist
+       AliHFEtools::ExchangeXYGraphErrors(g);  // now exchange including errors
+    }
+
+    return result;
+}
+
+//__________________________________________
+Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g)
+{
+    // exchanges x-values and y-values.
+    if(g==0) return kFALSE;
+    Int_t nbin=g->GetN();
+    Double_t x,y;
+    for(Int_t i = 0; i < nbin; i ++)
+    {
+        g->GetPoint(i,x,y);
+        g->SetPoint(i,y,x);
+    }
+
+    return kTRUE;
+}
+
+//__________________________________________
+Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g)
+{
+    // exchanges x-values and y-values and
+    // corresponding errors.
+    if(g==0) return kFALSE;
+    Int_t nbin=g->GetN();
+    Double_t x,y;
+    Double_t ex,ey;
+    for(Int_t i = 0; i < nbin; i ++)
+    {
+        g->GetPoint(i,x,y);
+        ex = g->GetErrorX(i);
+        ey = g->GetErrorY(i);
+        g->SetPoint(i,y,x);
+        g->SetPointError(i,ey,ex);
+    }
+
+    return kTRUE;
+
+}
+
+//__________________________________________
+TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
+{
+    // Creates a TH1D from TGraph g. The binwidth of the first bin has to
+    // specified. The others bins are calculated automatically. Supports also Graphs
+    // with non constant x steps. The axis of the Graph can be exchanged if
+    // exchange=kTRUE (modifies the Graph).
+
+
+    TH1D* result = 0;
+    if(g == 0)              return result;
+    if(firstBinWidth == -1) return result;
+    TString myname="";
+    myname = g->GetName();
+    myname += "_graph";
+    if(exchange) AliHFEtools::ExchangeXYGraph(g);
+
+    Int_t nBinX = g->GetN();
+    Double_t* x = g->GetX();
+    Double_t* y = g->GetY();
+
+    if(nBinX < 1) return result;
+
+    //------------------------------------------
+    // create the Matrix for the equation system
+    // and init the values
+
+    Int_t nDim = nBinX - 1;
+    TMatrixD a(nDim,nDim);
+    TMatrixD b(nDim,1);
+
+    Double_t* aA = a.GetMatrixArray();
+    Double_t* aB = b.GetMatrixArray();
+    memset(aA,0,nDim * nDim * sizeof(Double_t));
+    memset(aB,0,nDim * sizeof(Double_t));
+    //------------------------------------------
+
+    //------------------------------------------
+    // setup equation system
+    // width for 1st bin is given therefore
+    // we shift bin parameter (column) by one to the left
+    // to reduce the matrix size
+
+    Double_t* xAxis = new Double_t [nBinX + 1];
+    Double_t* binW  = new Double_t [nBinX ];
+    binW[0] = firstBinWidth;
+
+    aB[0] = x[1] - x[0]  - 0.5 * binW[0];
+    aA[0] = 0.5;
+
+    for(Int_t col = 1; col < nDim ; col ++)
+    {
+       Int_t row = col;
+       aB[col] = x[col + 1] - x[ col ];
+       aA[row * nDim + col - 1 ] = 0.5;
+       aA[row * nDim + col     ] = 0.5;
+    }
+    //------------------------------------------
+
+    //------------------------------------------
+    // solve the equations
+    a.Invert();
+    TMatrixD c = a * b;
+    //------------------------------------------
+
+    //------------------------------------------
+    // calculate the bin boundaries
+    xAxis[0] = x[0] - 0.5 * binW[0];
+    memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t));
+    for(Int_t col = 0; col < nBinX ; col ++) {
+       xAxis[col + 1] = x[col] + 0.5 * binW[col];
+    }
+    //------------------------------------------
+
+    //------------------------------------------
+    // setup the final hist
+    result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis);
+    for(Int_t i = 0; i < nBinX; i ++){
+       result->SetBinContent(i + 1, y[i]);
+    }
+    result->SetMarkerColor(markercolor);
+    result->SetMarkerStyle(markerstyle);
+    result->SetMarkerSize(markersize);
+    //------------------------------------------
+
+    delete [] xAxis;
+    delete [] binW;
+
+
+    return result;
+}
index b88469f..f548bcb 100644 (file)
@@ -26,6 +26,10 @@ class TParticle;
 class AliAODMCParticle;
 class AliPIDResponse;
 class AliVParticle;
+class TGraph;
+class TGraphErrors;
+class TH1D;
+class TString;
 
 class AliHFEtools : public TObject{
   public:
@@ -42,6 +46,10 @@ class AliHFEtools : public TObject{
     static AliPIDResponse *GetDefaultPID(Bool_t isMC = kTRUE, Bool_t isESD = kTRUE);
     static void DestroyDefaultPID();
     static void SetLogLevel(Int_t loglevel) { fgLogLevel = loglevel ;}
+    static TH1D* GraphErrorsToHist(TGraphErrors* g = 0, Double_t firstBinWidth = -1, Bool_t exchange=kFALSE, Int_t markerstyle=8, Int_t markercolor=2, Float_t markersize=0.7);
+    static TH1D* GraphToHist(TGraph* g = 0, Double_t firstBinWidth = -1, Bool_t exchange=kFALSE, Int_t markerstyle=8, Int_t markercolor=2, Float_t markersize=0.7);
+    static Bool_t ExchangeXYGraph(TGraph* g = 0);
+    static Bool_t ExchangeXYGraphErrors(TGraphErrors* g = 0);
 
   private:
       AliHFEtools(const AliHFEtools &);
index cbfda41..e402dfa 100644 (file)
@@ -224,15 +224,15 @@ void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
   Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
   nbins[kElectronLike] = 100;
   nbins[kNClustersLike] = 200;
-  nbins[kCentralityBin] = 12;
+  nbins[kCentralityBinLike] = 12;
   Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   binMin[kElectronLike] = 0.;      
   binMin[kNClustersLike] = 0.;
-  binMin[kCentralityBin] = -1.;
+  binMin[kCentralityBinLike] = -1.;
   binMax[kElectronLike] = 1.;
   binMax[kNClustersLike] = 200.;
-  binMax[kCentralityBin] = 11.;
+  binMax[kCentralityBinLike] = 11.;
 
   fHistos->CreateTHnSparse("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
   fHistos->BinLogAxis("fLikeTRD", kP);
@@ -247,13 +247,15 @@ void AliHFEtrdPIDqa::CreateQAHistogram(){
   Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
   nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
   nbins[kNClusters] = 200;
+  nbins[kCentralityBinQA] = 12;
   Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   binMin[kNonZeroTrackletCharge] = 0.;      
-  binMin[kNClusters] = 0.;      
+  binMin[kNClusters] = 0.;
+  binMin[kCentralityBinQA] = -1.;
   Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
   binMax[kNClusters] = 200.;
-
+  binMax[kCentralityBinQA] = 11.;
   fHistos->CreateTHnSparse("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
   fHistos->BinLogAxis("fQAtrack", kP);
 }
@@ -267,15 +269,18 @@ void AliHFEtrdPIDqa::CreatedEdxHistogram(){
   Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
   nbins[kdEdx] = 100;
   nbins[kNclusters] = 261;
-  nbins[kNonZeroSlices] = 9; 
+  nbins[kNonZeroSlices] = 9;
+  nbins[kCentralityBindEdx] = 12;
   Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   binMin[kdEdx] = 0.;     
   binMin[kNclusters] = 0;
   binMin[kNonZeroSlices] = 0.;
+  binMin[kCentralityBindEdx] = -1.;
   Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   binMax[kdEdx] = 10000.;
   binMax[kNclusters] = 260.;
   binMax[kNonZeroSlices] = 8.;
+  binMax[kCentralityBindEdx] = 11.;
 
   fHistos->CreateTHnSparse("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
   fHistos->BinLogAxis("fQAdEdx", kP);
@@ -292,14 +297,17 @@ void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
   nbins[kTPCdEdx] = 600;
   nbins[kTRDdEdxMethod1] = 1000;
   nbins[kTRDdEdxMethod2] = 1000;
+  nbins[kCentralityBinTruncMean] = 12;
   Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   binMin[kTPCdEdx] = 0.;      
   binMin[kTRDdEdxMethod1] = 0.;      
-  binMin[kTRDdEdxMethod2] = 0.;      
+  binMin[kTRDdEdxMethod2] = 0.;
+  binMin[kCentralityBinTruncMean] = -1.;
   Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
   binMax[kTPCdEdx] = 600;
   binMax[kTRDdEdxMethod1] = 20000.;
   binMax[kTRDdEdxMethod2] = 20000.;
+  binMax[kCentralityBinTruncMean] = 11.;
 
   fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
   fHistos->BinLogAxis("fTRDtruncMean", kP);
@@ -375,12 +383,13 @@ void AliHFEtrdPIDqa::FillTRDLikelihoods(const AliESDtrack * const track, Int_t s
   // p
   // ntracklets
   // Electron Likelihood
+  // Centrality
   quantities[kSpecies] = species;
   quantities[kP] = outerPars ? outerPars->P() : track->P();
   quantities[kNTracklets] = track->GetTRDntrackletsPID();
   quantities[kElectronLike] = likeEle;
   quantities[kNClustersLike] =  track->GetTRDncls();
-  quantities[kCentralityBin] = fCentralityBin;
+  quantities[kCentralityBinLike] = fCentralityBin;
   fHistos->Fill("fLikeTRD", quantities);
 
 }
@@ -455,14 +464,17 @@ void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t speci
     quantitiesdEdx[kdEdx] = fTotalChargeInSlice0 ? track->GetTRDslice(iplane, 0) : dEdxSum; // hack by mfasel: In the new reconstruction, the total charge is stored in the first slice, in the old reconstruction it has to be calculated from the slice charges.     
     if(dEdxSum) ntrackletsNonZero++;
     // Fill dEdx histogram
+    quantitiesdEdx[kCentralityBindEdx] = fCentralityBin;
     if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
   }
   quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
+  quantitiesQA[kCentralityBinQA] = fCentralityBin;
   fHistos->Fill("fQAtrack", quantitiesQA);
 
   quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
   quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, 0.6);
   quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, 0.6);
+  quantitiesTruncMean[kCentralityBinTruncMean] = fCentralityBin;
   fHistos->Fill("fTRDtruncMean", quantitiesTruncMean);
 }
 
@@ -609,8 +621,8 @@ void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets, Int_t nCentrality, Bool
   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, isGreaterEqual ? 7 : binTracklets);
 
   if(isPbPb){
-      Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBin)->FindBin(nCentrality);
-      hLikeTRD->GetAxis(kCentralityBin)->SetRange(binCentrality, isGreaterEqual ? 11 : binCentrality);
+      Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(nCentrality);
+      hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, isGreaterEqual ? 11 : binCentrality);
       /*
        new TCanvas;
        TH2 *test = hLikeTRD->Projection(kCentralityBin, kP);
@@ -637,7 +649,7 @@ void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets, Int_t nCentrality, Bool
   // Undo ranges
   hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
-  hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kCentralityBin)->GetNbins());
+  hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kCentralityBinLike)->GetNbins());
 
   // Prepare List for output
   Char_t *listname=Form("%dTracklets", nTracklets);
@@ -842,8 +854,8 @@ Double_t AliHFEtrdPIDqa::CalculateIntegratedPionEfficiency(UInt_t nTracklets, Do
   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
 
   if(icentrality!=1){
-      Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBin)->FindBin(icentrality);
-      hLikeTRD->GetAxis(kCentralityBin)->SetRange(binCentrality, binCentrality);
+      Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(icentrality);
+      hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, binCentrality);
   }
 
   Int_t pbinMin = hLikeTRD->GetAxis(kP)->FindBin(pmax),
index 656deb2..7c5ef34 100644 (file)
@@ -119,25 +119,28 @@ class AliHFEtrdPIDqa : public TNamed{
     enum QuantitiesLike_t{
       kElectronLike = 3,
       kNClustersLike = 4,
-      kCentralityBin = 5,
+      kCentralityBinLike = 5,
       kQuantitiesLike = 6
     };
     enum QuantitiesQAtrack_t{
       kNonZeroTrackletCharge = 3,
       kNClusters = 4,
-      kQuantitiesQA = 5
+      kCentralityBinQA = 5,
+      kQuantitiesQA = 6
     };
     enum QuantitiesdEdx_t{
       kNclusters = 3,
       kNonZeroSlices = 4,
       kdEdx = 5,
-      kQuantitiesdEdx = 6
+      kCentralityBindEdx = 6,
+      kQuantitiesdEdx = 7
     };
     enum QuantitiesTruncMean_t{
       kTPCdEdx = 3,
       kTRDdEdxMethod1 = 4,
       kTRDdEdxMethod2 = 5,
-      kQuantitiesTruncMean = 6
+      kCentralityBinTruncMean = 6,
+      kQuantitiesTruncMean = 7
     };
 
     void ProcessTrackESD(const AliESDtrack * const track, Int_t species);