Update of hfe code
authorrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Nov 2012 16:47:58 +0000 (16:47 +0000)
committerrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Nov 2012 16:47:58 +0000 (16:47 +0000)
27 files changed:
PWGHF/hfe/AliAnalysisTaskHFE.cxx
PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx
PWGHF/hfe/AliAnalysisTaskHFEFlow.h
PWGHF/hfe/AliHFECorrectSpectrumBase.cxx [new file with mode: 0644]
PWGHF/hfe/AliHFECorrectSpectrumBase.h [new file with mode: 0644]
PWGHF/hfe/AliHFEInclusiveSpectrum.cxx [new file with mode: 0644]
PWGHF/hfe/AliHFEInclusiveSpectrum.h [new file with mode: 0644]
PWGHF/hfe/AliHFEInclusiveSpectrumQA.cxx [new file with mode: 0644]
PWGHF/hfe/AliHFEInclusiveSpectrumQA.h [new file with mode: 0644]
PWGHF/hfe/AliHFENonPhotonicElectron.cxx
PWGHF/hfe/AliHFENonPhotonicElectron.h
PWGHF/hfe/AliHFEdebugTreeTask.cxx
PWGHF/hfe/AliHFEdebugTreeTaskAOD.cxx
PWGHF/hfe/AliHFEdebugTreeTaskAOD.h
PWGHF/hfe/AliHFEextraCuts.cxx
PWGHF/hfe/AliHFEextraCuts.h
PWGHF/hfe/AliHFEmcQA.cxx
PWGHF/hfe/AliHFEpidObject.cxx
PWGHF/hfe/AliHFEpidObject.h
PWGHF/hfe/AliHFEpidTOF.cxx
PWGHF/hfe/AliHFEpidTOF.h
PWGHF/hfe/AliHFEpidTPC.cxx
PWGHF/hfe/AliHFEpidTPC.h
PWGHF/hfe/AliHFEsignalCuts.cxx
PWGHF/hfe/AliHFEsmearDCA.cxx [new file with mode: 0644]
PWGHF/hfe/AliHFEsmearDCA.h [new file with mode: 0644]
PWGHF/hfe/AliHFEspectrum.cxx

index 7e7fb82..41a261d 100644 (file)
@@ -626,7 +626,11 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
     if(!mcH->InitOk()) return;
     if(!mcH->TreeK()) return;
     if(!mcH->TreeTR()) return;
-    if(GetPlugin(kNonPhotonicElectron) && fBackgroundSubtraction) fBackgroundSubtraction->SetMCEvent(fMCEvent);
+    // Background subtraction-------------------------------------------------------------------
+    if(GetPlugin(kNonPhotonicElectron) && fBackgroundSubtraction){
+      fBackgroundSubtraction->SetMCEvent(fMCEvent);
+    }
+    //------------------------------------------------------------------------------------------
   }
 
   if(IsAODanalysis() && HasMCData()){
@@ -963,6 +967,13 @@ void AliAnalysisTaskHFE::ProcessESD(){
 
   fCFM->SetRecEventInfo(fESD);
 
+  // Get Number of contributors to the primary vertex for multiplicity-dependent correction
+  Int_t ncontribVtx = 0;
+  const AliESDVertex *priVtx = fESD->GetPrimaryVertexTracks();
+  if(priVtx){
+    ncontribVtx = priVtx->GetNContributors();
+  }
+
   // minjung for IP QA(temporary ~ 2weeks)
   if(!fExtraCuts){
     fExtraCuts = new AliHFEextraCuts("hfetmpCuts","HFE tmp Cuts");
@@ -1226,6 +1237,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
     hfetrack.SetRecTrack(track);
     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
     hfetrack.SetCentrality(fCentralityF);
+    hfetrack.SetMulitplicity(ncontribVtx);
     if(IsPbPb()) hfetrack.SetPbPb();
     else hfetrack.SetPP();
     fPID->SetVarManager(fVarManager);
@@ -1235,8 +1247,10 @@ void AliAnalysisTaskHFE::ProcessESD(){
     // Background subtraction------------------------------------------------------------------------------------------
     if (GetPlugin(kNonPhotonicElectron)) {
       Int_t indexmother = -1;
-      Int_t mcsource = 1;
-      if(HasMCData()) mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
+      Int_t mcsource = -1;
+      if(HasMCData()){
+       mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
+      }
       fBackgroundSubtraction->LookAtNonHFE(itrack, track, fInputEvent, 1, fCentralityF, -1, mcsource, indexmother);
     }
     //-----------------------------------------------------------------------------------------------------------------
@@ -1534,6 +1548,18 @@ void AliAnalysisTaskHFE::ProcessAOD(){
 
   fCFM->SetRecEventInfo(fAOD);
 
+  if(!fExtraCuts){
+    fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
+  }
+  fExtraCuts->SetRecEventInfo(fAOD);
+
+  // Get Number of contributors to the primary vertex for multiplicity-dependent correction
+  Int_t ncontribVtx = 0;
+  AliAODVertex *priVtx = fAOD->GetPrimaryVertex();
+  if(priVtx){
+    ncontribVtx = priVtx->GetNContributors();
+  }
+
   // Look for kink mother
   Int_t numberofvertices = fAOD->GetNumberOfVertices();
   Double_t listofmotherkink[numberofvertices];
@@ -1562,6 +1588,7 @@ void AliAnalysisTaskHFE::ProcessAOD(){
   AliAODTrack *track = NULL;
   AliAODMCParticle *mctrack = NULL;
   Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
+  Double_t dataDca[6]; // [source, pT, dca, centrality]
   Int_t nElectronCandidates = 0;
   Int_t pid;
   Bool_t signal;
@@ -1644,6 +1671,53 @@ void AliAnalysisTaskHFE::ProcessAOD(){
       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
     }
 
+    if(HasMCData()){
+      Double_t hfeimpactR4all=0., hfeimpactnsigmaR4all=0.;
+      Int_t sourceDca =-1;
+      if(mctrack && (TMath::Abs(mctrack->GetPdgCode()) == 211)){
+        if(track->Pt()>4.){
+          fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
+          dataDca[0]=0; //pion
+          dataDca[1]=track->Pt();
+          dataDca[2]=hfeimpactR4all;
+          dataDca[3]=fCentralityF;
+          dataDca[4] = -1; // not store V0 for the moment
+          dataDca[5] = double(track->Charge());
+          fQACollection->Fill("Dca", dataDca);
+        }
+      }
+      else if(mctrack && (TMath::Abs(mctrack->GetPdgCode()) == 11)){ // to increas statistics for Martin
+        if(signal){
+          fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4all, hfeimpactnsigmaR4all);
+          if(fSignalCuts->IsCharmElectron(track)){
+            sourceDca=1;
+          }
+          else if(fSignalCuts->IsBeautyElectron(track)){
+            sourceDca=2;
+          }
+          else if(fSignalCuts->IsGammaElectron(track)){
+            sourceDca=3;
+          }
+          else if(fSignalCuts->IsNonHFElectron(track)){
+            sourceDca=4;
+          }
+          else if(fSignalCuts->IsJpsiElectron(track)){
+            sourceDca=5;
+          }
+          else {
+            sourceDca=6;
+          }
+          dataDca[0]=sourceDca;
+          dataDca[1]=track->Pt();
+          dataDca[2]=hfeimpactR4all;
+          dataDca[3]=fCentralityF;
+          dataDca[4] = -1; // not store V0 for the moment
+          dataDca[5] = double(track->Charge());
+          if(signal) fQACollection->Fill("Dca", dataDca);
+        }
+      }
+    }
+
     //printf("Will process to PID\n");
 
     // track accepted, do PID
@@ -1652,6 +1726,7 @@ void AliAnalysisTaskHFE::ProcessAOD(){
     hfetrack.SetRecTrack(track);
     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
     hfetrack.SetCentrality(fCentralityF);
+    hfetrack.SetMulitplicity(ncontribVtx); // for correction
     if(IsPbPb()) hfetrack.SetPbPb();
     else hfetrack.SetPP();
     fPID->SetVarManager(fVarManager);
@@ -1661,8 +1736,8 @@ void AliAnalysisTaskHFE::ProcessAOD(){
     // Background subtraction----------------------------------------------------------------------------------------------
     if (GetPlugin(kNonPhotonicElectron)) {
       Int_t indexmother = -1;
-      Int_t mcsource = 1;
-      if(HasMCData()) mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
+      Int_t mcsource = -1;
+      if(HasMCData())  mcsource = fBackgroundSubtraction->FindMother(mctrack->GetLabel(),indexmother);
       fBackgroundSubtraction->LookAtNonHFE(itrack, track, fInputEvent, 1, fCentralityF, -1,mcsource, indexmother);
     }
     //---------------------------------------------------------------------------------------------------------------------
@@ -1730,6 +1805,33 @@ void AliAnalysisTaskHFE::ProcessAOD(){
         fQACollection->Fill("PIDperformance", dataE);
       }
     }
+
+    if (GetPlugin(kDEstep)) {
+      if (!HasMCData()){
+        Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
+        fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
+        dataDca[0]=-1; //for data, don't know the origin
+        dataDca[1]=track->Pt();
+        dataDca[2]=hfeimpactR;
+        dataDca[3]=fCentralityF;
+        dataDca[4] = -1; // not store V0 for the moment
+        dataDca[5] = double(track->Charge());
+        fQACollection->Fill("Dca", dataDca);
+      }
+
+      // Fill Containers for impact parameter analysis
+      if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
+      if(signal) {
+        // Apply weight for background contamination after ip cut
+        if(fBackGroundFactorApply) {
+              fWeightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
+              if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
+              else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
+              // weightBackGround as special weight
+              fVarManager->FillContainer(fContainer, "hadronicBackground", 2, kFALSE, fWeightBackGround);
+        }
+      }
+    }
   }
   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
 }
index 76400be..32c58f3 100644 (file)
@@ -32,6 +32,9 @@
 #include "TParticle.h"\r
 #include "TF1.h"\r
 \r
+#include <TDirectory.h>\r
+#include <TTreeStream.h>\r
+\r
 #include "AliVEventHandler.h"\r
 #include "AliAnalysisTaskSE.h"\r
 #include "AliAnalysisManager.h"\r
@@ -109,6 +112,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :
   fMaxNumberOfIterations(100),\r
   fPrecisionPhi(0.001),\r
   fUseMCReactionPlane(kFALSE),\r
+  fSP(kFALSE),\r
   fMCPID(kFALSE),\r
   fNoPID(kFALSE),\r
   fChi2OverNDFCut(3.0),\r
@@ -176,7 +180,8 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() :
   fOppSignDeltaPhiMaps(0x0),\r
   fSameSignDeltaPhiMaps(0x0),\r
   fOppSignAngle(0x0),\r
-  fSameSignAngle(0x0)\r
+  fSameSignAngle(0x0),\r
+  fDebugStreamer(0)\r
 {\r
   // Constructor\r
 \r
@@ -219,6 +224,7 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :
   fMaxNumberOfIterations(100),\r
   fPrecisionPhi(0.001),\r
   fUseMCReactionPlane(kFALSE),\r
+  fSP(kFALSE),\r
   fMCPID(kFALSE),\r
   fNoPID(kFALSE),\r
   fChi2OverNDFCut(3.0),\r
@@ -286,7 +292,8 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) :
   fOppSignDeltaPhiMaps(0x0),\r
   fSameSignDeltaPhiMaps(0x0),\r
   fOppSignAngle(0x0),\r
-  fSameSignAngle(0x0)\r
+  fSameSignAngle(0x0),\r
+  fDebugStreamer(0)\r
 {\r
   //\r
   // named ctor\r
@@ -351,6 +358,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref
   fMaxNumberOfIterations(ref.fMaxNumberOfIterations),\r
   fPrecisionPhi(ref.fPrecisionPhi),\r
   fUseMCReactionPlane(ref.fUseMCReactionPlane),\r
+  fSP(ref.fSP),\r
   fMCPID(ref.fMCPID),\r
   fNoPID(ref.fNoPID),\r
   fChi2OverNDFCut(ref.fChi2OverNDFCut),\r
@@ -418,7 +426,8 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref
   fOppSignDeltaPhiMaps(NULL),\r
   fSameSignDeltaPhiMaps(NULL),\r
   fOppSignAngle(NULL),\r
-  fSameSignAngle(NULL)\r
+  fSameSignAngle(NULL),\r
+  fDebugStreamer(0)\r
 {\r
   //\r
   // Copy Constructor\r
@@ -468,6 +477,7 @@ void AliAnalysisTaskHFEFlow::Copy(TObject &o) const {
   target.fMaxNumberOfIterations = fMaxNumberOfIterations;\r
   target.fPrecisionPhi = fPrecisionPhi;\r
   target.fUseMCReactionPlane = fUseMCReactionPlane;\r
+  target.fSP = fSP;\r
   target.fMCPID = fMCPID;\r
   target.fNoPID = fNoPID;\r
   target.fChi2OverNDFCut = fChi2OverNDFCut;\r
@@ -497,7 +507,7 @@ void AliAnalysisTaskHFEFlow::Copy(TObject &o) const {
     target.fContamination[k] = fContamination[k];\r
     target.fv2contamination[k] = fv2contamination[k];\r
   }\r
\r
+  target.fDebugStreamer=fDebugStreamer;\r
 }\r
 //____________________________________________________________\r
 AliAnalysisTaskHFEFlow::~AliAnalysisTaskHFEFlow(){\r
@@ -518,7 +528,7 @@ AliAnalysisTaskHFEFlow::~AliAnalysisTaskHFEFlow(){
   if(fBackgroundSubtraction) delete fBackgroundSubtraction;\r
   //if(fPIDBackgroundqa) delete fPIDBackgroundqa;\r
   //if(fHFEVZEROEventPlane) delete fHFEVZEROEventPlane;\r
\r
+  if ( fDebugStreamer ) delete fDebugStreamer;\r
 \r
 }\r
 //________________________________________________________________________\r
@@ -718,6 +728,12 @@ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects()
   Double_t maxCos = 1.0;\r
   Double_t binLimCos[nBinsCos+1];\r
   for(Int_t i=0; i<=nBinsCos; i++) binLimCos[i]=(Double_t)minCos + (maxCos-minCos)/nBinsCos*(Double_t)i ;\r
+\r
+  // Int_t nBinsCosSP = 50;\r
+  // Double_t minCosSP = -100.0;\r
+  // Double_t maxCosSP = 100.0;\r
+  // Double_t binLimCosSP[nBinsCosSP+1];\r
+  // for(Int_t i=0; i<=nBinsCosSP; i++) binLimCosSP[i]=(Double_t)minCosSP + (maxCosSP-minCosSP)/nBinsCosSP*(Double_t)i ;\r
  \r
   Int_t nBinsC = 11;\r
   Double_t minC = 0.0;\r
@@ -1426,10 +1442,11 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   Float_t eventPlaneV0A = 0.0;\r
   Float_t eventPlaneV0C = 0.0;\r
   Float_t eventPlaneV0 = 0.0;\r
-  TVector2 *standardQ = 0x0;\r
+  TVector2 *qTPC = 0x0;\r
   TVector2 *qsub1a = 0x0;\r
   TVector2 *qsub2a = 0x0;\r
-\r
+  TVector2 qV0A,qV0C,qV0,*qAna;\r
+  \r
   // V0\r
 \r
   if(fHFEVZEROEventPlane && (!fAODAnalysis)){\r
@@ -1460,25 +1477,27 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   }\r
   else {\r
     \r
-    eventPlaneV0 = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0", fInputEvent,2));\r
-    AliDebug(2,Form("eventPlaneV0 %f",eventPlaneV0));\r
+    Double_t qVx, qVy;  //TR: info\r
+    eventPlaneV0 = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,10,2,qVx,qVy));\r
     if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
-    AliDebug(2,Form("eventPlaneV0 %f",eventPlaneV0));\r
-    eventPlaneV0A = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0A", fInputEvent,2));\r
+    qV0.Set(qVx,qVy);\r
+    eventPlaneV0A = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,8,2,qVx,qVy));\r
     if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi();\r
-    eventPlaneV0C = TVector2::Phi_0_2pi(vEPa->GetEventplane("V0C", fInputEvent,2));\r
+    qV0A.Set(qVx,qVy);\r
+    eventPlaneV0C = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,9,2,qVx,qVy));\r
     if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi();\r
+    qV0C.Set(qVx,qVy);\r
   \r
   }\r
 \r
   // TPC\r
 \r
-  standardQ = vEPa->GetQVector(); \r
+  qTPC = vEPa->GetQVector(); \r
   Double_t qx = -1.0;\r
   Double_t qy = -1.0;\r
-  if(standardQ) {\r
-    qx = standardQ->X();\r
-    qy = standardQ->Y();\r
+  if(qTPC) {\r
+    qx = qTPC->X();\r
+    qy = qTPC->Y();\r
   }  \r
   TVector2 qVectorfortrack;\r
   qVectorfortrack.Set(qx,qy);\r
@@ -1486,10 +1505,22 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
 \r
   // Choose the one used for v2\r
 \r
-  if(fVZEROEventPlane) eventPlanea = eventPlaneV0;\r
-  if(fVZEROEventPlaneA) eventPlanea = eventPlaneV0A;\r
-  if(fVZEROEventPlaneC) eventPlanea = eventPlaneV0C;\r
-  if(!fVZEROEventPlane) eventPlanea = eventPlaneTPC;\r
+  if(fVZEROEventPlane){ //TR: info\r
+    eventPlanea = eventPlaneV0;\r
+    qAna = &qV0;\r
+  }\r
+  if(fVZEROEventPlaneA){\r
+    eventPlanea = eventPlaneV0A;\r
+    qAna = &qV0A;\r
+  }\r
+  if(fVZEROEventPlaneC){\r
+    eventPlanea = eventPlaneV0C;\r
+    qAna = &qV0C;\r
+  }\r
+  if(!fVZEROEventPlane){\r
+    eventPlanea = eventPlaneTPC;\r
+    qAna = &qV0C;\r
+  }\r
 \r
   valuecossinephiep[0] = TMath::Cos(2*eventPlanea);\r
   valuecossinephiep[1] = TMath::Sin(2*eventPlanea);\r
@@ -1515,19 +1546,65 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
     diffsub1sub2asin = TMath::Sin(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));\r
   }\r
 \r
+  // if ( !fDebugStreamer ) {\r
+  //   //debug stream\r
+  //   TDirectory *backup = gDirectory;\r
+  //   fDebugStreamer = new TTreeSRedirector("TaskHFEFlowdebug.root");\r
+  //   if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer\r
+  // }     \r
+\r
+  // {\r
+\r
+  //   double v0nrom = TMath::Sqrt(qV0.X()*qV0.X()+qV0.Y()*qV0.Y());\r
+  //   double v0Anrom = TMath::Sqrt(qV0A.X()*qV0A.X()+qV0A.Y()*qV0A.Y());\r
+  //   double v0Cnrom = TMath::Sqrt(qV0C.X()*qV0C.X()+qV0C.Y()*qV0C.Y());\r
+  //   double sub1nrom = TMath::Sqrt(qsub1a->X()*qsub1a->X()+qsub1a->Y()*qsub1a->Y());\r
+  //   double sub2nrom = TMath::Sqrt(qsub2a->X()*qsub2a->X()+qsub2a->Y()*qsub2a->Y());\r
+\r
+  //   (* fDebugStreamer) << "UserExec" <<\r
+  //     "binct="<<binct<<\r
+  //     "qV0="<<v0nrom<<\r
+  //     "qV0A="<<v0Anrom<<\r
+  //     "qV0C="<<v0Cnrom<<\r
+  //     "qsub1a="<<sub1nrom<<\r
+  //     "qsub2a="<<sub2nrom<<\r
+  //     "\n";\r
+  // }\r
+\r
   // three sub events in case of VZEROA and VZEROC\r
-  diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C));\r
-  diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC));\r
-  diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC));\r
+  if(!fSP){\r
+    diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C));  //TR: \r
+    diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC));  //TR: \r
+    diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC));  //TR: \r
+  }\r
+  else{\r
+    if(fVZEROEventPlaneA){\r
+      diffsubasubb = qV0A.X()*qV0C.X()+qV0A.Y()*qV0C.Y();\r
+      diffsubasubc = qV0A.X()*qTPC->X()+qV0A.Y()*qTPC->Y();\r
+      diffsubbsubc = qV0C.X()*qTPC->X()+qV0C.Y()*qTPC->Y();\r
+    }\r
+    else if(fVZEROEventPlaneC){\r
+      diffsubasubb = qV0C.X()*qV0A.X()+qV0C.Y()*qV0A.Y();\r
+      diffsubasubc = qV0C.X()*qTPC->X()+qV0C.Y()*qTPC->Y();\r
+      diffsubbsubc = qV0A.X()*qTPC->X()+qV0A.Y()*qTPC->Y();\r
+    }\r
+  }\r
 \r
   diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneV0C));\r
   diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneTPC));\r
   diffsubbsubcsin = TMath::Sin(2.*(eventPlaneV0C - eventPlaneTPC));\r
   // three sub events in case of VZERO all\r
   if(fVZEROEventPlane && (!fVZEROEventPlaneA) && (!fVZEROEventPlaneC)) {\r
-    diffsubasubb = TMath::Cos(2.*(eventPlaneV0 - eventPlanesub1a));\r
-    diffsubasubc = TMath::Cos(2.*(eventPlaneV0 - eventPlanesub2a));\r
-    diffsubbsubc = TMath::Cos(2.*(eventPlanesub1a - eventPlanesub2a));\r
+    if(!fSP){\r
+      diffsubasubb = TMath::Cos(2.*(eventPlaneV0 - eventPlanesub1a));     //TR: \r
+      diffsubasubc = TMath::Cos(2.*(eventPlaneV0 - eventPlanesub2a));     //TR: \r
+      diffsubbsubc = TMath::Cos(2.*(eventPlanesub1a - eventPlanesub2a));  //TR: \r
+    }\r
+    else{\r
+      diffsubasubb = qV0.X()*qsub1a->X()+qV0.Y()*qsub1a->Y();         \r
+      diffsubasubc = qV0.X()*qsub2a->X()+qV0.Y()*qsub2a->Y();   \r
+      diffsubbsubc = qsub1a->X()*qsub2a->X()+qsub1a->Y()*qsub2a->Y();\r
+    }\r
     \r
     diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0 - eventPlanesub1a));\r
     diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0 - eventPlanesub2a));\r
@@ -1540,14 +1617,14 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
   ////////////////////////////////////////////////////////\r
   \r
   //if(!fVZEROEventPlane) {\r
-  // if(!standardQ) {\r
+  // if(!qTPC) {\r
   //  eventplanedefined = kFALSE;\r
-      //PostData(1, fListHist);\r
-      //return;\r
+  //PostData(1, fListHist);\r
+  //return;\r
   //}\r
   //}\r
 \r
-  if((!standardQ) || (!qsub1a) || (!qsub2a)) {\r
+  if((!qTPC) || (!qsub1a) || (!qsub2a)) {\r
     AliDebug(2,"No event plane");\r
     return;\r
   }\r
@@ -1623,16 +1700,16 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
     }\r
     else {\r
       valuensparsefbis[0] = diffsubasubb;\r
-      valuensparsefbis[1] = diffsubbsubc;\r
-      valuensparsefbis[2] = diffsubasubc;\r
-      fCosResabc->Fill(&valuensparsefbis[0]);\r
+      valuensparsefbis[1] = diffsubasubc;     //TR: WTF!!!\r
+      valuensparsefbis[2] = diffsubbsubc;\r
+      fCosResabc->Fill(&valuensparsefbis[0]); //TR: info\r
       valuensparsefbissin[0] = diffsubasubbsin;\r
       valuensparsefbissin[1] = diffsubbsubcsin;\r
       valuensparsefbissin[2] = diffsubasubcsin;\r
       if(fMonitorEventPlane) fSinResabc->Fill(&valuensparsefbissin[0]);\r
       if(fMonitorEventPlane) {\r
        fProfileCosResab->Fill(valuensparsefbis[3],valuensparsefbis[0]);\r
-       fProfileCosResac->Fill(valuensparsefbis[3],valuensparsefbis[1]);\r
+       fProfileCosResac->Fill(valuensparsefbis[3],valuensparsefbis[1]);  //TR: WTF!!!\r
        fProfileCosResbc->Fill(valuensparsefbis[3],valuensparsefbis[2]);\r
       }\r
     }\r
@@ -1761,16 +1838,16 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
        if(!(aodtrack->TestFilterBit(fFilter))) continue;  // Only process AOD tracks where the HFE is set\r
       }\r
     }\r
-    \r
+\r
     if(fApplyCut) {\r
 \r
       valuetrackingcuts[0] = track->Pt(); \r
       valuetrackingcuts[1] = 0;\r
-\r
\r
       // RecKine: ITSTPC cuts  \r
       if(!fHFECuts->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)) continue;\r
       if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);\r
-      \r
+\r
       // Reject kink mother\r
       if(fRejectKinkMother) {\r
        if(fAODAnalysis) {\r
@@ -1819,7 +1896,6 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
       if(fMonitorTrackCuts) fTrackingCuts->Fill(&valuetrackingcuts[0]);    \r
       \r
     }\r
-\r
     AliDebug(2,"Survived");\r
 \r
     /////////////////////////////////////////////////////////\r
@@ -1830,8 +1906,8 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
     //if(eventplanedefined && (!fVZEROEventPlane)) {\r
     if(!fVZEROEventPlane) {\r
       // Subtract the tracks from the event plane\r
-      Double_t qX = standardQ->X() - vEPa->GetQContributionX(track);  //Modify the components: subtract the track you want to look at with your analysis\r
-      Double_t qY = standardQ->Y() - vEPa->GetQContributionY(track);  //Modify the components: subtract the track you want to look at with your analysis\r
+      Double_t qX = qTPC->X() - vEPa->GetQContributionX(track);  //Modify the components: subtract the track you want to look at with your analysis\r
+      Double_t qY = qTPC->Y() - vEPa->GetQContributionY(track);  //Modify the components: subtract the track you want to look at with your analysis\r
       TVector2 newQVectorfortrack;\r
       newQVectorfortrack.Set(qX,qY);\r
       eventplanesubtracted = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; \r
@@ -1851,7 +1927,7 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
        else fillEventPlane = kFALSE;\r
       }\r
     }\r
-    \r
+\r
     ///////////////\r
     // MC\r
     //////////////\r
@@ -2028,7 +2104,7 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
     //\r
     valuensparseh[0] = TMath::Cos(2*TVector2::Phi_mpi_pi(phitrack-eventplanesubtracted));\r
     if(fillEventPlane) {\r
-      fCosPhiMaps->Fill(&valuensparseh[0]);\r
+      fCosPhiMaps->Fill(&valuensparseh[0]); //TR: fCosPhiQSum+=valuensparseh[0]*TMath:Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y()); fCosPhiQN++;\r
       if((valuefractioncont[1] >=0) && (valuefractioncont[1] < 11)){\r
        if(fContamination[((Int_t)valuefractioncont[1])]){\r
          Double_t weight = fContamination[((Int_t)valuefractioncont[1])]->Eval(track->P());\r
@@ -2046,7 +2122,9 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/)
        }     \r
       }\r
       if(fMonitorEventPlane) {\r
-       fProfileCosPhiMaps->Fill(valuensparsehprofile[1],valuensparsehprofile[2],valuensparseh[0]);\r
+       if(fSP)\r
+         valuensparseh[0] *= TMath::Sqrt(qAna->X()*qAna->X()+qAna->Y()*qAna->Y());\r
+       fProfileCosPhiMaps->Fill(valuensparsehprofile[1],valuensparsehprofile[2],valuensparseh[0]);  //TR: info\r
       }\r
     }\r
     \r
index 74940ae..f25cab4 100644 (file)
@@ -48,9 +48,9 @@ class TArrayI;
 class AliAODMCHeader;\r
 class TClonesArray;\r
 class AliHFENonPhotonicElectron;\r
+class TTreeSRedirector;\r
 \r
 class AliAnalysisTaskHFEFlow: public AliAnalysisTaskSE {\r
-  \r
 public:\r
 \r
   typedef enum{\r
@@ -122,6 +122,7 @@ public:
   void SetMaxNumberOfIterations(Int_t maxNumberOfIterations) { fMaxNumberOfIterations = maxNumberOfIterations; };\r
   void SetPrecisionPhi(Double_t precisionPhi) { fPrecisionPhi = precisionPhi;};\r
   void SetUseMCReactionPlane(Bool_t useMCReactionPlane) { fUseMCReactionPlane = useMCReactionPlane;};\r
+  void SetUseSP(Bool_t useSP) { fSP = useSP;}\r
   void SetMCPID(Bool_t mcPID) { fMCPID = mcPID;};\r
   void SetNoPID(Bool_t noPID) { fNoPID = noPID;};\r
 \r
@@ -180,6 +181,7 @@ private:
   Int_t     fMaxNumberOfIterations; // Max number of iteration for adding v2\r
   Double_t  fPrecisionPhi;  // precision phi for adding v2\r
   Bool_t    fUseMCReactionPlane; // use MC reaction plane\r
+  Bool_t    fSP;        // calculate using scalar product method (instead of event plane method)\r
 \r
   Bool_t    fMCPID; // MC PID for electrons\r
   Bool_t    fNoPID; // No PID for checks\r
@@ -300,7 +302,9 @@ private:
   THnSparseF *fOppSignAngle;         // ! Opening Angles\r
   THnSparseF *fSameSignAngle;        // ! Opening Angles\r
 \r
-  Int_t FindMother(Int_t tr, AliMCEvent *mcEvent, Int_t &indexmother);\r
+  TTreeSRedirector  *fDebugStreamer;               //!Debug streamer\r
+\r
+ Int_t FindMother(Int_t tr, AliMCEvent *mcEvent, Int_t &indexmother);\r
   Int_t CheckPdg(Int_t tr, AliMCEvent* mcEvent);\r
   Int_t IsMotherGamma(Int_t tr, AliMCEvent* mcEvent);\r
   Int_t IsMotherPi0(Int_t tr, AliMCEvent* mcEvent);\r
diff --git a/PWGHF/hfe/AliHFECorrectSpectrumBase.cxx b/PWGHF/hfe/AliHFECorrectSpectrumBase.cxx
new file mode 100644 (file)
index 0000000..fb33505
--- /dev/null
@@ -0,0 +1,514 @@
+
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// The following containers have to be set:
+//  - Correction framework container for real data
+//  - Correction framework container for MC (Efficiency Map)
+//  - Correction framework container for background coming from data
+//  - Correction framework container for background coming from MC
+//
+//  Author: 
+//            Raphaelle Bailhache <R.Bailhache@gsi.de>
+//            Markus Fasel <M.Fasel@gsi.de>
+//
+
+#include <TArrayD.h>
+#include <TH1.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TPad.h>
+#include <TH2D.h>
+#include <TF1.h>
+
+#include "AliPID.h"
+#include "AliCFContainer.h"
+#include "AliCFDataGrid.h"
+#include "AliCFEffGrid.h"
+#include "AliCFGridSparse.h"
+#include "AliCFUnfolding.h"
+#include "AliLog.h"
+
+#include "AliHFECorrectSpectrumBase.h"
+#include "AliHFEcuts.h"
+#include "AliHFEcontainer.h"
+#include "AliHFEtools.h"
+
+ClassImp(AliHFECorrectSpectrumBase)
+
+//____________________________________________________________
+AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const char *name):
+  TNamed(name, ""),
+  fCFContainers(new TObjArray(kDataContainerV0+1)),
+  fCorrelation(NULL),
+  fEfficiencyFunction(NULL),
+  fEtaSelected(kFALSE),
+  fSetSmoothing(kFALSE),
+  fNbDimensions(1),
+  fNEvents(0),
+  fStepMC(-1),
+  fStepTrue(-1),
+  fStepData(-1),
+  fStepBeforeCutsV0(-1),
+  fStepAfterCutsV0(-1),
+  fStepGuessedUnfolding(-1),
+  fNumberOfIterations(10),
+  fChargeChoosen(kAllCharge),
+  fTestCentralityLow(-1),
+  fTestCentralityHigh(-1)
+{
+  //
+  // Default constructor
+  //
+
+  memset(fEtaRange, 0, sizeof(Double_t) * 2);
+  memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
+}
+//____________________________________________________________
+AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const AliHFECorrectSpectrumBase &ref):
+  TNamed(ref),
+  fCFContainers(NULL),
+  fCorrelation(NULL),
+  fEfficiencyFunction(NULL),
+  fEtaSelected(ref.fEtaSelected),
+  fSetSmoothing(ref.fSetSmoothing),
+  fNbDimensions(ref.fNbDimensions),
+  fNEvents(ref.fNEvents),
+  fStepMC(ref.fStepMC),
+  fStepTrue(ref.fStepTrue),
+  fStepData(ref.fStepData),
+  fStepBeforeCutsV0(ref.fStepBeforeCutsV0),
+  fStepAfterCutsV0(ref.fStepAfterCutsV0),
+  fStepGuessedUnfolding(ref.fStepGuessedUnfolding),
+  fNumberOfIterations(ref.fNumberOfIterations),
+  fChargeChoosen(ref.fChargeChoosen),
+  fTestCentralityLow(ref.fTestCentralityLow),
+  fTestCentralityHigh(ref.fTestCentralityHigh)
+{
+  //
+  // Copy constructor
+  //
+  ref.Copy(*this);
+
+}
+//____________________________________________________________
+AliHFECorrectSpectrumBase &AliHFECorrectSpectrumBase::operator=(const AliHFECorrectSpectrumBase &ref){
+  //
+  // Assignment operator
+  //
+  if(this == &ref) 
+    ref.Copy(*this);
+  return *this;
+}
+//____________________________________________________________
+void AliHFECorrectSpectrumBase::Copy(TObject &o) const {
+  // 
+  // Copy into object o
+  //
+  AliHFECorrectSpectrumBase &target = dynamic_cast<AliHFECorrectSpectrumBase &>(o);
+  target.fCFContainers = fCFContainers;
+  target.fCorrelation = fCorrelation;
+  target.fEfficiencyFunction = fEfficiencyFunction;
+  target.fEtaSelected = fEtaSelected;
+  target.fSetSmoothing = fSetSmoothing;
+  target.fNbDimensions = fNbDimensions;
+  target.fNEvents = fNEvents;
+  target.fStepMC = fStepMC;
+  target.fStepTrue = fStepTrue;
+  target.fStepData = fStepData;
+  target.fStepBeforeCutsV0 = fStepBeforeCutsV0;
+  target.fStepAfterCutsV0 = fStepAfterCutsV0;
+  target.fStepGuessedUnfolding = fStepGuessedUnfolding;
+  target.fNumberOfIterations = fNumberOfIterations;
+  target.fChargeChoosen = fChargeChoosen;
+  target.fTestCentralityLow = fTestCentralityLow;
+  target.fTestCentralityHigh = fTestCentralityHigh;
+  target.fEtaRange[0] = fEtaRange[0];
+  target.fEtaRange[1] = fEtaRange[1];
+  target.fEtaRangeNorm[0] = fEtaRangeNorm[0];
+  target.fEtaRangeNorm[1] = fEtaRangeNorm[1];
+
+}
+
+//____________________________________________________________
+AliHFECorrectSpectrumBase::~AliHFECorrectSpectrumBase(){
+  //
+  // Destructor
+  //
+  if(fCFContainers) delete fCFContainers;
+}
+//__________________________________________________________________________________
+TGraphErrors *AliHFECorrectSpectrumBase::Normalize(THnSparse * const spectrum) const {
+  //
+  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
+  // Give the final pt spectrum to be compared
+  //
+
+  if(fNEvents > 0) {
+
+    TH1D* projection = spectrum->Projection(0);
+    CorrectFromTheWidth(projection);
+    TGraphErrors *graphError = NormalizeTH1(projection);
+    return graphError;
+  
+  }
+    
+  return 0x0;
+  
+
+}
+//__________________________________________________________________________________
+TGraphErrors *AliHFECorrectSpectrumBase::Normalize(AliCFDataGrid * const spectrum) const {
+  //
+  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
+  // Give the final pt spectrum to be compared
+  //
+  if(fNEvents > 0) {
+
+    TH1D* projection = (TH1D *) spectrum->Project(0);
+    CorrectFromTheWidth(projection);
+    TGraphErrors *graphError = NormalizeTH1(projection);
+
+    return graphError;
+    
+  }
+
+  return 0x0;
+  
+
+}
+//__________________________________________________________________________________
+TGraphErrors *AliHFECorrectSpectrumBase::NormalizeTH1(TH1 *input) const {
+  //
+  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
+  // Give the final pt spectrum to be compared
+  //
+  Double_t chargecoefficient = 0.5;
+  if(fChargeChoosen != 0) chargecoefficient = 1.0;
+
+  Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
+  printf("Normalizing Eta Range %f\n", etarange);
+  if(fNEvents > 0) {
+
+    TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
+    Double_t p = 0, dp = 0; Int_t point = 1;
+    Double_t n = 0, dN = 0;
+    Double_t nCorr = 0, dNcorr = 0;
+    Double_t errdN = 0, errdp = 0;
+    for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
+      point = ibin - input->GetXaxis()->GetFirst();
+      p = input->GetXaxis()->GetBinCenter(ibin);
+      dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
+      n = input->GetBinContent(ibin);
+      dN = input->GetBinError(ibin);
+      // New point
+      nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * 1./(2. * TMath::Pi() * p) * n;
+      errdN = 1./(2. * TMath::Pi() * p);
+      errdp = 1./(2. * TMath::Pi() * p*p) * n;
+      dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
+      
+      spectrumNormalized->SetPoint(point, p, nCorr);
+      spectrumNormalized->SetPointError(point, dp, dNcorr);
+    }
+    spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
+    spectrumNormalized->SetMarkerStyle(22);
+    spectrumNormalized->SetMarkerColor(kBlue);
+    spectrumNormalized->SetLineColor(kBlue);
+
+    return spectrumNormalized;
+    
+  }
+
+  return 0x0;
+  
+
+}
+//____________________________________________________________
+void AliHFECorrectSpectrumBase::SetContainer(AliCFContainer *cont, AliHFECorrectSpectrumBase::CFContainer_t type){
+  //
+  // Set the container for a given step to the 
+  //
+  if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
+  fCFContainers->AddAt(cont, type);
+}
+
+//____________________________________________________________
+AliCFContainer *AliHFECorrectSpectrumBase::GetContainer(AliHFECorrectSpectrumBase::CFContainer_t type){
+  //
+  // Get Correction Framework Container for given type
+  //
+  if(!fCFContainers) return NULL;
+  return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
+}
+//____________________________________________________________
+AliCFContainer *AliHFECorrectSpectrumBase::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
+  //
+  // 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()],
+           *varMax = new Double_t[container->GetNVar()];
+
+  Double_t *binLimits;
+  for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
+    
+    binLimits = new Double_t[container->GetNBins(ivar)+1];
+    container->GetBinLimits(ivar,binLimits);
+    varMin[ivar] = binLimits[0];
+    varMax[ivar] = binLimits[container->GetNBins(ivar)];
+    // source
+    if(ivar == 4){
+      if((source>= 0) && (source<container->GetNBins(ivar))) {
+             varMin[ivar] = binLimits[source];
+             varMax[ivar] = binLimits[source];
+      }     
+    }
+    // charge
+    if(ivar == 3) {
+      if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
+    }
+    // eta
+    if(ivar == 1){
+      for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++) 
+        AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic))); 
+      if(fEtaSelected){
+        varMin[ivar] = fEtaRange[0];
+        varMax[ivar] = fEtaRange[1];
+      }
+    }
+    if(fEtaSelected){
+      fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
+      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]));
+    }
+    // centrality
+    if(ivar == 5){
+       if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
+           varMin[ivar] = binLimits[centralitylow];
+           varMax[ivar] = binLimits[centralityhigh];
+
+           TAxis *axistest = container->GetAxis(5,0);
+           AliDebug(1, Form("Number of bin in centrality direction %d\n",axistest->GetNbins()));
+           AliDebug(1, Form("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]));
+           Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
+           Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
+           AliDebug(1, Form("Low centrality %f and high centrality %f\n",lowcentrality,highcentrality));
+       
+       }
+    }
+    
+    delete[] binLimits;
+    
+  }
+  
+  AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
+  delete[] varMin; delete[] varMax;
+
+  return k;
+
+}
+
+//_________________________________________________________________________
+THnSparseF *AliHFECorrectSpectrumBase::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions,Chargetype_t charge,Int_t centralitylow, Int_t centralityhigh) const {
+  //
+  // Slice correlation
+  //
+
+  Int_t ndimensions = correlationmatrix->GetNdimensions();
+  //printf("Number of dimension %d correlation map\n",ndimensions);
+  if(ndimensions < (2*nDim)) {
+    AliError("Problem in the dimensions");
+    return NULL;
+  }
+  
+  // Cut in centrality is centrality > -1
+  if((5+((Int_t)(ndimensions/2.))) < ndimensions) {
+    if((centralitylow >=0) && (centralityhigh >=0)) {
+      
+      TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
+      TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
+      
+      Int_t bins0 = axiscentrality0->GetNbins();
+      Int_t bins1 = axiscentrality1->GetNbins();
+      
+      AliDebug(1, Form("Number of centrality bins: %d and %d\n",bins0,bins1));
+      if(bins0 != bins1) {
+       AliError("Problem in the dimensions");
+       return NULL;
+      }
+      
+      if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
+       axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
+       axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
+       
+       Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
+       Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
+       Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
+       Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
+       AliDebug(1,Form("0 Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0));
+       AliDebug(1,Form("1 Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1));
+       
+      }    
+    }
+  }
+
+  // Cut in eta > -1
+  if(fEtaSelected){
+    if((1+((Int_t)(ndimensions/2.))) < ndimensions) {
+      
+      TAxis *axiseta0 = correlationmatrix->GetAxis(1);
+      TAxis *axiseta1 = correlationmatrix->GetAxis(1+((Int_t)(ndimensions/2.)));
+      
+      Int_t bins0 = axiseta0->GetNbins();
+      Int_t bins1 = axiseta1->GetNbins();
+      
+      AliDebug(1, Form("Number of eta bins: %d and %d\n",bins0,bins1));
+      if(bins0 != bins1) {
+       AliError("Problem in the dimensions");
+       return NULL;
+      }
+      
+      axiseta0->SetRangeUser(fEtaRange[0],fEtaRange[1]);
+      axiseta1->SetRangeUser(fEtaRange[0],fEtaRange[1]);
+       
+       Double_t loweta0 = axiseta0->GetBinLowEdge(axiseta0->FindBin(fEtaRange[0]));
+       Double_t higheta0 = axiseta0->GetBinUpEdge(axiseta0->FindBin(fEtaRange[1]));
+       Double_t loweta1 = axiseta1->GetBinLowEdge(axiseta1->FindBin(fEtaRange[0]));
+       Double_t higheta1 = axiseta1->GetBinUpEdge(axiseta1->FindBin(fEtaRange[1]));
+       AliInfo(Form("0 Low eta %f and high eta %f\n",loweta0,higheta0));
+       AliInfo(Form("1 Low eta %f and high eta %f\n",loweta1,higheta1));
+       
+    }    
+  }
+
+  // Cut in charge
+  if(charge != kAllCharge) {
+    if((3+((Int_t)(ndimensions/2.))) < ndimensions) {
+      
+      TAxis *axischarge0 = correlationmatrix->GetAxis(3);
+      TAxis *axischarge1 = correlationmatrix->GetAxis(3+((Int_t)(ndimensions/2.)));
+      
+      Int_t bins0 = axischarge0->GetNbins();
+      Int_t bins1 = axischarge1->GetNbins();
+      
+      AliDebug(1, Form("Number of charge bins: %d and %d\n",bins0,bins1));
+      if(bins0 != bins1) {
+       AliError("Problem in the dimensions");
+       return NULL;
+      }
+      
+      axischarge0->SetRangeUser(charge,charge);
+      axischarge1->SetRangeUser(charge,charge);
+      
+      Double_t lowcharge0 = axischarge0->GetBinLowEdge(axischarge0->FindBin(charge));
+      Double_t highcharge0 = axischarge0->GetBinUpEdge(axischarge0->FindBin(charge));
+      Double_t lowcharge1 = axischarge1->GetBinLowEdge(axischarge1->FindBin(charge));
+      Double_t highcharge1 = axischarge1->GetBinUpEdge(axischarge1->FindBin(charge));
+      AliInfo(Form("0 Low charge %f and high charge %f\n",lowcharge0,highcharge0));
+      AliInfo(Form("1 Low charge %f and high charge %f\n",lowcharge1,highcharge1));
+      
+    }
+  }
+  
+
+  Int_t ndimensionsContainer = (Int_t) ndimensions/2;
+  
+  Int_t *dim = new Int_t[nDim*2];
+  for(Int_t iter=0; iter < nDim; iter++){
+    dim[iter] = dimensions[iter];
+    dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
+  }
+    
+  THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
+
+  delete[] dim; 
+  return k;
+  
+}
+//___________________________________________________________________________
+void AliHFECorrectSpectrumBase::CorrectFromTheWidth(TH1D *h1) const {
+  //
+  // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
+  //
+
+  TAxis *axis = h1->GetXaxis();
+  Int_t nbinX = h1->GetNbinsX();
+
+  for(Int_t i = 1; i <= nbinX; i++) {
+
+    Double_t width = axis->GetBinWidth(i);
+    Double_t content = h1->GetBinContent(i);
+    Double_t error = h1->GetBinError(i); 
+    h1->SetBinContent(i,content/width); 
+    h1->SetBinError(i,error/width);
+  }
+
+}
+
+//___________________________________________________________________________
+void AliHFECorrectSpectrumBase::CorrectStatErr(AliCFDataGrid *backgroundGrid) const { 
+  //
+  // Correct statistical error
+  //
+
+  TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
+  Int_t nbinX = h1->GetNbinsX();
+  Int_t bins[1];
+  for(Long_t i = 1; i <= nbinX; i++) {
+    bins[0] = i;
+    Float_t content = h1->GetBinContent(i);
+    if(content>0){
+      Float_t error = TMath::Sqrt(content);
+      backgroundGrid->SetElementError(bins, error);
+    }
+  }
+}
+//_________________________________________________________________________
+TObject* AliHFECorrectSpectrumBase::GetSpectrum(const AliCFContainer * const c, Int_t step) {
+  AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
+  return data;
+}
+//_________________________________________________________________________
+TObject* AliHFECorrectSpectrumBase::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
+  // 
+  // Create efficiency grid and calculate efficiency
+  // of step to step0
+  //
+  TString name("eff");
+  name += step;
+  name+= step0;
+  AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
+  eff->CalculateEfficiency(step,step0);
+  return eff;
+}
diff --git a/PWGHF/hfe/AliHFECorrectSpectrumBase.h b/PWGHF/hfe/AliHFECorrectSpectrumBase.h
new file mode 100644 (file)
index 0000000..8ea5fe2
--- /dev/null
@@ -0,0 +1,134 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// For more information see the implementation file
+//
+#ifndef ALIHFECORRECTSPECTRUMBASE_H
+#define ALIHFECORRECTSPECTRUMBASE_H
+
+#ifndef ROOT_TNamed
+#include <TNamed.h>
+#endif
+
+class TGraphErrors;
+class TObject;
+class TH1;
+class TF1;
+class TList;
+class TObjArray;
+class AliCFContainer;
+class AliHFEcontainer;
+class AliCFDataGrid;
+class AliCFEffGrid;
+
+class AliHFECorrectSpectrumBase : public TNamed{
+  public:
+    enum CFContainer_t{
+      kDataContainer  = 0,
+      kBackgroundData = 1,
+      kMCContainerMC = 2,
+      kMCContainerESD = 3,
+      kMCContainerCharmMC = 4,
+      kMCWeightedContainerNonHFEESD =5,
+      kMCWeightedContainerConversionESD = 6,
+      kDataContainerV0 = 7
+    };
+
+    enum Chargetype_t{
+      kNegCharge = -1,
+      kPosCharge = 1,
+      kAllCharge = 0
+    };
+   
+    AliHFECorrectSpectrumBase(const char* name);
+    ~AliHFECorrectSpectrumBase();
+    
+
+    virtual Bool_t Init(const AliHFEcontainer */*datahfecontainer*/, const AliHFEcontainer */*mchfecontainer*/, const AliHFEcontainer */*v0hfecontainer*/, const AliHFEcontainer */*bghfecontainer*/) { return kTRUE;};
+    virtual Bool_t Correct(Bool_t /*subtractcontamination*/) { return kTRUE;};
+   
+    TGraphErrors *Normalize(THnSparse * const spectrum) const;
+    TGraphErrors *Normalize(AliCFDataGrid * const spectrum) const;
+    TGraphErrors *NormalizeTH1(TH1 *input) const;
+    void CorrectFromTheWidth(TH1D *h1) const;
+    void CorrectStatErr(AliCFDataGrid *backgroundGrid) const;
+    
+    void SetCorrelation(THnSparseF * const correlation) {fCorrelation = correlation; };
+    void SetContainer(AliCFContainer *cont, AliHFECorrectSpectrumBase::CFContainer_t type);
+    void SetEfficiencyFunction(TF1 *efficiencyFunction) { fEfficiencyFunction = efficiencyFunction; };
+    
+    void SetNumberOfEvents(Int_t nEvents) { fNEvents = nEvents; };
+    void SetMCEffStep(Int_t step) { fStepMC = step; };
+    void SetMCTruthStep(Int_t step) { fStepTrue = step; };
+    void SetStepToCorrect(Int_t step) { fStepData = step; };
+    void SetStepBeforeCutsV0(Int_t step) { fStepBeforeCutsV0 = step; };
+    void SetStepAfterCutsV0(Int_t step) { fStepAfterCutsV0 = step; };
+    void SetNbDimensions(Int_t nbDimensions) { fNbDimensions = nbDimensions; };
+    void SetChargeChoosen(Chargetype_t chargechoosen) {fChargeChoosen = chargechoosen; };
+    void SetEtaRange(Double_t etamin, Double_t etamax) { fEtaRange[0] = etamin; fEtaRange[1] = etamax; fEtaSelected = kTRUE; }
+    void SetSmoothing(Bool_t setSmoothing) {fSetSmoothing = setSmoothing;};
+    void SetTestOneBinCentrality(Double_t centralitymin, Double_t centralitymax) { fTestCentralityLow = centralitymin; fTestCentralityHigh = centralitymax;}
+    void SetStepGuessedUnfolding(Int_t stepGuessedUnfolding) { fStepGuessedUnfolding = stepGuessedUnfolding; };
+    void SetNumberOfIteration(Int_t numberOfIteration) { fNumberOfIterations = numberOfIteration; };
+    
+    
+
+ protected:
+    AliHFECorrectSpectrumBase(const AliHFECorrectSpectrumBase &ref);
+    AliHFECorrectSpectrumBase &operator=(const AliHFECorrectSpectrumBase &ref);
+    virtual void Copy(TObject &o) const;
+    AliCFContainer *GetContainer(AliHFECorrectSpectrumBase::CFContainer_t contt);
+    AliCFContainer *GetSlicedContainer(AliCFContainer *cont, Int_t ndim, Int_t *dimensions,Int_t source=-1,Chargetype_t charge=kAllCharge,Int_t centralitylow=-1, Int_t centralityhigh=-1);
+    THnSparseF *GetSlicedCorrelation(THnSparseF *correlationmatrix,Int_t nDim, Int_t *dimensions,Chargetype_t charge=kAllCharge,Int_t centralitylow=-1, Int_t centralityhigh=-1) const;
+    TObject* GetSpectrum(const AliCFContainer * const c, Int_t step);
+    TObject* GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0);
+
+    TObjArray *fCFContainers;     // List of Correction Framework Containers
+    THnSparseF *fCorrelation;     // Correlation Matrices
+    TF1 *fEfficiencyFunction;     // Efficiency Function
+   
+    Bool_t fEtaSelected;              // Switch for eta selection
+    Bool_t fSetSmoothing;             // Set smoothing
+
+    Int_t fNbDimensions;          // Number of dimensions for the correction
+    Int_t fNEvents;               // Number of Events
+    Int_t fStepMC;                // MC step (for unfolding)
+    Int_t fStepTrue;              // MC step of the final spectrum
+    Int_t fStepData;              // Data Step (various applications)
+    Int_t fStepBeforeCutsV0;      // Before cuts V0
+    Int_t fStepAfterCutsV0;       // After cuts V0
+    Int_t fStepGuessedUnfolding;  // Step for first guessed unfolding
+    Int_t fNumberOfIterations;    // Number of iterations
+    Chargetype_t fChargeChoosen;         // Select positive or negative electrons
+
+    Double_t fEtaRange[2];        // Eta range 
+    Double_t fEtaRangeNorm[2];    // Eta range used in the normalization
+
+    Int_t fTestCentralityLow;     // To test one bin in centrality only
+    Int_t fTestCentralityHigh;    // To test one bin in centrality only
+      
+
+
+  private:
+   
+   
+    ClassDef(AliHFECorrectSpectrumBase, 1) 
+};
+#endif
+
diff --git a/PWGHF/hfe/AliHFEInclusiveSpectrum.cxx b/PWGHF/hfe/AliHFEInclusiveSpectrum.cxx
new file mode 100644 (file)
index 0000000..4e92fa9
--- /dev/null
@@ -0,0 +1,615 @@
+
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// The following containers have to be set:
+//  - Correction framework container for real data
+//  - Correction framework container for MC (Efficiency Map)
+//  - Correction framework container for background coming from data
+//  - Correction framework container for background coming from MC
+//
+//  Author: 
+//            Raphaelle Bailhache <R.Bailhache@gsi.de>
+//            Markus Fasel <M.Fasel@gsi.de>
+//
+
+#include <TArrayD.h>
+#include <TH1.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TPad.h>
+#include <TH2D.h>
+#include <TF1.h>
+
+#include "AliPID.h"
+#include "AliCFContainer.h"
+#include "AliCFDataGrid.h"
+#include "AliCFEffGrid.h"
+#include "AliCFGridSparse.h"
+#include "AliCFUnfolding.h"
+#include "AliLog.h"
+
+#include "AliHFEInclusiveSpectrum.h"
+#include "AliHFEInclusiveSpectrumQA.h"
+#include "AliHFEcuts.h"
+#include "AliHFEcontainer.h"
+#include "AliHFEtools.h"
+
+ClassImp(AliHFEInclusiveSpectrum)
+
+//____________________________________________________________
+AliHFEInclusiveSpectrum::AliHFEInclusiveSpectrum(const char *name):
+  AliHFECorrectSpectrumBase(name),
+  fQA(NULL)
+{
+  //
+  // Default constructor
+  //
+
+  fQA = new AliHFEInclusiveSpectrumQA("AliHFEInclusiveSpectrumQA");
+}
+//____________________________________________________________
+AliHFEInclusiveSpectrum::AliHFEInclusiveSpectrum(const AliHFEInclusiveSpectrum &ref):
+  AliHFECorrectSpectrumBase(ref),
+  fQA(ref.fQA)
+{
+  //
+  // Copy constructor
+  //
+  ref.Copy(*this);
+
+}
+//____________________________________________________________
+AliHFEInclusiveSpectrum &AliHFEInclusiveSpectrum::operator=(const AliHFEInclusiveSpectrum &ref){
+  //
+  // Assignment operator
+  //
+  if(this == &ref) 
+    ref.Copy(*this);
+  return *this;
+}
+//____________________________________________________________
+void AliHFEInclusiveSpectrum::Copy(TObject &o) const {
+  // 
+  // Copy into object o
+  //
+  AliHFEInclusiveSpectrum &target = dynamic_cast<AliHFEInclusiveSpectrum &>(o);
+  target.fQA = fQA;
+
+
+}
+//____________________________________________________________
+AliHFEInclusiveSpectrum::~AliHFEInclusiveSpectrum(){
+  //
+  // Destructor
+  //
+  if(fQA) delete fQA;
+  
+}
+//____________________________________________________________
+Bool_t AliHFEInclusiveSpectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer */*bghfecontainer*/){
+  //
+  // Init what we need for the correction:
+  //
+  // Raw spectrum, hadron contamination
+  // MC efficiency maps, correlation matrix
+  // V0 efficiency if wanted
+  //
+  // This for a given dimension.
+  //
+  //
+
+  Int_t kNdim = 3;
+
+  Int_t dims[kNdim];
+  switch(fNbDimensions){
+  case 1:   dims[0] = 0;
+    break;
+  case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
+    break;
+  case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
+    break;
+  default:
+    AliError("Container with this number of dimensions not foreseen (yet)");
+    return kFALSE;
+  };
+  
+  //
+  // Data container: raw spectrum + hadron contamination  
+  //  
+  AliCFContainer *datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
+  AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
+  if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
+  AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
+  SetContainer(datacontainerD,AliHFECorrectSpectrumBase::kDataContainer);
+  SetContainer(contaminationcontainerD,AliHFECorrectSpectrumBase::kBackgroundData);
+  // QA 
+  Int_t dimqa = datacontainer->GetNVar();
+  Int_t dimsqa[dimqa];
+  for(Int_t i = 0; i < dimqa; i++) dimsqa[i] = i;
+  AliCFContainer *datacontainerDQA = GetSlicedContainer(datacontainer, dimqa, dimsqa, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  fQA->AddResultAt(datacontainerDQA,AliHFEInclusiveSpectrumQA::kDataProjection);
+
+  //
+  // MC container: ESD/MC efficiency maps + MC/MC efficiency maps 
+  //
+  AliCFContainer *mccontaineresd = 0x0;
+  AliCFContainer *mccontainermc = 0x0;
+  mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
+  mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
+  if((!mccontaineresd) || (!mccontainermc)) return kFALSE;  
+  Int_t source = -1;
+  AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
+  SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerMC);
+  SetContainer(mccontaineresdD,AliHFECorrectSpectrumBase::kMCContainerESD);
+
+  //
+  // Correlation matrix
+  //  
+  THnSparseF *mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
+  if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
+  if(!mccorrelation) return kFALSE;
+  THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  if(!mccorrelationD) {
+    printf("No correlation\n");
+    return kFALSE;
+  }
+  SetCorrelation(mccorrelationD);
+  // QA
+  fQA->AddResultAt(mccorrelation,AliHFEInclusiveSpectrumQA::kCMProjection);
+  //
+  // V0 container Electron, pt eta phi
+  //
+  if(v0hfecontainer) {
+    AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
+    if(!containerV0) return kFALSE;
+    AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+    if(!containerV0Electron) return kFALSE;
+    SetContainer(containerV0Electron,AliHFECorrectSpectrumBase::kDataContainerV0);
+  }
+
+  //
+  fQA->DrawProjections();
+
+  
+  return kTRUE;
+}
+//____________________________________________________________
+Bool_t AliHFEInclusiveSpectrum::Correct(Bool_t subtractcontamination){
+  //
+  // Correct the spectrum for efficiency and unfolding
+  // with both method and compare
+  //
+  
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(1111);
+  gStyle->SetPadBorderMode(0);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetPadLeftMargin(0.13);
+  gStyle->SetPadRightMargin(0.13);
+
+  printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
+
+  ///////////////////////////
+  // Check initialization
+  ///////////////////////////
+
+  if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
+    AliInfo("You have to init before");
+    return kFALSE;
+  }
+  
+  if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
+    AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
+    return kFALSE;
+  }
+  SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
+    
+  AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
+  //////////////////////////////////
+  // Subtract hadron background
+  /////////////////////////////////
+  AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
+  if(subtractcontamination) {
+    dataspectrumaftersubstraction = SubtractBackground();
+    dataGridAfterFirstSteps = dataspectrumaftersubstraction;
+  }
+
+  ////////////////////////////////////////////////
+  // Correct for TPC efficiency from V0 if any
+  ///////////////////////////////////////////////
+  AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
+  AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
+  if(dataContainerV0){
+    dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
+    dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;  
+  }
+
+  //////////////////////////////////////////////////////////////////////////////
+  // Correct for efficiency parametrized (if TPC efficiency is parametrized)
+  /////////////////////////////////////////////////////////////////////////////
+  AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
+  if(fEfficiencyFunction){
+    dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
+    dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;  
+  }
+    
+  ///////////////
+  // Unfold
+  //////////////
+  THnSparse *correctedspectrum = Unfold(dataGridAfterFirstSteps);
+  if(!correctedspectrum){
+    AliError("No corrected spectrum\n");
+    return kFALSE;
+  }
+  
+  /////////////////////
+  // Simply correct
+  ////////////////////
+  AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
+
+  // QA final results
+  TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
+  TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
+  fQA->AddResultAt(correctedspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
+  fQA->AddResultAt(alltogetherspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
+  fQA->DrawResult();
+
+  return kTRUE;
+}
+
+//____________________________________________________________
+AliCFDataGrid* AliHFEInclusiveSpectrum::SubtractBackground(){
+  //
+  // Apply background subtraction
+  //
+
+  // Raw spectrum
+  AliCFContainer *dataContainer = GetContainer(kDataContainer);
+  if(!dataContainer){
+    AliError("Data Container not available");
+    return NULL;
+  }
+  printf("Step data: %d\n",fStepData);
+  AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
+
+  AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
+  dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction"); 
+
+  // Background Estimate
+  AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
+  if(!backgroundContainer){
+    AliError("MC background container not found");
+    return NULL;
+  }
+  Int_t stepbackground = 1; 
+  AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
+
+  // Subtract 
+  spectrumSubtracted->Add(backgroundGrid,-1.0);
+  
+  // QA
+  TH1D *subtractedspectrum = (TH1D *) spectrumSubtracted->Project(0);
+  CorrectFromTheWidth(subtractedspectrum);
+  TH1D *rawspectrum = (TH1D *) dataspectrumbeforesubstraction->Project(0);
+  CorrectFromTheWidth(rawspectrum);
+  fQA->AddResultAt(subtractedspectrum,AliHFEInclusiveSpectrumQA::kAfterSC);
+  fQA->AddResultAt(rawspectrum,AliHFEInclusiveSpectrumQA::kBeforeSC);
+  fQA->DrawSubtractContamination();
+
+  return spectrumSubtracted;
+}
+
+//____________________________________________________________
+AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Apply TPC pid efficiency correction from parametrisation
+  //
+
+  // Data in the right format
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    dataGrid = bgsubpectrum;
+  }
+  else {
+    
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      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");
+    return NULL;
+  }
+  AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
+  dataContainerbis->Add(dataContainerbis,-1.0);
+
+
+  Int_t* coord = new Int_t[nbdimensions];
+  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) {
+    
+    Double_t value = h->GetBinContent(i, coord);
+    //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
+    //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
+    
+    // Get the bin co-ordinates given an coord
+    for (Int_t j = 0; j < nbdimensions; ++j)
+      points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
+
+    if (!fEfficiencyFunction->IsInside(points))
+         continue;
+    TF1::RejectPoint(kFALSE);
+
+    // Evaulate function at points
+    Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
+    //printf("Value efficiency is %f\n",valueEfficiency);
+
+    if(valueEfficiency > 0.0) {
+      h->SetBinContent(coord,value/valueEfficiency);
+      dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
+    }
+    Double_t error = h->GetBinError(i);
+    h->SetBinError(coord,error/valueEfficiency);
+    dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
+
+   
+  } 
+
+  delete[] coord;
+  delete[] points;
+
+  AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
+
+ // QA
+  TH1D *afterE = (TH1D *) resultt->Project(0);
+  CorrectFromTheWidth(afterE);
+  TH1D *beforeE = (TH1D *) dataGrid->Project(0);
+  CorrectFromTheWidth(beforeE);
+  fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterPE);
+  fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforePE);
+  fQA->AddResultAt(fEfficiencyFunction,AliHFEInclusiveSpectrumQA::kPEfficiency);
+  fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kParametrized);
+  
+  return resultt;
+
+}
+//____________________________________________________________
+AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Apply TPC pid efficiency correction from V0
+  //
+
+  AliCFContainer *v0Container = GetContainer(kDataContainerV0);
+  if(!v0Container){
+    AliError("V0 Container not available");
+    return NULL;
+  }
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
+  efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
+
+  // Data in the right format
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    dataGrid = bgsubpectrum;
+  }
+  else {
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
+  } 
+
+  // Correct
+  AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
+  result->ApplyEffCorrection(*efficiencyD);
+
+  // QA
+  TH1D *afterE = (TH1D *) result->Project(0);
+  CorrectFromTheWidth(afterE);
+  TH1D *beforeE = (TH1D *) dataGrid->Project(0);
+  CorrectFromTheWidth(beforeE);
+  TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
+  fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterV0);
+  fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeV0);
+  fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kV0Efficiency);
+  fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kV0);
+
+  return result;
+
+}
+//____________________________________________________________
+THnSparse *AliHFEInclusiveSpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Return the unfolded spectrum
+  //
+
+  AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return NULL;
+  }
+
+  if(!fCorrelation){
+    AliError("No Correlation map available");
+    return NULL;
+  }
+
+  // Data 
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    dataGrid = bgsubpectrum;
+  }
+  else {
+
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
+  } 
+  
+  // Guessed
+  AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
+  THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
+  efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
+
+  // Unfold 
+  
+  AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);
+  if(fSetSmoothing) unfolding.UseSmoothing();
+  unfolding.Unfold();
+
+  // Results
+  THnSparse* result = unfolding.GetUnfolded();
+  THnSparse* residual = unfolding.GetEstMeasured();
+  
+  // QA
+  TH1D *residualh = (TH1D *) residual->Projection(0);
+  TH1D *beforeE = (TH1D *) dataGrid->Project(0);
+  TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
+  TH1D *afterE = (TH1D *) result->Projection(0);
+  CorrectFromTheWidth(residualh);
+  CorrectFromTheWidth(beforeE);
+  CorrectFromTheWidth(afterE);
+  fQA->AddResultAt(residualh,AliHFEInclusiveSpectrumQA::kResidualU);
+  fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterU);
+  fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeU);
+  fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kUEfficiency);
+  fQA->DrawUnfolding();
+
+  return (THnSparse *) result->Clone();
+
+}
+//____________________________________________________________
+AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Apply unfolding and efficiency correction together to bgsubspectrum
+  //
+
+  AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return NULL;
+  }
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
+  efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
+
+  // Data in the right format
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    dataGrid = bgsubpectrum;
+  }
+  else {
+    
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
+  } 
+
+  // Correct
+  AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
+  result->ApplyEffCorrection(*efficiencyD);
+
+  // QA
+  TH1D *afterE = (TH1D *) result->Project(0);
+  CorrectFromTheWidth(afterE);
+  TH1D *beforeE = (TH1D *) dataGrid->Project(0);
+  CorrectFromTheWidth(beforeE);
+  TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
+  fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterMCE);
+  fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeMCE);
+  fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kMCEfficiency);
+  fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kMC);
+
+  
+  return result;
+
+}
+//____________________________________________________________
+void AliHFEInclusiveSpectrum::WriteResults(const char *filename)
+{
+  //
+  // Write results
+  //
+
+  AliCFContainer *dataContainer = GetContainer(kDataContainer);
+  AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
+  TObject *unfolded = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
+  TObject *correctedspectrum = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
+
+  TFile *file = TFile::Open(filename,"recreate");
+  if(dataContainer) dataContainer->Write("data");
+  if(mcContainer) mcContainer->Write("mcefficiency");
+  if(fCorrelation) fCorrelation->Write("correlationmatrix");
+  if(unfolded) unfolded->Write("unfoldedspectrum");
+  if(correctedspectrum) correctedspectrum->Write("correctedspectrum");
+  if(fQA) fQA->Write("QAResults");
+  file->Close();
+
+}
+
diff --git a/PWGHF/hfe/AliHFEInclusiveSpectrum.h b/PWGHF/hfe/AliHFEInclusiveSpectrum.h
new file mode 100644 (file)
index 0000000..88fe135
--- /dev/null
@@ -0,0 +1,67 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// For more information see the implementation file
+//
+#ifndef ALIHFEINCLUSIVESPECTRUM_H
+#define ALIHFEINCLUSIVESPECTRUM_H
+
+#include "AliHFECorrectSpectrumBase.h"
+
+
+class TGraphErrors;
+class TObject;
+class TH1;
+class TF1;
+class TList;
+class TObjArray;
+class AliCFContainer;
+class AliHFEcontainer;
+class AliCFDataGrid;
+class AliCFEffGrid;
+class AliHFEInclusiveSpectrumQA;
+
+class AliHFEInclusiveSpectrum : public AliHFECorrectSpectrumBase{
+  public:
+  
+  AliHFEInclusiveSpectrum(const char* name);
+    ~AliHFEInclusiveSpectrum();
+    
+
+    virtual Bool_t Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer=0x0, const AliHFEcontainer */*bghfecontainer*/=0x0);
+    virtual Bool_t Correct(Bool_t subtractcontamination=kTRUE);
+   
+    AliCFDataGrid *SubtractBackground();
+    AliCFDataGrid *CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum = 0x0);
+    AliCFDataGrid *CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum = 0x0);
+    THnSparse *Unfold(AliCFDataGrid* const bgsubpectrum = 0x0);
+    AliCFDataGrid *CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum = 0x0);
+
+    void WriteResults(const char *filename);
+   
+ private:
+    AliHFEInclusiveSpectrum(const AliHFEInclusiveSpectrum &ref);
+    AliHFEInclusiveSpectrum &operator=(const AliHFEInclusiveSpectrum &ref);
+    virtual void Copy(TObject &o) const;
+    AliHFEInclusiveSpectrumQA *fQA; // QA
+   
+    ClassDef(AliHFEInclusiveSpectrum, 1) 
+};
+#endif
+
diff --git a/PWGHF/hfe/AliHFEInclusiveSpectrumQA.cxx b/PWGHF/hfe/AliHFEInclusiveSpectrumQA.cxx
new file mode 100644 (file)
index 0000000..63d181c
--- /dev/null
@@ -0,0 +1,562 @@
+
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// The following containers have to be set:
+//  - Correction framework container for real data
+//  - Correction framework container for MC (Efficiency Map)
+//  - Correction framework container for background coming from data
+//  - Correction framework container for background coming from MC
+//
+//  Author: 
+//            Raphaelle Bailhache <R.Bailhache@gsi.de>
+//
+
+#include <TArrayD.h>
+#include <TH1.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TObject.h>
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TPad.h>
+#include <TH2D.h>
+#include <TF1.h>
+
+#include "AliPID.h"
+#include "AliCFContainer.h"
+#include "AliCFDataGrid.h"
+#include "AliCFEffGrid.h"
+#include "AliCFGridSparse.h"
+#include "AliCFUnfolding.h"
+#include "AliLog.h"
+
+#include "AliHFEInclusiveSpectrumQA.h"
+#include "AliHFECorrectSpectrumBase.h"
+#include "AliHFEcuts.h"
+#include "AliHFEcontainer.h"
+#include "AliHFEtools.h"
+
+ClassImp(AliHFEInclusiveSpectrumQA)
+
+const Char_t *AliHFEInclusiveSpectrumQA::fgkNameCanvas[AliHFEInclusiveSpectrumQA::kNTypeEfficiency] = {
+  "V0Efficiency",
+  "MCEfficiency",
+  "ParametrizedEfficiency"
+};
+
+//____________________________________________________________
+AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA():
+  TNamed(),
+  fPtMax(7.0),
+  fListOfResult(),
+  fWriteToFile(kTRUE)
+{
+  //
+  // Default constructor
+  //
+
+  fListOfResult = new TObjArray(kNResults);
+  fListOfResult->SetName("ListOfResults");
+
+}
+//____________________________________________________________
+AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA(const char *name):
+  TNamed(name, ""),
+  fPtMax(7.0),
+  fListOfResult(),
+  fWriteToFile(kTRUE)
+{
+  //
+  // Default constructor
+  //
+
+  fListOfResult = new TObjArray(kNResults);
+  fListOfResult->SetName("ListOfResults");
+
+}
+
+//____________________________________________________________
+AliHFEInclusiveSpectrumQA::~AliHFEInclusiveSpectrumQA(){
+  //
+  // Destructor
+  //
+  if(fListOfResult) delete fListOfResult;
+}
+//____________________________________________________________
+void AliHFEInclusiveSpectrumQA::AddResultAt(TObject *obj,Int_t index)
+{
+  //
+  // Init what we need for the correction:
+  //
+
+  if(fListOfResult) fListOfResult->AddAt(obj,index);
+
+}
+//____________________________________________________________
+TObject *AliHFEInclusiveSpectrumQA::GetResult(Int_t index)
+{
+  //
+  // Get result
+  //
+
+  if(fListOfResult) return fListOfResult->UncheckedAt(index);
+  else return 0x0;
+
+}
+//____________________________________________________________
+void AliHFEInclusiveSpectrumQA::DrawProjections() const
+{
+  //
+  // get spectrum for beauty 2nd method
+  //
+  //
+  AliCFContainer *data = (AliCFContainer *) fListOfResult->UncheckedAt(kDataProjection);
+  THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
+  if(!data || !correlation) return;
+
+  Int_t ndimcont = data->GetNVar();
+  Int_t ndimcor = correlation->GetNdimensions();
+  Int_t charge = 3;
+  Int_t centrality = 5;
+  Int_t eta = 1;
+
+  TCanvas * canvas = new TCanvas("Projections","Projections",1000,700);
+  Int_t n = 0;
+  if(charge < ndimcont) n++;
+  if(centrality < ndimcont) n++;
+  if(eta < ndimcont) n++;
+  canvas->Divide(2,n);
+  Int_t counter = 1;
+
+  if(charge < ndimcont) {
+   
+    canvas->cd(counter);
+    TH1 *checkcharge = (TH1 *) data->Project(data->GetNStep()-1,charge);
+    checkcharge->Draw();
+    counter++;
+    canvas->cd(counter);
+    TH2F* projectioncharge = (TH2F *) correlation->Projection(charge,charge+((Int_t)(ndimcor/2.)));
+    projectioncharge->Draw("colz");
+    counter++;  
+
+  }
+
+  if(centrality < ndimcont) {
+    canvas->cd(counter);
+    TH1 *checkcentrality = (TH1 *) data->Project(data->GetNStep()-1,centrality);
+    checkcentrality->Draw();
+    counter++;
+    canvas->cd(counter);
+    TH2F *projectioncentrality = (TH2F *) correlation->Projection(centrality,centrality+((Int_t)(ndimcor/2.)));
+    projectioncentrality->Draw("colz");
+    counter++;  
+  }
+
+  if(eta < ndimcont) {
+    canvas->cd(counter);
+    TH1 *checketa = (TH1 *) data->Project(data->GetNStep()-1,eta);
+    checketa->Draw();
+    counter++;
+    canvas->cd(counter);
+    TH2D* projectioneta = (TH2D *) correlation->Projection(eta,eta+((Int_t)(ndimcor/2.)));
+    projectioneta->Draw("colz");
+  }
+
+
+}
+//____________________________________________________________
+void AliHFEInclusiveSpectrumQA::DrawSubtractContamination() const
+{
+  //
+  // get spectrum for beauty 2nd method
+  //
+  //
+  TH1D *measuredTH1Daftersubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSC);
+  TH1D *measuredTH1Dbeforesubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSC);
+  if(!measuredTH1Daftersubstraction || !measuredTH1Dbeforesubstraction) return;
+
+  SetStyle();
+
+  TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
+  cbackgroundsubtraction->Divide(2,1);
+  cbackgroundsubtraction->cd(1);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  measuredTH1Daftersubstraction->SetStats(0);
+  measuredTH1Daftersubstraction->SetTitle("");
+  measuredTH1Daftersubstraction->GetYaxis()->SetTitleOffset(1.5);
+  measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  measuredTH1Daftersubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  measuredTH1Daftersubstraction->SetMarkerStyle(25);
+  measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
+  measuredTH1Daftersubstraction->SetLineColor(kBlack);
+  measuredTH1Dbeforesubstraction->SetStats(0);
+  measuredTH1Dbeforesubstraction->SetTitle("");
+  measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  measuredTH1Dbeforesubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
+  measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
+  measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
+  measuredTH1Daftersubstraction->Draw();
+  measuredTH1Dbeforesubstraction->Draw("same");
+  TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
+  legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
+  legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
+  legsubstraction->SetFillStyle(0);
+  legsubstraction->SetLineStyle(0);
+  legsubstraction->SetLineColor(0);
+  legsubstraction->Draw("same");
+  cbackgroundsubtraction->cd(2);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
+  ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
+  ratiomeasuredcontamination->SetTitle("");
+  ratiomeasuredcontamination->GetYaxis()->SetTitleOffset(1.5);
+  ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
+  ratiomeasuredcontamination->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  ratiomeasuredcontamination->GetYaxis()->SetRangeUser(0.8,1.2);
+  ratiomeasuredcontamination->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  ratiomeasuredcontamination->Sumw2();
+  ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
+  ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
+  ratiomeasuredcontamination->SetStats(0);
+  ratiomeasuredcontamination->SetMarkerStyle(26);
+  ratiomeasuredcontamination->SetMarkerColor(kBlack);
+  ratiomeasuredcontamination->SetLineColor(kBlack);
+  for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
+    ratiomeasuredcontamination->SetBinError(k+1,0.0);
+  }
+  ratiomeasuredcontamination->Draw("P");
+  if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.png");
+
+}
+
+
+//____________________________________________________________
+void AliHFEInclusiveSpectrumQA::DrawCorrectWithEfficiency(Int_t typeeff) const
+{
+  //
+  // Correct the spectrum for efficiency and unfolding
+  // with both method and compare
+  //
+  
+  TH1D *afterE = 0x0;
+  TH1D *beforeE = 0x0;
+  TH1D *efficiencyDproj = 0x0;
+  TF1 *efficiencyparametrized = 0x0;
+
+  if(typeeff== kV0) {
+    afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterV0);
+    beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeV0);
+    efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kV0Efficiency);
+  }
+  if(typeeff== kMC) {
+    afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
+    beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeMCE);
+    efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kMCEfficiency);
+  }
+ if(typeeff== kParametrized) {
+    afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterPE);
+    beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforePE);
+    efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiency);
+  }
+
+ if((typeeff==kV0 || typeeff==kMC) && (!afterE || !beforeE || !efficiencyDproj)) return;
+ if(typeeff==kParametrized && (!afterE || !beforeE || !efficiencyparametrized)) return;
+
+  SetStyle();
+
+  TCanvas * cEfficiency = new TCanvas(AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],1000,700);
+  cEfficiency->Divide(2,1);
+  cEfficiency->cd(1);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  afterE->SetStats(0);
+  afterE->SetTitle("");
+  afterE->GetYaxis()->SetTitleOffset(1.5);
+  afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  afterE->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  afterE->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  afterE->SetMarkerStyle(25);
+  afterE->SetMarkerColor(kBlack);
+  afterE->SetLineColor(kBlack);
+  beforeE->SetStats(0);
+  beforeE->SetTitle("");
+  beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  beforeE->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  beforeE->SetMarkerStyle(24);
+  beforeE->SetMarkerColor(kBlue);
+  beforeE->SetLineColor(kBlue);
+  afterE->Draw();
+  beforeE->Draw("same");
+  TLegend *legefficiency = new TLegend(0.4,0.6,0.89,0.89);
+  legefficiency->AddEntry(beforeE,"Before Efficiency correction","p");
+  legefficiency->AddEntry(afterE,"After Efficiency correction","p");
+  legefficiency->SetFillStyle(0);
+  legefficiency->SetLineStyle(0);
+  legefficiency->SetLineColor(0);
+  legefficiency->Draw("same");
+  cEfficiency->cd(2);
+  gPad->SetTicks();
+  if((typeeff==kV0 || typeeff==kMC)) {
+    efficiencyDproj->SetTitle("");
+    efficiencyDproj->SetStats(0);
+    efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
+    efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
+    efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
+    efficiencyDproj->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+    efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
+    efficiencyDproj->SetMarkerStyle(25);
+    efficiencyDproj->Draw();
+  }
+  if(typeeff==kParametrized) {
+    efficiencyparametrized->Draw();
+  }
+  
+  if(fWriteToFile) {
+    if(typeeff==kV0) cEfficiency->SaveAs("EfficiencyV0.png");
+    if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMC.png");
+    if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrized.png");
+  }
+
+}
+
+//____________________________________________________________
+void AliHFEInclusiveSpectrumQA::DrawUnfolding() const
+{
+  //
+  // Draw unfolding
+  //
+  TH1D *measuredspectrumD = (TH1D *) fListOfResult->UncheckedAt(kBeforeU);
+  TH1D *residualspectrumD = (TH1D *) fListOfResult->UncheckedAt(kResidualU);
+  TH1D *efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kUEfficiency);
+  THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
+  
+  if(!measuredspectrumD || !residualspectrumD || !efficiencyDproj || !correlation) return;
+  
+  Int_t ndimcor = (Int_t) correlation->GetNdimensions()/2.;
+
+  SetStyle();
+
+  TCanvas * cunfolding = new TCanvas("unfolding","unfolding",1000,700);
+  cunfolding->Divide(2,2);
+  cunfolding->cd(1);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  residualspectrumD->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  residualspectrumD->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  residualspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  residualspectrumD->SetStats(0);
+  residualspectrumD->SetTitle("");
+  residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  residualspectrumD->SetMarkerStyle(26);
+  residualspectrumD->SetMarkerColor(kBlue);
+  residualspectrumD->SetLineColor(kBlue);
+  residualspectrumD->Sumw2();
+  residualspectrumD->Draw("P");
+  measuredspectrumD->SetStats(0);
+  measuredspectrumD->SetTitle("");  
+  measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  measuredspectrumD->SetMarkerStyle(25);
+  measuredspectrumD->SetMarkerColor(kBlack);
+  measuredspectrumD->SetLineColor(kBlack);
+  measuredspectrumD->Draw("same");
+  TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
+  legres->AddEntry(residualspectrumD,"Residual","p");
+  legres->AddEntry(measuredspectrumD,"Measured","p");
+  legres->SetFillStyle(0);
+  legres->SetLineStyle(0);
+  legres->SetLineColor(0);
+  legres->Draw("same");
+  cunfolding->cd(2);
+  gPad->SetTicks();
+  TH1D* ratioresidual = (TH1D*)residualspectrumD->Clone();
+  ratioresidual->SetName("ratioresidual");
+  ratioresidual->SetTitle("");
+  ratioresidual->GetYaxis()->SetRangeUser(0.6,1.4);
+  ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
+  ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  ratioresidual->Divide(measuredspectrumD);
+  ratioresidual->SetStats(0);
+  ratioresidual->Draw();
+  cunfolding->cd(3);
+  gPad->SetTicks();
+  efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
+  efficiencyDproj->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
+  efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
+  efficiencyDproj->SetStats(0);
+  efficiencyDproj->SetTitle("");
+  efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
+  efficiencyDproj->SetMarkerStyle(26);
+  efficiencyDproj->SetMarkerColor(kBlue);
+  efficiencyDproj->SetLineColor(kBlue);
+  efficiencyDproj->Sumw2();
+  efficiencyDproj->Draw("P");
+  cunfolding->cd(4);
+  TH2F *projectioncorr = (TH2F *) correlation->Projection(0,ndimcor);
+  projectioncorr->GetYaxis()->SetTitle("p^{ESD}_{T} [GeV/c]");
+  projectioncorr->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
+  projectioncorr->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  projectioncorr->GetYaxis()->SetRangeUser(0.0,fPtMax);
+  projectioncorr->SetStats(0);
+  projectioncorr->SetTitle("");
+  projectioncorr->Draw("colz");
+
+  if(fWriteToFile){
+    cunfolding->SaveAs("Unfolding.png");
+  }
+  
+}
+//____________________________________________________________
+void AliHFEInclusiveSpectrumQA::DrawResult()
+{
+  //
+  // Draw Results
+  //
+  TGraphErrors* correctedspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultUnfolded);
+  TGraphErrors* alltogetherspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt( kFinalResultDirectEfficiency);
+  if(!correctedspectrumD || !alltogetherspectrumD) return;
+  
+  SetStyle();
+
+  TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
+  ccorrected->Divide(2,1);
+  ccorrected->cd(1);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  correctedspectrumD->SetTitle("");
+  correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  correctedspectrumD->GetYaxis()->SetRangeUser(0.000001,100.0);
+  correctedspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  correctedspectrumD->SetMarkerStyle(26);
+  correctedspectrumD->SetMarkerColor(kBlue);
+  correctedspectrumD->SetLineColor(kBlue);
+  correctedspectrumD->Draw("AP");
+  alltogetherspectrumD->SetTitle("");
+  alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
+  alltogetherspectrumD->SetMarkerStyle(25);
+  alltogetherspectrumD->SetMarkerColor(kBlack);
+  alltogetherspectrumD->SetLineColor(kBlack);
+  alltogetherspectrumD->Draw("P");
+  TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
+  legcorrected->AddEntry(correctedspectrumD,"Unfolded","p");
+  legcorrected->AddEntry(alltogetherspectrumD,"Direct corrected","p");
+  legcorrected->SetFillStyle(0);
+  legcorrected->SetLineStyle(0);
+  legcorrected->SetLineColor(0);
+  legcorrected->Draw("same");
+  ccorrected->cd(2);
+  gPad->SetTicks();
+  TH1D* ratiocorrected = DivideSpectra(correctedspectrumD,alltogetherspectrumD);
+  ratiocorrected->SetName("ratiocorrected");
+  ratiocorrected->SetTitle("");
+  ratiocorrected->GetYaxis()->SetTitleOffset(1.5);
+  ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
+  ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  ratiocorrected->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  ratiocorrected->GetYaxis()->SetRangeUser(0.4,1.4);
+  ratiocorrected->SetStats(0);
+  ratiocorrected->Draw();
+  if(fWriteToFile)ccorrected->SaveAs("CorrectedResults.png");
+
+}
+//____________________________________________________________
+TH1D *AliHFEInclusiveSpectrumQA::DivideSpectra(TGraphErrors *ga, TGraphErrors *gb) 
+{
+  //
+  // Divide Spectra
+  //
+
+  TH1D *afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
+  if(!afterE) return 0x0;
+
+  TH1D *histoB = (TH1D*) afterE->Clone();
+  histoB->Sumw2();
+  histoB->SetName("ratio");
+  TH1D *histoa = (TH1D*) afterE->Clone();
+  histoa->Sumw2();
+  histoa->SetName("a");
+  TH1D *histob = (TH1D*) afterE->Clone();
+  histob->Sumw2();
+  histob->SetName("b");
+  
+  double xa,ya,xb,yb,eya,eyb;
+  Int_t npointsa = ga->GetN();
+  Int_t npointsb = gb->GetN();
+  if(npointsa != npointsb) {
+    printf("Problem the two spectra have not the same number of points\n");
+    return 0x0;
+  }
+  for(Int_t k = 0; k < npointsa; k++){
+    ga->GetPoint(k,xa,ya);
+    gb->GetPoint(k,xb,yb);
+    //
+    Double_t centerhisto = histoa->GetBinCenter(k+1);
+    //
+    if((TMath::Abs(xa-xb) > 0.0001) || (TMath::Abs(xa-centerhisto) > 0.0001)) {
+      printf("Problem not the same x axis\n");
+      return 0x0;
+    }
+    histoa->SetBinContent(k+1,ya);
+    histob->SetBinContent(k+1,yb);
+    //
+    eya = ga->GetErrorY(k);
+    eyb = gb->GetErrorY(k);
+    //
+    histoa->SetBinError(k+1,eya);
+    histob->SetBinError(k+1,eyb);
+   
+  }
+  
+  histoB->Sumw2();
+  histoB->Divide(histoa,histob,1.0,1.0,"B");
+
+  return histoB;  
+
+}
+//__________________________________________
+void AliHFEInclusiveSpectrumQA::SetStyle() const
+{
+  //
+  // Set style
+  //
+
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(1111);
+  gStyle->SetPadBorderMode(0);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetPadLeftMargin(0.13);
+  gStyle->SetPadRightMargin(0.13);
+
+}
diff --git a/PWGHF/hfe/AliHFEInclusiveSpectrumQA.h b/PWGHF/hfe/AliHFEInclusiveSpectrumQA.h
new file mode 100644 (file)
index 0000000..9bad0f4
--- /dev/null
@@ -0,0 +1,105 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// For more information see the implementation file
+//
+#ifndef ALIHFEINCLUSIVESPECTRUMQA_H
+#define ALIHFEINCLUSIVESPECTRUMQA_H
+
+#ifndef ROOT_TNamed
+#include <TNamed.h>
+#endif
+
+
+class TObjArray;
+class TGraphErrors;
+class TObject;
+
+
+class AliHFEInclusiveSpectrumQA : public TNamed{
+ public:
+
+  enum Results_t{
+    kDataProjection = 0,
+    kCMProjection = 1,
+    kBeforeSC = 2,
+    kAfterSC = 3,
+    kBeforeV0 = 4,
+    kAfterV0 = 5,
+    kV0Efficiency = 6,
+    kBeforePE = 7,
+    kAfterPE = 8,
+    kPEfficiency = 9,
+    kBeforeMCE = 10,
+    kAfterMCE = 11,
+    kMCEfficiency = 12,
+    kBeforeU = 13,
+    kAfterU = 14,
+    kResidualU = 15,
+    kUEfficiency = 16,
+    kFinalResultUnfolded = 17,
+    kFinalResultDirectEfficiency = 18,
+    kNResults = 19
+  };
+
+
+  enum EfficiencyCorrection_t{
+    kV0 = 0,
+    kMC = 1,
+    kParametrized = 2,
+    kNTypeEfficiency = 3
+  };
+   
+  void AddResultAt(TObject *obj,Int_t index);
+  TObject *GetResult(Int_t index);
+  
+  void DrawProjections() const;
+  void DrawSubtractContamination() const;
+  void DrawCorrectWithEfficiency(Int_t typeeff) const;
+  void DrawUnfolding() const;
+  void DrawResult();
+  
+  void SetStyle() const;
+  void SetWriteToFile(Bool_t writetofile) {fWriteToFile=writetofile; };
+  void SetPtMax(Double_t ptmax) {fPtMax = ptmax; };
+
+
+  TH1D *DivideSpectra(TGraphErrors *ga, TGraphErrors *gb);
+  
+
+  AliHFEInclusiveSpectrumQA();
+  AliHFEInclusiveSpectrumQA(const char* name);
+  ~AliHFEInclusiveSpectrumQA();
+  
+  protected:
+  
+  private:
+
+  static const Char_t* fgkNameCanvas[kNTypeEfficiency];     // Name of canvas
+
+    AliHFEInclusiveSpectrumQA(const AliHFEInclusiveSpectrumQA &);
+    AliHFEInclusiveSpectrumQA &operator=(const AliHFEInclusiveSpectrumQA &);
+    Double_t   fPtMax;             // Pt max to plot
+    TObjArray *fListOfResult;     // ListOfResults
+    Bool_t fWriteToFile;           // Write plots to eps files
+   
+    ClassDef(AliHFEInclusiveSpectrumQA, 1) 
+};
+#endif
+
index 7594deb..9e14524 100644 (file)
@@ -70,19 +70,19 @@ AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const char *name, const Cha
   ,fFilter(-1)
   ,fChi2OverNDFCut(3.0)
   ,fMaxDCA(3.0)
-  ,fMaxOpeningTheta(0.02)
-  ,fMaxOpeningPhi(0.1)
+//  ,fMaxOpeningTheta(0.02)
+//  ,fMaxOpeningPhi(0.1)
   ,fMaxOpening3D(TMath::Pi())
   ,fMaxInvMass(0.6)
   ,fSetMassConstraint(kFALSE)
   ,fArraytrack(NULL)
   ,fCounterPoolBackground(0)
   ,fListOutput(NULL)
-  ,fMCSourceee(NULL)
-  ,fUSignee(NULL)
-  ,fLSignee(NULL)
-  ,fUSignAngle(NULL)
-  ,fLSignAngle(NULL)
+  ,fMCSource(NULL)
+  ,fUSign(NULL)
+  ,fLSign(NULL)
+//  ,fUSignAngle(NULL)
+//  ,fLSignAngle(NULL)
 {
   //
   // Constructor
@@ -94,7 +94,7 @@ AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const char *name, const Cha
 //________________________________________________________________________
 AliHFENonPhotonicElectron::AliHFENonPhotonicElectron()
   :TNamed()
-  ,fMCEvent(0)
+  ,fMCEvent(NULL)
   ,fAODArrayMCInfo(NULL)
   ,fHFEBackgroundCuts(NULL)
   ,fPIDBackground(0x0)
@@ -104,19 +104,19 @@ AliHFENonPhotonicElectron::AliHFENonPhotonicElectron()
   ,fFilter(-1)
   ,fChi2OverNDFCut(3.0)
   ,fMaxDCA(3.0)
-  ,fMaxOpeningTheta(0.02)
-  ,fMaxOpeningPhi(0.1)
-  ,fMaxOpening3D(TMath::TwoPi())
+//  ,fMaxOpeningTheta(0.02)
+//  ,fMaxOpeningPhi(0.1)
+  ,fMaxOpening3D(TMath::Pi())
   ,fMaxInvMass(0.6)
   ,fSetMassConstraint(kFALSE)
   ,fArraytrack(NULL)
   ,fCounterPoolBackground(0)
   ,fListOutput(NULL)
-  ,fMCSourceee(NULL)
-  ,fUSignee(NULL)
-  ,fLSignee(NULL)
-  ,fUSignAngle(NULL)
-  ,fLSignAngle(NULL)
+  ,fMCSource(NULL)
+  ,fUSign(NULL)
+  ,fLSign(NULL)
+//  ,fUSignAngle(NULL)
+//  ,fLSignAngle(NULL)
 {
   //
   // Constructor
@@ -138,19 +138,19 @@ AliHFENonPhotonicElectron::AliHFENonPhotonicElectron(const AliHFENonPhotonicElec
   ,fFilter(ref.fFilter)
   ,fChi2OverNDFCut(ref.fChi2OverNDFCut)
   ,fMaxDCA(ref.fMaxDCA)
-  ,fMaxOpeningTheta(ref.fMaxOpeningTheta)
-  ,fMaxOpeningPhi(ref.fMaxOpeningPhi)
+//  ,fMaxOpeningTheta(ref.fMaxOpeningTheta)
+//  ,fMaxOpeningPhi(ref.fMaxOpeningPhi)
   ,fMaxOpening3D(ref.fMaxOpening3D)
   ,fMaxInvMass(ref.fMaxInvMass)
   ,fSetMassConstraint(ref.fSetMassConstraint)
   ,fArraytrack(NULL)
   ,fCounterPoolBackground(0)
   ,fListOutput(ref.fListOutput)
-  ,fMCSourceee(ref.fMCSourceee)
-  ,fUSignee(ref.fUSignee)
-  ,fLSignee(ref.fLSignee)
-  ,fUSignAngle(ref.fUSignAngle)
-  ,fLSignAngle(ref.fLSignAngle)
+  ,fMCSource(ref.fMCSource)
+  ,fUSign(ref.fUSign)
+  ,fLSign(ref.fLSign)
+//  ,fUSignAngle(ref.fUSignAngle)
+//  ,fLSignAngle(ref.fLSignAngle)
 {
   //
   // Copy Constructor
@@ -243,54 +243,40 @@ void AliHFENonPhotonicElectron::Init()
 
   const Int_t nDimMCSource=3;
   Int_t nBinMCSource[nDimMCSource] = {nBinsC,nBinsPt,nBinsSource};
-  fMCSourceee = new THnSparseF("fMCSourceee","fMCSourceee",nDimMCSource,nBinMCSource);
-  fMCSourceee->SetBinEdges(0,binLimC);
-  fMCSourceee->SetBinEdges(1,binLimPt);
-  fMCSourceee->SetBinEdges(2,binLimSource);
-/*  fMCSourceee->GetAxis(0)->SetTitle("Multiplicity (arb.unit)");
-  fMCSourceee->GetAxis(1)->SetTitle("#it{p}_{mc,T} (GeV/c)");
-  fMCSourceee->GetAxis(2)->SetTitle("MC");*/
-  fMCSourceee->Sumw2();
-  AliDebug(2,"AliHFENonPhotonicElectron: fMCSourceee");
+  fMCSource = new THnSparseF("fMCSource","fMCSource",nDimMCSource,nBinMCSource);
+  fMCSource->SetBinEdges(0,binLimC);
+  fMCSource->SetBinEdges(1,binLimPt);
+  fMCSource->SetBinEdges(2,binLimSource);
+  fMCSource->Sumw2();
+  AliDebug(2,"AliHFENonPhotonicElectron: fMCSource");
 
   // ee invariant mass Unlike Sign
   const Int_t nDimUSign=6;
   Int_t nBinUSign[nDimUSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource,nBinsAngle};
-  fUSignee = new THnSparseF("fUSignee","fUSignee",nDimUSign,nBinUSign);
-  fUSignee->SetBinEdges(0,binLimPhi);
-  fUSignee->SetBinEdges(1,binLimC);
-  fUSignee->SetBinEdges(2,binLimPt);
-  fUSignee->SetBinEdges(3,binLimInvMass);
-  fUSignee->SetBinEdges(4,binLimSource);
-  fUSignee->SetBinEdges(5,binLimAngle);
-/*  fUSignee->GetAxis(0)->SetTitle("#Delta #phi");
-  fUSignee->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fUSignee->GetAxis(2)->SetTitle("#it{p}_{T} (GeV/c)");
-  fUSignee->GetAxis(3)->SetTitle("m_{e^{#pm}e^{#mp}}} (GeV/#it{c}^{2}");
-  fUSignee->GetAxis(4)->SetTitle("MC");
-  fUSignee->GetAxis(5)->SetTitle("#alpha");*/
-  fUSignee->Sumw2();
-  AliDebug(2,"AliHFENonPhotonicElectron: fUSignee");
+  fUSign = new THnSparseF("fUSign","fUSign",nDimUSign,nBinUSign);
+  fUSign->SetBinEdges(0,binLimPhi);
+  fUSign->SetBinEdges(1,binLimC);
+  fUSign->SetBinEdges(2,binLimPt);
+  fUSign->SetBinEdges(3,binLimInvMass);
+  fUSign->SetBinEdges(4,binLimSource);
+  fUSign->SetBinEdges(5,binLimAngle);
+  fUSign->Sumw2();
+  AliDebug(2,"AliHFENonPhotonicElectron: fUSign");
 
   // ee invariant mass Like Sign
   const Int_t nDimLSign=6;
   Int_t nBinLSign[nDimLSign] = {nBinsPhi,nBinsC,nBinsPt,nBinsInvMass,nBinsSource,nBinsAngle};
-  fLSignee = new THnSparseF("fLSignee","fLSignee",nDimLSign,nBinLSign);
-  fLSignee->SetBinEdges(0,binLimPhi);
-  fLSignee->SetBinEdges(1,binLimC);
-  fLSignee->SetBinEdges(2,binLimPt);
-  fLSignee->SetBinEdges(3,binLimInvMass);
-  fLSignee->SetBinEdges(4,binLimSource);
-  fLSignee->SetBinEdges(5,binLimAngle);
-/*  fLSignee->GetAxis(0)->SetTitle("#Delta #phi");
-  fLSignee->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fLSignee->GetAxis(2)->SetTitle("#it{p}_{T} (GeV/c)");
-  fLSignee->GetAxis(3)->SetTitle("m_{e^{#pm}e^{#pm}}} (GeV/#it{c}^{2}");
-  fLSignee->GetAxis(4)->SetTitle("MC");
-  fLSignee->GetAxis(5)->SetTitle("#alpha");*/
-  fLSignee->Sumw2();
-  AliDebug(2,"AliHFENonPhotonicElectron: fLSignee");
-
+  fLSign = new THnSparseF("fLSign","fLSign",nDimLSign,nBinLSign);
+  fLSign->SetBinEdges(0,binLimPhi);
+  fLSign->SetBinEdges(1,binLimC);
+  fLSign->SetBinEdges(2,binLimPt);
+  fLSign->SetBinEdges(3,binLimInvMass);
+  fLSign->SetBinEdges(4,binLimSource);
+  fLSign->SetBinEdges(5,binLimAngle);
+  fLSign->Sumw2();
+  AliDebug(2,"AliHFENonPhotonicElectron: fLSign");
+
+/*
   // ee angle Unlike Sign
   const Int_t nDimUSignAngle=3;
   Int_t nBinUSignAngle[nDimUSignAngle] = {nBinsAngle,nBinsC,nBinsSource};
@@ -298,9 +284,6 @@ void AliHFENonPhotonicElectron::Init()
   fUSignAngle->SetBinEdges(0,binLimAngle);
   fUSignAngle->SetBinEdges(1,binLimC);
   fUSignAngle->SetBinEdges(2,binLimSource);
-/*  fUSignAngle->GetAxis(0)->SetTitle("#alpha");
-  fUSignAngle->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fUSignAngle->GetAxis(2)->SetTitle("MC");*/
   fUSignAngle->Sumw2();
   AliDebug(2,"AliHFENonPhotonicElectron: fUSignAngle");
 
@@ -311,18 +294,15 @@ void AliHFENonPhotonicElectron::Init()
   fLSignAngle->SetBinEdges(0,binLimAngle);
   fLSignAngle->SetBinEdges(1,binLimC);
   fLSignAngle->SetBinEdges(2,binLimSource);
-/*  fLSignAngle->GetAxis(0)->SetTitle("#alpha");
-  fLSignAngle->GetAxis(1)->SetTitle("Multiplicity (arb.unit)");
-  fLSignAngle->GetAxis(2)->SetTitle("MC");*/
   fLSignAngle->Sumw2();
   AliDebug(2,"AliHFENonPhotonicElectron: fLSignAngle");
+*/
 
-
-  fListOutput->Add(fMCSourceee);
-  fListOutput->Add(fUSignee);
-  fListOutput->Add(fLSignee);
-  fListOutput->Add(fUSignAngle);
-  fListOutput->Add(fLSignAngle);
+  fListOutput->Add(fMCSource);
+  fListOutput->Add(fUSign);
+  fListOutput->Add(fLSign);
+//  fListOutput->Add(fUSignAngle);
+//  fListOutput->Add(fLSignAngle);
 
 }
 
@@ -385,7 +365,10 @@ Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent,
 
     if(aodeventu)
     {
-      // AOD
+      /**                              **
+       *       AOD Analysis             *
+       **                              **/
+
       AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
       if(aodtrack)
       {
@@ -421,7 +404,10 @@ Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent,
     }
     else
     {
-      // ESD
+      /**                              **
+       *       ESD Analysis             *
+       **                              **/
+
       AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
       if(esdtrack)
       {
@@ -435,8 +421,8 @@ Int_t AliHFENonPhotonicElectron::FillPoolAssociatedTracks(AliVEvent *inputEvent,
       // PID track cuts
       AliHFEpidObject hfetrack2;
 
-      if(!aodeventu) hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
-      else hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
+      if(!aodeventu)   hfetrack2.SetAnalysisType(AliHFEpidObject::kESDanalysis);
+      else             hfetrack2.SetAnalysisType(AliHFEpidObject::kAODanalysis);
 
       hfetrack2.SetRecTrack(track);
       if(binct>-1)
@@ -494,12 +480,12 @@ Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1,
   AliDebug(2,Form("process track %d",iTrack1));
 
   //Set Fill-Arrays for THnSparse
-  Double_t value_MCSource[3] = { binct, track1->Pt(), source};         //Centrality    Pt              Source
-  Double_t valueee[6] = { deltaphi, binct, track1->Pt(), -1, source,-1};       //DeltaPhi      Centrality      Pt      InvariantMass   Source
-  Double_t valueAngle[3] = {   -1,     binct,          source};                //Angle         Centrality      Source
+  Double_t valueMCSource[3]    = { -1, -1, source};                                            //Centrality    Pt              Source
+  Double_t valueSign[6]                = { deltaphi, binct, track1->Pt(), -1, source, -1};             //DeltaPhi      Centrality      Pt      InvariantMass   Source
+//  Double_t valueAngle[3]     = { -1, binct, source};                                         //Angle         Centrality      Source
 
-  Bool_t USignPhotonic = kFALSE;
-  Bool_t LSignPhotonic = kFALSE;
+  Bool_t uSignPhotonic = kFALSE;
+  Bool_t lSignPhotonic = kFALSE;
   Bool_t hasdcaT1 = kFALSE;
   Bool_t hasdcaT2 = kFALSE;
 
@@ -542,8 +528,10 @@ Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1,
   AliKFParticle *ktrack2;
   AliKFVertex primV(*(vEvent->GetPrimaryVertex()));
 
-  // FILL MCsource TODO: correct?
-  fMCSourceee->Fill(&value_MCSource[0],weight);
+  //! FILL MCsource TODO: Use MC information!!!
+  /**  if(fAODArrayMCInfo)     valueMCSource[3] = {fAODArrayMCInfo};
+       if(fMCEvent)            valueMCSource[3] = {fMCEvent};          */
+  fMCSource->Fill(&valueMCSource[0],weight);
 
 
   for(Int_t idex = 0; idex < fCounterPoolBackground; idex++)
@@ -565,8 +553,8 @@ Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1,
     fCharge2 = track2->Charge();               //Charge from track2
 
     // Reset the MC info
-    valueAngle[2] = source;
-    valueee[4] = source;
+    //valueAngle[2] = source;
+    valueSign[4] = source;
 
     // track cuts and PID already done
 
@@ -582,21 +570,21 @@ Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1,
        {
          if(source == kElectronfromconversion)
          {
-           valueAngle[2] = kElectronfromconversionboth;
-           valueee[4] = kElectronfromconversionboth;
+           //valueAngle[2] = kElectronfromconversionboth;
+           valueSign[4] = kElectronfromconversionboth;
            numberfound++;
          }
 
          if(source == kElectronfrompi0)
          {
-           valueAngle[2] = kElectronfrompi0both;
-           valueee[4] = kElectronfrompi0both;
+           //valueAngle[2] = kElectronfrompi0both;
+           valueSign[4] = kElectronfrompi0both;
          }
 
          if(source == kElectronfrometa)
          {
-           valueAngle[2] = kElectronfrometaboth;
-           valueee[4] = kElectronfrometaboth;
+           //valueAngle[2] = kElectronfrometaboth;
+           valueSign[4] = kElectronfrometaboth;
          }
        }
       }
@@ -628,27 +616,27 @@ Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1,
       invmassESD  = mother.M();
       angleESD    = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
 
-      valueAngle[0] = angleESD;
-      valueee[5] = angleESD;
+      //valueAngle[0] = angleESD;
+      valueSign[5] = angleESD;
 
-      if((fCharge1*fCharge2)>0.0)      fLSignAngle->Fill(&valueAngle[0],weight);
-      else                             fUSignAngle->Fill(&valueAngle[0],weight);
+      //if((fCharge1*fCharge2)>0.0)    fLSignAngle->Fill(&valueAngle[0],weight);
+      //else                           fUSignAngle->Fill(&valueAngle[0],weight);
 
       if(angleESD > fMaxOpening3D) continue;                            //! Cut on Opening Angle
-      valueee[3] = invmassESD;
+      valueSign[3] = invmassESD;
 
-      if((fCharge1*fCharge2)>0.0)      fLSignee->Fill(&valueee[0],weight);
-      else                             fUSignee->Fill(&valueee[0],weight);
+      if((fCharge1*fCharge2)>0.0)      fLSign->Fill(&valueSign[0],weight);
+      else                             fUSign->Fill(&valueSign[0],weight);
 
       if(invmassESD > fMaxInvMass) continue;                           //! Cut on Invariant Mass
 
-      if((fCharge1*fCharge2)>0.0)      LSignPhotonic=kTRUE;
-      else                             USignPhotonic=kTRUE;
+      if((fCharge1*fCharge2)>0.0)      lSignPhotonic=kTRUE;
+      else                             uSignPhotonic=kTRUE;
     }
     else
     {
       /**                              *
-       *       AliKF-Analysis          *
+       *       AOD-AliKF-Analysis      *
        **                              */
 
       fPDGtrack1 = 11;
@@ -683,29 +671,29 @@ Int_t AliHFENonPhotonicElectron::LookAtNonHFE(Int_t iTrack1, AliVTrack *track1,
       recoGamma.GetMass(invmassAOD,width);
       angleAOD = ktrack1->GetAngle(*ktrack2);
 
-      valueAngle[0] = angleAOD;
-      valueee[5] = angleAOD;
+      //valueAngle[0] = angleAOD;
+      valueSign[5] = angleAOD;
 
-      if((fCharge1*fCharge2)>0.0)      fLSignAngle->Fill(&valueAngle[0],weight);
-      else                             fUSignAngle->Fill(&valueAngle[0],weight);
+      //if((fCharge1*fCharge2)>0.0)    fLSignAngle->Fill(&valueAngle[0],weight);
+      //else                           fUSignAngle->Fill(&valueAngle[0],weight);
 
       if(angleAOD > fMaxOpening3D) continue;                           //! Cut on Opening Angle
 
-      valueee[3] = invmassAOD;
+      valueSign[3] = invmassAOD;
 
-      if((fCharge1*fCharge2)>0.0)      fLSignee->Fill(&valueee[0],weight);
-      else                             fUSignee->Fill(&valueee[0],weight);
+      if((fCharge1*fCharge2)>0.0)      fLSign->Fill(&valueSign[0],weight);
+      else                             fUSign->Fill(&valueSign[0],weight);
 
       if(invmassAOD > fMaxInvMass) continue;                           //! Cut on Invariant Mass
 
-      if((fCharge1*fCharge2)>0.0)      LSignPhotonic=kTRUE;
-      else                             USignPhotonic=kTRUE;
+      if((fCharge1*fCharge2)>0.0)      lSignPhotonic=kTRUE;
+      else                             uSignPhotonic=kTRUE;
     }
   }
 
-  if( USignPhotonic &&  LSignPhotonic) taggedphotonic = 6;
-  if(!USignPhotonic &&  LSignPhotonic) taggedphotonic = 4;
-  if( USignPhotonic && !LSignPhotonic) taggedphotonic = 2;
+  if( uSignPhotonic &&  lSignPhotonic) taggedphotonic = 6;
+  if(!uSignPhotonic &&  lSignPhotonic) taggedphotonic = 4;
+  if( uSignPhotonic && !lSignPhotonic) taggedphotonic = 2;
 
   return taggedphotonic;
 }
index 641d19e..6858619 100644 (file)
@@ -66,17 +66,17 @@ class AliHFENonPhotonicElectron : public TNamed {
 
   void  SetMaxInvMass          (Double_t MaxInvMass)           { fMaxInvMass           = MaxInvMass; };
   void  SetMaxOpening3D                (Double_t MaxOpening3D)         { fMaxOpening3D         = MaxOpening3D; };
-  void  SetMaxOpeningTheta     (Double_t MaxOpeningTheta)      { fMaxOpeningTheta      = MaxOpeningTheta; };
-  void  SetMaxOpeningPhi       (Double_t MaxOpeningPhi)        { fMaxOpeningPhi        = MaxOpeningPhi; };
+//  void  SetMaxOpeningTheta   (Double_t MaxOpeningTheta)      { fMaxOpeningTheta      = MaxOpeningTheta; };
+//  void  SetMaxOpeningPhi     (Double_t MaxOpeningPhi)        { fMaxOpeningPhi        = MaxOpeningPhi; };
   void  SetAlgorithmMA         (Bool_t algorithmMA)            { fAlgorithmMA          = algorithmMA; };
   void  SetMassConstraint      (Bool_t MassConstraint)         { fSetMassConstraint    = MassConstraint; };
 
   TList      *GetListOutput()       const { return fListOutput; };
-  THnSparseF *GetMCSourceeeHisto() const { return fMCSourceee; };
-  THnSparseF *GetUSigneeHisto()    const { return fUSignee; };
-  THnSparseF *GetLSigneeHisto()    const { return fLSignee; };
-  THnSparseF *GetUSignAngleHisto() const { return fUSignAngle; };
-  THnSparseF *GetLSignAngleHisto() const { return fLSignAngle; };
+  THnSparseF *GetMCSourceHisto() const { return fMCSource; };
+  THnSparseF *GetUSignHisto()    const { return fUSign; };
+  THnSparseF *GetLSignHisto()    const { return fLSign; };
+//  THnSparseF *GetUSignAngleHisto() const { return fUSignAngle; };
+//  THnSparseF *GetLSignAngleHisto() const { return fLSignAngle; };
 
   void     Init                                ();
   void     InitRun                     (const AliVEvent *inputEvent,const AliPIDResponse *pidResponse);
@@ -93,7 +93,7 @@ class AliHFENonPhotonicElectron : public TNamed {
 
 
  private:
-  AliMCEvent           *fMCEvent;                      //! MC event esd
+  AliMCEvent           *fMCEvent;                      //! MC event ESD
   TClonesArray         *fAODArrayMCInfo;               //! MC info particle AOD
   AliESDtrackCuts      *fHFEBackgroundCuts;            // HFE background cuts
   AliHFEpid            *fPIDBackground;                // PID background cuts
@@ -103,22 +103,22 @@ class AliHFENonPhotonicElectron : public TNamed {
   UInt_t                fFilter;                       // filter AOD status
   Double_t              fChi2OverNDFCut;               // Limit chi2
   Double_t              fMaxDCA;                       // Limit dca
-  Double_t              fMaxOpeningTheta;              // Limit opening angle in theta
-  Double_t              fMaxOpeningPhi;                // Limit opening angle in phi
+//  Double_t            fMaxOpeningTheta;              // Limit opening angle in theta
+//  Double_t            fMaxOpeningPhi;                // Limit opening angle in phi
   Double_t              fMaxOpening3D;                 // Limit opening 3D
   Double_t              fMaxInvMass;                   // Limit invariant mass
   Bool_t                fSetMassConstraint;            // Set mass constraint
   TArrayI              *fArraytrack;                   //! list of tracks
   Int_t                         fCounterPoolBackground;        // number of tracks
   TList                        *fListOutput;                   // List of histos
-  THnSparseF           *fMCSourceee;                   //! centrality, pt, Source MC
-  THnSparseF           *fUSignee;                      //! delta phi, c, pt, inv, source
-  THnSparseF           *fLSignee;                      //! delta phi, c, pt, inv, source
-  THnSparseF           *fUSignAngle;                   //! angle, c, source
-  THnSparseF           *fLSignAngle;                   //! angle, c, source
+  THnSparseF           *fMCSource;                     //! centrality, pt, Source MC
+  THnSparseF           *fUSign;                        //! delta phi, c, pt, inv, source
+  THnSparseF           *fLSign;                        //! delta phi, c, pt, inv, source
+//  THnSparseF         *fUSignAngle;                   //! angle, c, source
+//  THnSparseF         *fLSignAngle;                   //! angle, c, source
 
 
-  AliHFENonPhotonicElectron(const AliHFENonPhotonicElectron&); // not implemented
+  AliHFENonPhotonicElectron(const AliHFENonPhotonicElectron&ref); 
 
   ClassDef(AliHFENonPhotonicElectron, 1); //!example of analysis
 };
index f474474..ebf85bf 100644 (file)
@@ -194,6 +194,7 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){
                 << "Centrality="              << centrality
                 << "VertexZ="                 << vtx[2]
                 << "NumberOfContributors="    << ncontrib
+                << "run="                     << run
                 << "\n";
 
   // Common variables
index a96b687..33f53de 100644 (file)
 #include <TBits.h>
 #include <TString.h>
 #include <TArrayI.h>
-#include <TTree.h>
 
 #include "AliAnalysisManager.h"
 #include "AliCentrality.h"
 #include "AliAODTrack.h"
 #include "AliAODEvent.h"
+#include "AliAODPid.h"
 #include "AliHFEcuts.h"
 #include "AliHFEextraCuts.h"
 #include "AliInputEventHandler.h"
@@ -60,35 +60,7 @@ AliHFEdebugTreeTaskAOD::AliHFEdebugTreeTaskAOD():
   fNclustersTPCPID(0),
   fNclustersITS(2),
   fFilename("HFEtree.root"),
-  fDebugstream(kFALSE),
-  fDebugTree(NULL),
-  fDebugTreee(NULL),
-  fCentrality(-1.),
-  fRun(-1),
-  fDoublec(0),
-  fMomentum(0.),
-  fMomentumTPC(0.),
-  fTransverseMomentum(0.),
-  fEta(0.),
-  fPhi(0.),
-  fCharge(-1),
-  fNClustersTPCall(0),
-  fNClustersTPCPID(0),
-  fNClustersTPCshared(0),
-  fNCrossedRowsTPC(0),
-  fClusterRatioTPCall(0.),
-  fNClustersITS(0),
-  fStatusL0(-1),
-  fStatusL1(-1),
-  fSigmaTOF(0.),
-  fSigmaTPC(0.),
-  fDcaxy(0.),
-  fDcaz(0.),
-  fFilter2(0),
-  fFilter4(0),
-  fSource(-1),
-  fEr(0.0),
-  fSignal(0.)
+  fDebugTree(NULL)
 {
 
 }
@@ -106,45 +78,15 @@ AliHFEdebugTreeTaskAOD::AliHFEdebugTreeTaskAOD(const char *name):
   fNclustersTPCPID(0),
   fNclustersITS(2),
   fFilename("HFEtree.root"),
-  fDebugstream(kFALSE),
-  fDebugTree(NULL),
-  fDebugTreee(NULL),
-  fCentrality(-1.),
-  fRun(-1),
-  fDoublec(0),
-  fMomentum(0.),
-  fMomentumTPC(0.),
-  fTransverseMomentum(0.),
-  fEta(0.),
-  fPhi(0.),
-  fCharge(-1),
-  fNClustersTPCall(0),
-  fNClustersTPCPID(0),
-  fNClustersTPCshared(0),
-  fNCrossedRowsTPC(0),
-  fClusterRatioTPCall(0.),
-  fNClustersITS(0),
-  fStatusL0(-1),
-  fStatusL1(-1),
-  fSigmaTOF(0.),
-  fSigmaTPC(0.),
-  fDcaxy(0.),
-  fDcaz(0.),
-  fFilter2(0),
-  fFilter4(0),
-  fSource(-1),
-  fEr(0.0),
-  fSignal(0.)
+  fDebugTree(NULL)
 {
   fTPCpid = new AliHFEpidTPC("QAtpcPID");
-  DefineOutput(1, TTree::Class());
 }
 
 AliHFEdebugTreeTaskAOD::~AliHFEdebugTreeTaskAOD(){
 
     if(fDebugTree) delete fDebugTree;
     if(fTPCpid) delete fTPCpid;
-    //if(fDebugTreee) delete fDebugTreee;
 }
 
 void AliHFEdebugTreeTaskAOD::UserCreateOutputObjects(){
@@ -153,40 +95,7 @@ void AliHFEdebugTreeTaskAOD::UserCreateOutputObjects(){
   //
 
   //printf("test\n");
-  if(fDebugstream) {
-    fDebugTree = new TTreeSRedirector(fFilename.Data());
-  }
-  else {
-    // other possibility
-    fDebugTreee =  new TTree("PIDdebug","PIDdebug");
-    fDebugTreee->Branch("centrality",&fCentrality);
-    fDebugTreee->Branch("run",&fRun);
-    fDebugTreee->Branch("doublec",&fDoublec);
-    fDebugTreee->Branch("p",&fMomentum);
-    fDebugTreee->Branch("ptpc",&fMomentumTPC);
-    fDebugTreee->Branch("pt",&fTransverseMomentum);
-    fDebugTreee->Branch("eta",&fEta);
-    fDebugTreee->Branch("phi",&fPhi);
-    fDebugTreee->Branch("charge",&fCharge);
-    fDebugTreee->Branch("nclustersTPCall",&fNClustersTPCall);
-    fDebugTreee->Branch("nclustersTPCPID",&fNClustersTPCPID);
-    fDebugTreee->Branch("nclustersTPCshared",&fNClustersTPCshared);
-    fDebugTreee->Branch("nCrossedRowsTPC",&fNCrossedRowsTPC);
-    fDebugTreee->Branch("clusterRatioTPCall",&fClusterRatioTPCall);
-    fDebugTreee->Branch("nclustersITS",&fNClustersITS);
-    fDebugTreee->Branch("statusL0",&fStatusL0);
-    fDebugTreee->Branch("statusL1",&fStatusL1);
-    fDebugTreee->Branch("sigmaTOF",&fSigmaTOF);
-    fDebugTreee->Branch("sigmaTPC",&fSigmaTPC);
-    fDebugTreee->Branch("dcaxy",&fDcaxy);
-    fDebugTreee->Branch("dcaz",&fDcaz);
-    fDebugTreee->Branch("filter2",&fFilter2);
-    fDebugTreee->Branch("filter4",&fFilter4);
-    fDebugTreee->Branch("source",&fSource);
-    fDebugTreee->Branch("eR",&fEr);
-    fDebugTreee->Branch("signal",&fSignal);
-    PostData(1,fDebugTreee);
-  }
+  fDebugTree = new TTreeSRedirector(fFilename.Data());
 
  // printf("testa\n");
   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
@@ -236,6 +145,8 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     return;
   }
 
+  AliAODTrack copyTrack;
+
   // MC info
   Bool_t mcthere = kTRUE;
   AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
@@ -271,9 +182,12 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
   }
   fExtraCuts->SetRecEventInfo(fInputEvent);
 
+
   // Get run number
-  fRun = fInputEvent->GetRunNumber();
-  
+  Int_t run = fInputEvent->GetRunNumber();
+
   // Derive trigger 
   UInt_t trigger = fInputHandler->IsEventSelected();
   Bool_t isMBTrigger = trigger & AliVEvent::kMB;
@@ -288,10 +202,19 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
   Double_t ncontrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
 
   // Get centrality
-  fCentrality = -1.;
+  Float_t centrality = -1.;
   AliCentrality *hicent = fInputEvent->GetCentrality();
-  fCentrality = hicent->GetCentralityPercentile("V0M");
-  
+  centrality = hicent->GetCentralityPercentile("V0M");
+
+  // Store event selection variables
+  (*fDebugTree) << "EventDebug"
+                << "Centrality="              << centrality
+                << "VertexZ="                 << vtx[2]
+                << "NumberOfContributors="    << ncontrib
+                << "run="                 << run
+                << "\n";
+
+
 
   // Look for kink mother
   Int_t numberofvertices = aodE->GetNumberOfVertices();
@@ -313,38 +236,112 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     }
   }
   //printf("Number of kink mother in the events %d\n",numberofmotherkink);
+
+  //
+  // Loop on MC tracks only
+  //
+  AliAODMCParticle *mctrack = NULL;
+  // Monte-Carlo info
+  Double_t eR,vx,vy,vz;
+  Double_t chargemc, etamc, phimc, momentummc, transversemomentummc;
+  Int_t source,pdg,signal;
+  if(mcthere) {
+    for(Int_t itrack = 0; itrack < fAODArrayMCInfo->GetEntriesFast(); itrack++) {
+      mctrack = (AliAODMCParticle *)(fAODArrayMCInfo->At(itrack));
+      if(!mctrack) continue;
+      signal = 0;
+      if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) signal = 1;
+      // Kinematics
+      chargemc = mctrack->Charge() > 0. ? 1. : -1.;
+      momentummc = mctrack->P() * chargemc;
+      transversemomentummc = mctrack->Pt() * chargemc;
+      etamc = mctrack->Eta();
+      phimc = mctrack->Phi();
+      pdg = mctrack->GetPdgCode();
+      
+      // Get Production Vertex in radial direction
+      vx = mctrack->Xv();
+      vy = mctrack->Yv(); 
+      vz = mctrack->Zv(); 
+      eR = TMath::Sqrt(vx*vx+vy*vy);
+
+      // Get Mother PDG code of the particle
+      Int_t motherPdg = 0;
+      Int_t motherlabel = TMath::Abs(mctrack->GetMother());
+      if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
+        AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
+        if(mother) motherPdg = mother->GetPdgCode();
+      }
+      
+      // derive source
+      source = 5;
+      if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
+      else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
+      else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
+      else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
+      else if(TMath::Abs(pdg) == 11) source = 4;
+      else source = 5;
+
+      (*fDebugTree) << "MCDebug"
+                    << "centrality="          << centrality
+                    << "MBtrigger="           << isMBTrigger 
+                    << "CentralTrigger="      << isCentralTrigger
+                    << "SemicentralTrigger="  << isSemicentralTrigger
+                    << "EMCALtrigger="        << isEMCALTrigger
+                    << "run="                 << run
+                    << "p="                   << momentummc
+                    << "pt="                  << transversemomentummc
+                    << "eta="                 << etamc
+                    << "phi="                 << phimc
+                    << "pdg="                 << pdg
+                   << "px="                  << vx
+                   << "py="                  << vy
+                   << "pz="                  << vz
+                    << "ProductionVertex="    << eR
+                    << "motherPdg="           << motherPdg
+                    << "source="              << source
+                    << "\n";
+    }
+  }
   
   // Common variables
-  //Double_t charge, eta, phi, momentum, momentumTPC, transversemomentum;
-
+  Double_t charge, eta, phi, momentum, momentumTPC, transversemomentum;
+  
   //
   // Loop on reconstructed tracks
   //
 
   TArrayI *arraytrack = new TArrayI(fInputEvent->GetNumberOfTracks());
+  Int_t counterdc=0;
   
   AliAODTrack *track = 0x0;
-  AliAODMCParticle *mctrack = NULL;
   for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
     // fill the tree
     track = dynamic_cast<AliAODTrack *>(fInputEvent->GetTrack(itrack));
     if(!track) continue;
     // Cut track (Only basic track cuts)
-//    printf("testv\n");
+    // printf("testv\n");
     if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
     //
     //printf("testu\n");
-    fSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
-    fSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
+    Double_t nSigmaTOF = pid->NumberOfSigmasTOF(track, AliPID::kElectron);
+    Double_t nSigmaTPC = pid->NumberOfSigmasTPC(track, AliPID::kElectron);
     Double_t tPCdEdx = track->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.;
+    // Eta correction
+    copyTrack.~AliAODTrack();
+    new(&copyTrack) AliAODTrack(*track);
+    if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kAODanalysis);
+    if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kAODanalysis);
+    Double_t nSigmaTPCcorr = pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron);
+    Double_t tPCdEdxcorr = copyTrack.GetDetPid() ? copyTrack.GetDetPid()->GetTPCsignal() : 0.;
    
     // Kinematics
-    fCharge = track->Charge() > 0 ? 1. : -1.;
-    fEta = track->Eta();
-    fPhi = track->Phi();
-    fMomentum = track->P() * fCharge;
-    fTransverseMomentum = track->Pt() * fCharge;
-    fMomentumTPC = track->GetDetPid() ? track->GetDetPid()->GetTPCmomentum() : track->P();
+    charge = track->Charge() > 0 ? 1. : -1.;
+    eta = track->Eta();
+    phi = track->Phi();
+    momentum = track->P() * charge;
+    transversemomentum = track->Pt() * charge;
+    momentumTPC = track->GetDetPid() ? track->GetDetPid()->GetTPCmomentum() : track->P();
   
     // status
     ULong_t status = track->GetStatus();
@@ -352,40 +349,63 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) itsrefit = 1;
     Int_t tpcrefit=0;
     if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) tpcrefit = 1;
+    Int_t tofpid=0;
+    if((status & AliESDtrack::kTOFpid) == AliESDtrack::kTOFpid) tofpid = 1;
 
     // ITS number of clusters
-    fNClustersITS = (Int_t) track->GetITSNcls();
+    UChar_t nclustersITS = track->GetITSNcls();
     // TPC number of clusters (different definitions)
     UChar_t nclustersTPCfit = track->GetTPCNcls();
-    //UChar_t nclustersTPCall = 0;
+    UChar_t nclustersTPCall = 0;
     const TBits &clusterTPC = track->GetTPCClusterMap();
-    fNClustersTPCall = (Int_t) clusterTPC.CountBits();
-    fNClustersTPCPID = (Int_t) track->GetTPCsignalN();
+    nclustersTPCall = clusterTPC.CountBits();
+    UChar_t nclustersTPCPID = track->GetTPCsignalN();
     UChar_t nfindableTPC =  track->GetTPCNclsF();
     Double_t clusterRatioTPCfit = 0.0;
     if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCfit = static_cast<Double_t>(nclustersTPCfit)/static_cast<Double_t>(nfindableTPC);
-    //Double_t clusterRatioTPCall = 0.0;
-    if((static_cast<Double_t>(nfindableTPC))>0.0) fClusterRatioTPCall = static_cast<Float_t>(fNClustersTPCall)/static_cast<Float_t>(nfindableTPC);
-    fNClustersTPCshared = 0;
-    fNCrossedRowsTPC = (Int_t) track->GetTPCNCrossedRows();
+    Double_t clusterRatioTPCall = 0.0;
+    if((static_cast<Double_t>(nfindableTPC))>0.0) clusterRatioTPCall = static_cast<Double_t>(nclustersTPCall)/static_cast<Double_t>(nfindableTPC);
+    UChar_t nclustersTPCshared = 0;
+    Float_t ncrossedRowsTPC = track->GetTPCNCrossedRows();
     const TBits &sharedTPC = track->GetTPCSharedMap();
-    for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) fNClustersTPCshared++;
+    for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
     // TRD clusters and tracklets
     UChar_t nclustersTRD = track->GetTRDncls();
     UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID();
-    Int_t   nslicesTRD = track->GetNumberOfTRDslices();
+    Int_t   nslicesTRD = track->GetDetPid()->GetTRDnSlices();
     Int_t   chi2TRD = track->GetTRDchi2();
+    AliAODPid* aodpid= track->GetDetPid();
+    Double_t* arraytrdsignals;
+    arraytrdsignals=aodpid->GetTRDsignal();
+    Int_t nslicetemp=0;
+    Int_t trdlayer[6];
+    for(Int_t iplane = 0; iplane < 6; iplane++){
+       trdlayer[iplane]=0;
+    }
+
+    for(Int_t iplane = 0; iplane < 6; iplane++){
+       nslicetemp=0;
+       for(Int_t b=(iplane*8);b<((iplane*8)+8);b++)
+       {
+           if(ntrackletsTRDPID>0)
+           {
+               if(arraytrdsignals[b]>0.001) nslicetemp++;
+           }
+       }
+       if(nslicetemp > 0) trdlayer[iplane]=1;
+    }
+
     // ITS and TRD acceptance maps
     UChar_t itsPixel = track->GetITSClusterMap();
-    fStatusL0 = 0;
-    if(TESTBIT(itsPixel, 0)) fStatusL0 = 1; 
-    fStatusL1 = 0;
-    if(TESTBIT(itsPixel, 1)) fStatusL1 = 1; 
+    Bool_t statusL0 = kFALSE;
+    if(TESTBIT(itsPixel, 0)) statusL0 = kTRUE; 
+    Bool_t statusL1 = kFALSE;
+    if(TESTBIT(itsPixel, 1)) statusL1 = kTRUE; 
 
     // HFE DCA
-    fDcaxy = -999.;
-    fDcaz = -999.;
-    fExtraCuts->GetImpactParameters((AliVTrack *)track,fDcaxy,fDcaz);
+    Float_t dcaxy = -999.;
+    Float_t dcaz = -999.;
+    fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
 
     // Kink
     Int_t kink = 0;
@@ -404,15 +424,16 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     Int_t id = track->GetID();
 
     // Double counted
-    fDoublec = 0;
-    for(Int_t l=0; l < itrack; l++){
+    Int_t doublec = 0;
+    for(Int_t l=0; l < counterdc; l++){
       Int_t iTrack2 = arraytrack->At(l);
-      if(iTrack2==id) fDoublec=1;
+      if(iTrack2==id) doublec=1;
     }
     //printf("Doublec %d\n",doublec);
 
     // Add the id at this place
-    arraytrack->AddAt(id,itrack);
+    arraytrack->AddAt(id,counterdc);
+    counterdc++;
 
     // Filter
     Int_t filter[20];
@@ -425,9 +446,9 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
     }
     Int_t filter0 = filter[0];
     Int_t filter1 = filter[1];
-    fFilter2 = filter[2];
+    Int_t filter2 = filter[2];
     Int_t filter3 = filter[3];
-    fFilter4 = filter[4];
+    Int_t filter4 = filter[4];
     Int_t filter5 = filter[5];
     Int_t filter6 = filter[6];
     Int_t filter7 = filter[7];
@@ -448,19 +469,13 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
 
     //printf("track\n");
 
-    // Monte-Carlo info
-    Double_t vx,vy,vz;
-    fEr = 0.0;
-    Double_t chargemc, etamc, phimc, momentummc, transversemomentummc;
-    Int_t pdg;
-    fSignal = 0;
-    fSource = 0;
+    signal = 0;
     if(mcthere){
       Int_t label = TMath::Abs(track->GetLabel());
       if(label && label < fAODArrayMCInfo->GetEntriesFast())
         mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(label));
       if(!mctrack) continue;
-      if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) fSignal = 1;
+      if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) signal = 1;
       // Kinematics
       chargemc = mctrack->Charge() > 0. ? 1. : -1.;
       momentummc = mctrack->P() * chargemc;
@@ -473,7 +488,7 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
       vx = mctrack->Xv();
       vy = mctrack->Yv(); 
       vz = mctrack->Zv(); 
-      fEr = TMath::Sqrt(vx*vx+vy*vy);
+      eR = TMath::Sqrt(vx*vx+vy*vy);
       
       // Get Mother PDG code of the particle
       Int_t motherPdg = 0;
@@ -484,100 +499,105 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
       }
       
       // derive source
-      fSource = 5;
-      if(fSignalCuts->IsCharmElectron(mctrack)) fSource = 0;
-      else if(fSignalCuts->IsBeautyElectron(mctrack)) fSource = 1;
-      else if(fSignalCuts->IsGammaElectron(mctrack)) fSource = 2;
-      else if(fSignalCuts->IsNonHFElectron(mctrack)) fSource = 3;
-      else if(TMath::Abs(pdg) == 11) fSource = 4;
-      else fSource = 5;
+      source = 5;
+      if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
+      else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
+      else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
+      else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
+      else if(TMath::Abs(pdg) == 11) source = 4;
+      else source = 5;
       
     }
 
     
     // Fill Tree
     //printf("Fill\n");
-    if(fDebugstream) {
-      (*fDebugTree) << "PIDdebug"
-                   << "centrality="          << fCentrality
-                   << "MBtrigger="           << isMBTrigger 
-                   << "CentralTrigger="      << isCentralTrigger
-                   << "SemicentralTrigger="  << isSemicentralTrigger
-                   << "EMCALtrigger="        << isEMCALTrigger
-                   << "run="                 << fRun
-                   << "eventnb="             << eventnb
-                   << "vx="                  << vtx[0]
-                   << "vy="                  << vtx[1]
-                   << "vz="                  << vtx[2]
-                   << "ncontrib="            << ncontrib
-                   << "id="                  << id
-                   << "dc="                  << fDoublec
-                   << "p="                   << fMomentum
-                   << "ptpc="                << fMomentumTPC
-                   << "pt="                  << fTransverseMomentum
-                   << "eta="                 << fEta
-                   << "phi="                 << fPhi
-                   << "charge="              << fCharge
-                   << "itsrefit="            << itsrefit
-                   << "tpcrefit="            << tpcrefit
-                   << "nclustersTPC="        << nclustersTPCfit
-                   << "nclustersTPCall="     << fNClustersTPCall
-                   << "nclustersTPCPID="     << fNClustersTPCPID
-                   << "nclustersTPCshared="  << fNClustersTPCshared
-                   << "ncrossedRowsTPC="     << fNCrossedRowsTPC
-                   << "clusterRatioTPC="     << clusterRatioTPCfit
-                   << "clusterRatioTPCall="  << fClusterRatioTPCall
-                   << "nclustersITS="        << fNClustersITS
-                   << "nclustersTRD="        << nclustersTRD
-                   << "ntrackletsTRD="       << ntrackletsTRDPID
-                   << "nslicesTRD="          << nslicesTRD
-                   << "chi2TRD="             << chi2TRD
-                   << "statusITS0="          << fStatusL0
-                   << "statusITS1="          << fStatusL1
-                   << "TOFsigmaEl="          << fSigmaTOF
-                   << "TPCsigmaEl="          << fSigmaTPC
-                   << "TPCdEdx="             << tPCdEdx
-                   << "dcaR="                << fDcaxy
-                   << "dcaZ="                << fDcaz
-                   << "kinkdaughter="        << kink
-                   << "kinkmother="          << kinkmotherpass
-                   << "nbofmotherkink="      << numberofmotherkink
-                   << "filter0="             << filter0
-                   << "filter1="             << filter1
-                   << "filter2="             << fFilter2
-                   << "filter3="             << filter3
-                   << "filter4="             << fFilter4
-                   << "filter5="             << filter5
-                   << "filter6="             << filter6
-                   << "filter7="             << filter7
-                   << "filter8="             << filter8
-                   << "filter9="             << filter9
-                   << "filter10="            << filter10
-                   << "filter11="            << filter11
-                   << "filter12="            << filter12
-                   << "filter13="            << filter13
-                   << "filter14="            << filter14
-                   << "filter15="            << filter15
-                   << "filter16="            << filter16
-                   << "filter17="            << filter17
-                   << "filter18="            << filter18
-                   << "filter19="            << filter19
-                   << "mcp="                 << momentummc
-                   << "mcpt="                << transversemomentummc
-                   << "mceta="               << etamc
-                   << "mcphi="               << phimc
-                   << "mcpdg="               << pdg
-                   << "source="              << fSource
-                   << "px="                  << vx
-                   << "py="                  << vy
-                   << "pz="                  << vz
-                   << "eR="                  << fEr
-                   << "mccharge="            << chargemc
-                   << "signal="              << fSignal
-                   << "\n";
-    } else {
-      if((fFilter2==1) || (fFilter4==1)) fDebugTreee->Fill();
-    }
+    (*fDebugTree) << "PIDdebug"
+                  << "centrality="          << centrality
+                  << "MBtrigger="           << isMBTrigger 
+                  << "CentralTrigger="      << isCentralTrigger
+                  << "SemicentralTrigger="  << isSemicentralTrigger
+                  << "EMCALtrigger="        << isEMCALTrigger
+                 << "run="                 << run
+                 << "eventnb="             << eventnb
+                 << "vx="                  << vtx[0]
+                  << "vy="                  << vtx[1]
+                  << "vz="                  << vtx[2]
+                 << "ncontrib="            << ncontrib
+                 << "id="                  << id
+                 << "dc="                  << doublec
+                  << "p="                   << momentum
+                 << "ptpc="                << momentumTPC
+                 << "pt="                  << transversemomentum
+                 << "eta="                 << eta
+                  << "phi="                 << phi
+                 << "itsrefit="            << itsrefit
+                 << "tpcrefit="            << tpcrefit
+                 << "tofpid="              << tofpid
+                 << "nclustersTPC="        << nclustersTPCfit
+                 << "nclustersTPCall="     << nclustersTPCall
+                  << "nclustersTPCPID="     << nclustersTPCPID
+                  << "nclustersTPCshared="  << nclustersTPCshared
+                  << "ncrossedRowsTPC="     << ncrossedRowsTPC
+                  << "clusterRatioTPC="     << clusterRatioTPCfit
+                 << "clusterRatioTPCall="  << clusterRatioTPCall
+                  << "nclustersITS="        << nclustersITS
+                 << "nclustersTRD="        << nclustersTRD
+                 << "ntrackletsTRD="       << ntrackletsTRDPID
+                 << "nslicesTRD="          << nslicesTRD
+                  << "trd0="                << trdlayer[0]
+                  << "trd1="                << trdlayer[1]
+                  << "trd2="                << trdlayer[2]
+                  << "trd3="                << trdlayer[3]
+                  << "trd4="                << trdlayer[4]
+                  << "trd5="                << trdlayer[5]
+                  << "chi2TRD="             << chi2TRD
+                 << "statusITS0="          << statusL0
+                  << "statusITS1="          << statusL1
+                 << "TOFsigmaEl="          << nSigmaTOF
+                  << "TPCsigmaEl="          << nSigmaTPC
+                 << "TPCsigmaElcorr="      << nSigmaTPCcorr
+                  << "TPCdEdx="             << tPCdEdx
+                  << "TPCdEdxcorr="         << tPCdEdxcorr
+                 << "dcaR="                << dcaxy
+                  << "dcaZ="                << dcaz
+                 << "kinkdaughter="        << kink
+                 << "kinkmother="          << kinkmotherpass
+                 << "nbofmotherkink="      << numberofmotherkink
+                 << "filter0="             << filter0
+                 << "filter1="             << filter1
+                 << "filter2="             << filter2
+                 << "filter3="             << filter3
+                 << "filter4="             << filter4
+                 << "filter5="             << filter5
+                 << "filter6="             << filter6
+                 << "filter7="             << filter7
+                 << "filter8="             << filter8
+                 << "filter9="             << filter9
+                 << "filter10="            << filter10
+                 << "filter11="            << filter11
+                 << "filter12="            << filter12
+                 << "filter13="            << filter13
+                 << "filter14="            << filter14
+                 << "filter15="            << filter15
+                 << "filter16="            << filter16
+                 << "filter17="            << filter17
+                 << "filter18="            << filter18
+                 << "filter19="            << filter19
+                 << "mcp="                 << momentummc
+                  << "mcpt="                << transversemomentummc
+                 << "mceta="               << etamc
+                  << "mcphi="               << phimc
+                  << "mcpdg="               << pdg
+                 << "source="              << source
+                 << "px="                  << vx
+                  << "py="                  << vy
+                  << "pz="                  << vz
+                 << "eR="                  << eR
+                 << "mccharge="            << chargemc
+                 << "signal="              << signal
+                 << "\n";
+
  
     //printf("after\n");
 
@@ -586,8 +606,6 @@ void AliHFEdebugTreeTaskAOD::UserExec(Option_t *){
   arraytrack->~TArrayI();
   fEventNumber++;
 
-  if(!fDebugstream) PostData(1,fDebugTreee);
-
 
 }
 
index 293b7d4..9876111 100644 (file)
@@ -29,7 +29,6 @@ class AliAODMCHeader;
 class TClonesArray;
 class AliHFEsignalCuts;
 class AliHFEextraCuts;
-class TTree;
 
 class AliHFEdebugTreeTaskAOD : public AliAnalysisTaskSE{
   public:
@@ -46,7 +45,6 @@ class AliHFEdebugTreeTaskAOD : public AliAnalysisTaskSE{
     void SetMinNclustersTPC(Int_t mincl) { fNclustersTPC = mincl; };
     void SetMinNclustersTPCPID(Int_t mincl) { fNclustersTPCPID = mincl; };
     void SetMinNclustersITS(Int_t mincl) { fNclustersITS = mincl; };
-    void SetDebugStream(Bool_t on) { fDebugstream = on; };
     AliHFEpidTPC *GetTPCResponse() { return fTPCpid; }
     
   private:
@@ -64,40 +62,9 @@ class AliHFEdebugTreeTaskAOD : public AliAnalysisTaskSE{
     Int_t fNclustersTPCPID;           // Min Number of clusters for TPC PID
     Int_t fNclustersITS;              // Min Number of clusters in ITS
     TString fFilename;                // file name for the debug tree
-    Bool_t  fDebugstream;             // to choose the way
     TTreeSRedirector *fDebugTree;     // Debug Tree
-    TTree *fDebugTreee;               // Debug Tree
-
-    Float_t fCentrality;              // variable
-    Int_t   fRun;                     // run
-    Int_t   fDoublec;                 // double counted
-    Float_t fMomentum;                // Momentum
-    Float_t fMomentumTPC;             // Momentum TPC
-    Float_t fTransverseMomentum;      // Transverse Momentum
-    Float_t fEta;                     // Eta
-    Float_t fPhi;                     // Phi
-    Int_t   fCharge;                  // charge
-    Int_t   fNClustersTPCall;         // Nb of TPC clusters TPC all
-    Int_t   fNClustersTPCPID;         // Nb of TPC clusters TPC PID
-    Int_t   fNClustersTPCshared;      // Nb of TPC clusters shared
-    Int_t   fNCrossedRowsTPC;         // Nb of cross row TPC
-    Float_t fClusterRatioTPCall;      // cls ratio TPC all
-    Int_t   fNClustersITS;            // Nb of ITS clusters
-    Int_t   fStatusL0;                // status L0
-    Int_t   fStatusL1;                // status L1
-    Float_t fSigmaTOF;                // Sigma TOF
-    Float_t fSigmaTPC;                // Sigma TPC
-    Float_t fDcaxy;                   // Dcaxy
-    Float_t fDcaz;                    // Dcaz
-    Int_t   fFilter2;                 // filter2
-    Int_t   fFilter4;                 // filter4
-    Int_t   fSource;                  // source
-    Float_t fEr;                      // er
-    Float_t fSignal;                  // signal
-
-    
   
-    ClassDef(AliHFEdebugTreeTaskAOD, 2)
+    ClassDef(AliHFEdebugTreeTaskAOD, 1)
 };
 #endif
 
index 7a96fe7..df8a654 100644 (file)
@@ -30,6 +30,7 @@
 #include <TString.h>
 #include <TMath.h>
 
+#include "AliAODEvent.h"
 #include "AliAODTrack.h"
 #include "AliAODPid.h"
 #include "AliAODVertex.h"
@@ -44,6 +45,7 @@
 #include "AliVParticle.h"
 #include "AliVertexerTracks.h"
 #include "AliVVertex.h"
+#include "AliExternalTrackParam.h"
 
 #include "AliHFEextraCuts.h"
 
@@ -801,42 +803,38 @@ void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Flo
   //
   // Get impact parameter
   //
+
+  const Double_t kBeampiperadius=3.;
+  Double_t dcaD[2]={-999.,-999.},
+           covD[3]={-999.,-999.,-999.};
   TClass *type = track->IsA();
   if(type == AliESDtrack::Class()){
     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
     esdtrack->GetImpactParameters(radial, z);
   }
   else if(type == AliAODTrack::Class()){
-    AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
-    if(aodtrack){
-      Double_t xyz[3];
-      aodtrack->XYZAtDCA(xyz);
-      //printf("xyz[0] %f, xyz[1] %f, xyz[2] %f\n",xyz[0], xyz[1],xyz[2]);
-      if(xyz[0]==-999. || xyz[1]==-999. ||  xyz[2]==-999.){
-       if(!fEvent) {
-         z = xyz[2];
-         radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
-         //printf("No event\n");
-       }
-       else {
-         // Propagate the global track to the DCA.
-         const AliVVertex *vAOD = fEvent->GetPrimaryVertex();
-         Double_t posAtDCA[2] = {-999,-999};
-         Double_t covar[3] = {-999,-999,-999};
-         AliAODTrack *copyaodtrack = new AliAODTrack(*aodtrack);
-         if(copyaodtrack && copyaodtrack->PropagateToDCA(vAOD,fEvent->GetMagneticField(),100.,posAtDCA,covar)) {
-           z = posAtDCA[1];
-           radial = posAtDCA[0];
-           //printf("Propagate z %f and radial %f\n",z,radial);
-         }
-         delete copyaodtrack;
-       }
-      }
-      else {
-       z = xyz[2];
-       radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
-      }
+
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
     }
+
+    //Case ESD track: take copy constructor
+    AliAODTrack *aodtrack = NULL;
+    AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
+    if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
+
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+    AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
+    if(etp.PropagateToDCA(vtxAODSkip,aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
+      radial = dcaD[0];
+      z = dcaD[1];
+    }
+    //if(vtxAODSkip) delete vtxAODSkip;
+    if(aodtrack) delete aodtrack;
   }
 }
 //______________________________________________________
@@ -907,7 +905,6 @@ Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
 
 }
 
-
 //______________________________________________________
 Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
   //
@@ -930,85 +927,73 @@ Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
 }
 //______________________________________________________
 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
-       //
-       // Get HFE impact parameter (with recalculated primary vertex)
-       //
-       dcaxy=0;
-       dcansigmaxy=0;
+  //
+  // Get HFE impact parameter (with recalculated primary vertex)
+  //
+  dcaxy=0;
+  dcansigmaxy=0;
   if(!fEvent){
     AliDebug(1, "No Input event available\n");
     return;
   }
-  const Double_t kBeampiperadius=3.;
   TString type = track->IsA()->GetName();
-  Double_t dcaD[2]={-999.,-999.};
-  Double_t covD[3]={-999.,-999.,-999.};
-  Float_t dcaF[2]={-999.,-999.};
-  Float_t covF[3]={-999.,-999.,-999.};
-  
-  AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
-  if(!esdevent) {
-    AliDebug(1, "No esd event available\n");
-    return;
-  }
-  const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
-  if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
-   // recalculate primary vertex for peri. and pp
-   AliVertexerTracks vertexer(fEvent->GetMagneticField());
-   vertexer.SetITSMode();
-   vertexer.SetMinClusters(4);
-        Int_t skipped[2];
-   skipped[0] = track->GetID();
-   vertexer.SetSkipTracks(1,skipped);
-   
- //----diamond constraint-----------------------------
-   vertexer.SetConstraintOn();
-   Float_t diamondcovxy[3];
-   esdevent->GetDiamondCovXY(diamondcovxy);
-   Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
-   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
-   AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
-   vertexer.SetVtxStart(diamond);
-   delete diamond; diamond=NULL;
- //----------------------------------------------------
-    
-   vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
-
-   // Getting the DCA
-   // Propagation always done on a working copy to not disturb the track params of the original track
-   AliESDtrack *esdtrack = NULL;
-   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
-     // Case ESD track: take copy constructor
-     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
-     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
-   } else {
-    // Case AOD track: take different constructor
-    esdtrack = new AliESDtrack(track);
-   }
-   if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
-    // protection
-    dcaxy = dcaD[0];
-    if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
-    if(!covD[0]) dcansigmaxy = -49.;
-   }
-   delete esdtrack;
-   delete vtxESDSkip;
-  } 
-  else{
-   AliESDtrack *esdtrack = NULL;
-   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
-    // Case ESD track: take copy constructor
+  const Double_t kBeampiperadius=3.;
+  Double_t dcaD[2]={-999.,-999.},
+           covD[3]={-999.,-999.,-999.};
+
+  if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
+    //case of ESD tracks
+    AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
+    if(!esdevent) {
+      AliDebug(1, "No esd event available\n");
+      return;
+    }
+
+    //case ESD track: take copy constructor
+    AliESDtrack *esdtrack = NULL;
     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
-   } else {
-    // Case AOD track: take different constructor
-    esdtrack = new AliESDtrack(track);
-   }
-   if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF); 
-   dcaxy = dcaF[0];
-   if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
-   if(!covF[0]) dcansigmaxy = -49.;
-   delete esdtrack;
+
+    const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+    if(!vtxESDSkip) return;
+    if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+    }
+
+    if(!vtxESDSkip) return;
+    if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
+      dcaxy = dcaD[0];
+      if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+    }
+    //delete vtxESDSkip;
+    delete esdtrack;
+  }
+  else {
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
+    }
+
+    //Case ESD track: take copy constructor
+    AliAODTrack *aodtrack = NULL;
+    AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
+    if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
+
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+    if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+    } 
+    if(!vtxAODSkip) return;
+    AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
+    if(etp.PropagateToDCA(vtxAODSkip,aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
+      dcaxy = dcaD[0];
+      if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+    }
+    //if(vtxAODSkip) delete vtxAODSkip;
+    if(aodtrack) delete aodtrack;
   }
 }
 
@@ -1023,69 +1008,56 @@ void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2],
   }
   const Double_t kBeampiperadius=3.;
   TString type = track->IsA()->GetName();
-  Float_t dcaF[2]={-999.,-999.};
-  Float_t covF[3]={-999.,-999.,-999.};
   
-  AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
-  if(!esdevent) {
-    AliDebug(1, "No esd event available\n");
-    return;
-  }
-  const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
-  if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
-   // recalculate primary vertex for peri. and pp
-   AliVertexerTracks vertexer(fEvent->GetMagneticField());
-   vertexer.SetITSMode();
-   vertexer.SetMinClusters(4);
-        Int_t skipped[2];
-   skipped[0] = track->GetID();
-   vertexer.SetSkipTracks(1,skipped);
-   
- //----diamond constraint-----------------------------
-   vertexer.SetConstraintOn();
-   Float_t diamondcovxy[3];
-   esdevent->GetDiamondCovXY(diamondcovxy);
-   Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
-   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
-   AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
-   vertexer.SetVtxStart(diamond);
-   delete diamond; diamond=NULL;
- //----------------------------------------------------
-    
-   vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
-
-   // Getting the DCA
-   // Propagation always done on a working copy to not disturb the track params of the original track
-   AliESDtrack *esdtrack = NULL;
-   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
-     // Case ESD track: take copy constructor
-     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
-     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
-   } else {
-    // Case AOD track: take different constructor
-    esdtrack = new AliESDtrack(track);
-   }
-   if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
-   delete esdtrack;
-   delete vtxESDSkip;
-  } 
-  else{
-   AliESDtrack *esdtrack = NULL;
-   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
-    // Case ESD track: take copy constructor
+  if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
+    //case of ESD tracks
+    AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
+    if(!esdevent) {
+      AliDebug(1, "No esd event available\n");
+      return;
+    }
+
+    //case ESD track: take copy constructor
+    AliESDtrack *esdtrack = NULL;
     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
-   } else {
-    // Case AOD track: take different constructor
-    esdtrack = new AliESDtrack(track);
-   }
-   if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF); 
-   dcaD[0] = dcaF[0];
-   dcaD[1] = dcaF[1];
-   covD[0] = covF[0];
-   covD[1] = covF[1];
-   covD[2] = covF[2];
-   delete esdtrack;
+
+    const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+    if(!vtxESDSkip) return;
+    if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+    }
+    if(!vtxESDSkip) return;
+
+    if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
+
+    //delete vtxESDSkip;
+    delete esdtrack;
+  }
+  else {
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
+    }
+
+    //Case ESD track: take copy constructor
+    AliAODTrack *aodtrack = NULL;
+    AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
+    if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
+
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+    if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+    }
+    if(!vtxAODSkip) return;
+    AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
+    etp.PropagateToDCA(vtxAODSkip,aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD);
+
+    //delete vtxAODSkip;
+    delete aodtrack;
   }
 }
 
@@ -1158,3 +1130,101 @@ Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
   }
   return patternOK;
 }
+
+//---------------------------------------------------------------------------
+const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(AliESDEvent *esdevent, AliVTrack *track) {
+  //
+  // This method returns a primary vertex without the daughter tracks of the 
+  // candidate and it recalculates the impact parameters and errors for ESD tracks.
+  // 
+  // The output vertex is created with "new". 
+  //
+
+  const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+
+  AliVertexerTracks vertexer(fEvent->GetMagneticField());
+  vertexer.SetITSMode();
+  vertexer.SetMinClusters(4);
+  Int_t skipped[2];
+  skipped[0] = track->GetID();
+  vertexer.SetSkipTracks(1,skipped);
+
+  //diamond constraint
+  vertexer.SetConstraintOn();
+  Float_t diamondcovxy[3];
+  esdevent->GetDiamondCovXY(diamondcovxy);
+  Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
+  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
+  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
+  vertexer.SetVtxStart(diamond);
+  delete diamond; diamond=NULL;
+
+  vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
+
+  return vtxESDSkip;
+}
+
+//---------------------------------------------------------------------------
+AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod, AliVTrack *track) {
+  //
+  // This method returns a primary vertex without the daughter tracks of the 
+  // candidate and it recalculates the impact parameters and errors for AOD tracks.
+  // The output vertex is created with "new". 
+  //
+
+  AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
+  if(!vtxAOD) return 0;
+  TString title=vtxAOD->GetTitle();
+  if(!title.Contains("VertexerTracks")) return 0;
+
+  AliVertexerTracks *vertexer = new AliVertexerTracks(aod->GetMagneticField());
+
+  vertexer->SetITSMode();
+  vertexer->SetMinClusters(3);
+  vertexer->SetConstraintOff();
+
+  if(title.Contains("WithConstraint")) {
+    Float_t diamondcovxy[3];
+    aod->GetDiamondCovXY(diamondcovxy);
+    Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
+    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
+    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
+    vertexer->SetVtxStart(diamond);
+    delete diamond; diamond=NULL;
+  }
+  Int_t skipped[2]; for(Int_t i=0;i<2;i++) skipped[i]=-1;
+  Int_t id = (Int_t)track->GetID();
+  if(!(id<0)) skipped[0] = id;
+
+  /*// exclude tracks with global constrained parameters
+  Int_t nTracks=aod->GetNumberOfTracks();
+  for(Int_t i=0; i<nTracks; i++){
+    t = aod->GetTrack(i);
+    if(t->TestFilterMask(512)){
+      id = (Int_t)t->GetID();
+      skipped[nTrksToSkip++] = id;
+    }
+  }*/
+
+  vertexer->SetSkipTracks(1,skipped);
+  AliESDVertex *vtxESDNew = vertexer->FindPrimaryVertex(aod);
+
+  delete vertexer; vertexer=NULL;
+
+  if(!vtxESDNew) return 0;
+  if(vtxESDNew->GetNContributors()<=0) {
+    delete vtxESDNew; vtxESDNew=NULL;
+    return 0;
+  }
+
+  // convert to AliAODVertex
+  Double_t pos[3],cov[6],chi2perNDF;
+  vtxESDNew->GetXYZ(pos); // position
+  vtxESDNew->GetCovMatrix(cov); //covariance matrix
+  chi2perNDF = vtxESDNew->GetChi2toNDF();
+  delete vtxESDNew; vtxESDNew=NULL;
+
+  AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
+
+  return vtxAODNew;
+}
index 02b0354..7b1220a 100644 (file)
@@ -30,6 +30,10 @@ class TList;
 class AliVEvent;
 class AliVParticle;
 class AliVTrack;
+class AliVVertex;
+class AliAODVertex;
+class AliAODEvent;
+class AliESDEvent;
 
 class AliHFEextraCuts: public AliCFCutBase{
   public:
@@ -96,6 +100,8 @@ class AliHFEextraCuts: public AliCFCutBase{
     void GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy); // temporary moved from protected to publich for IP QA 
     void GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]);
     void GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z);
+    const AliVVertex* RemoveDaughtersFromPrimaryVtx(AliESDEvent *esdevent, AliVTrack *track);
+    AliAODVertex* RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod, AliVTrack *track);
     Int_t GetITSstatus(const AliVTrack * const track, Int_t layer) const;
     Bool_t CheckITSstatus(Int_t itsStatus) const;
     Bool_t CheckITSpattern(const AliVTrack *const track) const;
index 6653532..88245d7 100644 (file)
@@ -1939,7 +1939,7 @@ Int_t AliHFEmcQA::GetElecSource(const AliAODMCParticle * const mcpart)
     // J/spi
     //
     Int_t jLabel = partMother->GetMother();
-    if ((jLabel<0) || (jLabel>=fMCArray->GetEntriesFast())){
+    if ((jLabel>=0) && (jLabel<fMCArray->GetEntriesFast())){
       if((mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(TMath::Abs(jLabel))))){
        Int_t grandMaPDG = mctrack->GetPdgCode();
        if((int(abs(grandMaPDG)/100.)%10) == kBeauty || (int(abs(grandMaPDG)/1000.)%10) == kBeauty) {
index 12ce8ce..5e0d177 100644 (file)
@@ -33,6 +33,7 @@ AliHFEpidObject &AliHFEpidObject::operator=(const AliHFEpidObject &ref){
     fAnalysisType = ref.fAnalysisType;
     fAbInitioPID = ref.fAbInitioPID;
     fCentrality = ref.fCentrality;
+    fMultiplicity = ref.fMultiplicity;
   }
   return *this;
 }
index 6eff82a..8dbfc14 100644 (file)
@@ -35,6 +35,7 @@ class AliHFEpidObject{
       fAnalysisType(kESDanalysis),
       fAbInitioPID(-1),
       fCentrality(99),
+      fMultiplicity(0),
       fIsPbPb(kFALSE)         // Default: pp
       {
       }
@@ -43,6 +44,7 @@ class AliHFEpidObject{
       fAnalysisType(ref.fAnalysisType),
       fAbInitioPID(ref.fAbInitioPID),
       fCentrality(ref.fCentrality),
+      fMultiplicity(ref.fMultiplicity),
       fIsPbPb(ref.fIsPbPb)
       {
       }
@@ -54,12 +56,14 @@ class AliHFEpidObject{
     void SetAnalysisType(AnalysisType_t type) { fAnalysisType = type; }
     void SetAbInitioPID(Int_t abInitioPID) { fAbInitioPID = abInitioPID; }
     void SetCentrality(Int_t centrality) { fCentrality = centrality; }
+    void SetMulitplicity(Double_t mult) { fMultiplicity = mult; }
     void SetPbPb() { fIsPbPb = kTRUE; }
     void SetPP() { fIsPbPb = kFALSE; }
 
     const AliVTrack *GetRecTrack() const { return fkRecTrack; }
     Int_t GetAbInitioPID() const { return fAbInitioPID; }
     Int_t GetCentrality() const { return fCentrality; }
+    Double_t GetMultiplicity() const { return fMultiplicity; }
     Bool_t IsAODanalysis() const { return fAnalysisType == static_cast<UChar_t>(kAODanalysis); }
     Bool_t IsESDanalysis() const { return fAnalysisType == static_cast<UChar_t>(kESDanalysis); }
     Bool_t IsPbPb() const { return fIsPbPb; }
@@ -69,6 +73,7 @@ class AliHFEpidObject{
     UChar_t fAnalysisType;              // Analysis Mode (ESD or AOD)
     Int_t fAbInitioPID;                 // AbInitio PID
     Int_t fCentrality;                  // Centrality Information
+    Double_t fMultiplicity;             // Multiplicity information
     Bool_t fIsPbPb;                     // Collision type
 };
 #endif
index f22faa4..dc8b67b 100644 (file)
@@ -148,6 +148,21 @@ Int_t AliHFEpidTOF::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager
   return pdg;
 }
 //___________________________________________________________________
+void AliHFEpidTOF::SetTOFnSigmaBand(Float_t lower, Float_t upper)
+{
+  //
+  // Lower and higher cut independant of the centrality
+  //
+
+  for(Int_t k=0; k < 12; k++) {
+    fSigmaBordersTOFLower[k] = lower;
+    fSigmaBordersTOFUpper[k] = upper;
+  }
+
+  SetBit(kSigmaBand, kTRUE);
+
+}
+//___________________________________________________________________
 void AliHFEpidTOF::SetTOFnSigmaBandCentrality(Float_t lower, Float_t upper, Int_t centralityBin)
 {
   //
index 185b45e..51111d5 100644 (file)
@@ -30,7 +30,7 @@ class AliHFEpidTOF : public AliHFEpidBase{
     virtual Int_t     IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *piqa) const;
   
     void SetTOFnSigma(Float_t nSigma) { fNsigmaTOF = nSigma; };
-    void SetTOFnSigmaBand(Float_t lower, Float_t upper) { fSigmaBordersTOFLower[0] = lower; fSigmaBordersTOFUpper[0] = upper; SetBit(kSigmaBand, kTRUE); }
+    void SetTOFnSigmaBand(Float_t lower, Float_t upper);
     void SetTOFnSigmaBandCentrality(Float_t lower, Float_t upper, Int_t centralityBin); 
 
   protected:
index f8678bf..9ed1c1e 100644 (file)
@@ -50,6 +50,7 @@ AliHFEpidTPC::AliHFEpidTPC() :
   AliHFEpidBase()
   , fLineCrossingsEnabled(0)
   , fkEtaCorrection(NULL)
+  , fkCentralityCorrection(NULL)
   , fHasCutModel(kFALSE)
   , fUseOnlyOROC(kFALSE)
   , fNsigmaTPC(3)
@@ -75,6 +76,7 @@ AliHFEpidTPC::AliHFEpidTPC(const char* name) :
   AliHFEpidBase(name)
   , fLineCrossingsEnabled(0)
   , fkEtaCorrection(NULL)
+  , fkCentralityCorrection(NULL)
   , fHasCutModel(kFALSE)
   , fUseOnlyOROC(kFALSE)
   , fNsigmaTPC(3)
@@ -98,6 +100,7 @@ AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
   AliHFEpidBase("")
   , fLineCrossingsEnabled(0)
   , fkEtaCorrection(NULL)
+  , fkCentralityCorrection(NULL)
   , fHasCutModel(ref.fHasCutModel)
   , fUseOnlyOROC(ref.fUseOnlyOROC)
   , fNsigmaTPC(2)
@@ -127,7 +130,8 @@ void AliHFEpidTPC::Copy(TObject &o) const{
   //
   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
 
-  target.fkEtaCorrection =fkEtaCorrection;
+  target.fkEtaCorrection = fkEtaCorrection;
+  target.fkCentralityCorrection = fkCentralityCorrection;
   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
   target.fHasCutModel = fHasCutModel;
   target.fUseOnlyOROC = fUseOnlyOROC;
@@ -178,7 +182,7 @@ Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager
   const AliVTrack *rectrack;
   AliESDtrack esdtrack;
   AliAODTrack aodtrack;
-  if(fUseOnlyOROC && !fkEtaCorrection) {
+  /*if(fUseOnlyOROC && !(fkEtaCorrection || fkCentralityCorrection)) {
     if(track->IsESDanalysis()){
       esdtrack.~AliESDtrack();
       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
@@ -191,18 +195,25 @@ Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager
       rectrack = &aodtrack;
     }
   }
-  else if(fkEtaCorrection){
+  else if(fkEtaCorrection || fkCentralityCorrection){*/
+  if(fkEtaCorrection || fkCentralityCorrection){
     // Correction available
     // apply it on copy
     if(track->IsESDanalysis()){
       esdtrack.~AliESDtrack();
       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
-      ApplyEtaCorrection(&esdtrack, anatype);
+      if(track->IsPbPb() && HasCentralityCorrection())
+        ApplyCentralityCorrection(&esdtrack, track->GetMultiplicity(), anatype);
+      if(HasEtaCorrection())
+        ApplyEtaCorrection(&esdtrack, anatype);
       rectrack = &esdtrack;
     } else {
       aodtrack.~AliAODTrack();
       new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
-      ApplyEtaCorrection(&aodtrack, anatype);
+      if(track->IsPbPb() && HasCentralityCorrection())
+        ApplyCentralityCorrection(&aodtrack, track->GetMultiplicity(), anatype);
+      if(HasEtaCorrection())
+        ApplyEtaCorrection(&aodtrack, anatype);
       rectrack = &aodtrack;
     }
   } else {
@@ -311,13 +322,32 @@ void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::Analysi
            eta = track->Eta();
   if(anatype == AliHFEpidObject::kESDanalysis){
     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
-    esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
+    if(fkEtaCorrection->Eval(eta)>0.0) esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
   } else {
     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
     AliAODPid *pid = aodtrack->GetDetPid();
-    if(pid) pid->SetTPCsignal(original);
+    if(pid && fkEtaCorrection->Eval(eta)>0.0) pid->SetTPCsignal(original/fkEtaCorrection->Eval(eta));
   }
 }
+
+//___________________________________________________________________
+void AliHFEpidTPC::ApplyCentralityCorrection(AliVTrack *track, Double_t centralityEstimator, AliHFEpidObject::AnalysisType_t anatype) const{
+  //
+  // Apply correction for the eta dependence
+  // N.B. This correction has to be applied on a copy track
+  //
+  AliDebug(1, Form("Applying correction function %s\n", fkCentralityCorrection->GetName()));
+  Double_t original = track->GetTPCsignal();
+  if(anatype == AliHFEpidObject::kESDanalysis){
+    AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+    if(fkCentralityCorrection->Eval(centralityEstimator)>0.0) esdtrack->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
+  } else {
+    AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
+    AliAODPid *pid = aodtrack->GetDetPid();
+    if(pid && fkCentralityCorrection->Eval(centralityEstimator)>0.0) pid->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator));
+  }
+}
+
 //___________________________________________________________________
 void AliHFEpidTPC::UseOROC(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
   //
index 69de5a1..6210a89 100644 (file)
@@ -61,10 +61,13 @@ class AliHFEpidTPC : public AliHFEpidBase{
     void SetLowerSigmaCutDefault(const TF1 * const model) { fkLowerSigmaCut[0] = model; fHasCutModel = kTRUE; }
     void SetLowerSigmaCutCentrality(const TF1 * const model, Int_t centralityBin) { if(centralityBin < 11) fkLowerSigmaCut[centralityBin+1] = model; fHasCutModel = kTRUE; }
     void SetEtaCorrection(const TF1 *const param) { fkEtaCorrection = param; }
+    void SetCentralityCorrection(const TF1 *const param){ fkCentralityCorrection = param; }
     Bool_t HasEtaCorrection() const { return fkEtaCorrection != NULL; }
+    Bool_t HasCentralityCorrection() const { return fkCentralityCorrection != NULL; } 
 
     Double_t GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const;
     void ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const;
+    void ApplyCentralityCorrection(AliVTrack *track, Double_t centralityEstimator, AliHFEpidObject::AnalysisType_t anatype) const;
     void UseOROC(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const;
 
   protected:
@@ -83,6 +86,7 @@ class AliHFEpidTPC : public AliHFEpidBase{
     const TF1 *fkUpperSigmaCut[12];                         // Upper Sigma Cut
     const TF1 *fkLowerSigmaCut[12];                         // Lower Sigma Cut
     const TF1 *fkEtaCorrection;                             // Correction for the eta dependence
+    const TF1 *fkCentralityCorrection;                      // Correction for the centrality dependence
     Bool_t fHasCutModel;                                    // Has cut model functions
     Bool_t fUseOnlyOROC;                                    // Use only OROC
     Float_t fPAsigCut[2];                                   // Momentum region where to perform asymmetric sigma cut
index f4e4834..9a15643 100644 (file)
@@ -117,7 +117,7 @@ void AliHFEsignalCuts::SetMCAODInfo(TClonesArray *mcarray){
 }
 
 //____________________________________________________________
-Bool_t AliHFEsignalCuts::IsSelected(TObject *o) {
+Bool_t AliHFEsignalCuts::IsSelected(TObject *o){
   //
   // Define signal as electron coming from charm or beauty
   // @TODO: Implement setter so that also either of them can be defined
diff --git a/PWGHF/hfe/AliHFEsmearDCA.cxx b/PWGHF/hfe/AliHFEsmearDCA.cxx
new file mode 100644 (file)
index 0000000..9a49e54
--- /dev/null
@@ -0,0 +1,424 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Get improved dca info. by ITS upgrade implemented by the ALICE Heavy Flavour Electron Group
+//
+// Authors:
+//   MinJung Kweon <minjung@physi.uni-heidelberg.de>
+//
+#include <TBits.h>
+#include <TClass.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TList.h>
+#include <TString.h>
+#include <TMath.h>
+#include <TGraph.h>
+#include <TFile.h>
+
+#include "AliAODTrack.h"
+#include "AliAODPid.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
+#include "AliESDtrack.h"
+#include "AliLog.h"
+#include "AliMCParticle.h"
+#include "AliVEvent.h"
+#include "AliVTrack.h"
+#include "AliVParticle.h"
+#include "AliVertexerTracks.h"
+#include "AliVVertex.h"
+#include "AliExternalTrackParam.h"
+#include "AliMCEvent.h"
+
+#include "AliHFEsmearDCA.h"
+
+ClassImp(AliHFEsmearDCA)
+
+
+//______________________________________________________
+AliHFEsmearDCA::AliHFEsmearDCA(const Char_t */*name*/, const char *resfileCurURI, const char *resfileUpgURI, const Char_t */*title*/):
+   fEvent(NULL),
+   fMCEvent(NULL),
+   fD0ZResPCur  (0),
+   fD0ZResKCur  (0),
+   fD0ZResPiCur (0),
+   fD0ZResECur  (0),
+   fD0RPResPCur (0),
+   fD0RPResKCur (0),
+   fD0RPResPiCur(0),
+   fD0RPResECur (0),
+   fPt1ResPCur  (0),
+   fPt1ResKCur  (0),
+   fPt1ResPiCur (0),
+   fPt1ResECur  (0),
+   fD0ZResPUpg  (0),
+   fD0ZResKUpg  (0),
+   fD0ZResPiUpg (0),
+   fD0ZResEUpg  (0),
+   fD0RPResPUpg (0),
+   fD0RPResKUpg (0),
+   fD0RPResPiUpg(0),
+   fD0RPResEUpg (0),
+   fPt1ResPUpg  (0),
+   fPt1ResKUpg  (0),
+   fPt1ResPiUpg (0),
+   fPt1ResEUpg  (0)
+{
+  //
+  // Constructor to be used to create the task.
+  // The the URIs specify the resolution files to be used. 
+  // They are expected to contain TGraphs with the resolutions
+  // for the current and the upgraded ITS (see code for details).
+  // One may also specify for how many tracks debug information
+  // is written to the output.
+  //
+  TFile *resfileCur=TFile::Open(resfileCurURI);
+  fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
+  fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
+  fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
+  fD0RPResECur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResE")));
+  fD0ZResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP"  )));
+  fD0ZResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK"  )));
+  fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
+  fD0ZResECur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResE")));
+  fPt1ResPCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP"  )));
+  fPt1ResKCur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK"  )));
+  fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
+  fPt1ResECur  =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResE" )));
+  fD0RPResPCur ->SetName("D0RPResPCur" );
+  fD0RPResKCur ->SetName("D0RPResKCur" );
+  fD0RPResPiCur->SetName("D0RPResPiCur");
+  fD0RPResECur ->SetName("D0RPResECur");
+  fD0ZResPCur  ->SetName("D0ZResPCur"  ); 
+  fD0ZResKCur  ->SetName("D0ZResKCur"  );
+  fD0ZResPiCur ->SetName("D0ZResPiCur" );
+  fD0ZResECur  ->SetName("D0ZResECur" );
+  fPt1ResPCur  ->SetName("Pt1ResPCur"  );
+  fPt1ResKCur  ->SetName("Pt1ResKCur"  );
+  fPt1ResPiCur ->SetName("Pt1ResPiCur" );
+  fPt1ResECur  ->SetName("Pt1ResECur" );
+  delete resfileCur;
+  TFile *resfileUpg=TFile::Open(resfileUpgURI );
+  fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
+  fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
+  fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
+  fD0RPResEUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResE")));
+  fD0ZResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP"  )));
+  fD0ZResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK"  )));
+  fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
+  fD0ZResEUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResE" )));
+  fPt1ResPUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP"  )));
+  fPt1ResKUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK"  )));
+  fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
+  fPt1ResEUpg  =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResE" )));
+  fD0RPResPUpg ->SetName("D0RPResPUpg" );
+  fD0RPResKUpg ->SetName("D0RPResKUpg" );
+  fD0RPResPiUpg->SetName("D0RPResPiUpg");
+  fD0RPResEUpg ->SetName("D0RPResEUpg");
+  fD0ZResPUpg  ->SetName("D0ZResPUpg"  );
+  fD0ZResKUpg  ->SetName("D0ZResKUpg"  );
+  fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
+  fD0ZResEUpg  ->SetName("D0ZResEUpg" );
+  fPt1ResPUpg  ->SetName("Pt1ResPUpg"  );
+  fPt1ResKUpg  ->SetName("Pt1ResKUpg"  );
+  fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
+  fPt1ResEUpg  ->SetName("Pt1ResEUpg" );
+  delete resfileUpg;
+
+}
+
+//______________________________________________________
+AliHFEsmearDCA::AliHFEsmearDCA(const AliHFEsmearDCA &c):
+   TObject(c),
+   fEvent(c.fEvent),
+   fMCEvent(c.fMCEvent),
+   fD0ZResPCur  (c.fD0ZResPCur),
+   fD0ZResKCur  (c.fD0ZResKCur),
+   fD0ZResPiCur (c.fD0ZResPiCur),
+   fD0ZResECur  (c.fD0ZResECur),
+   fD0RPResPCur (c.fD0RPResPCur),
+   fD0RPResKCur (c.fD0RPResKCur),
+   fD0RPResPiCur(c.fD0RPResPiCur),
+   fD0RPResECur (c.fD0RPResECur),
+   fPt1ResPCur  (c.fPt1ResPCur),
+   fPt1ResKCur  (c.fPt1ResKCur),
+   fPt1ResPiCur (c.fPt1ResPiCur),
+   fPt1ResECur  (c.fPt1ResECur),
+   fD0ZResPUpg  (c.fD0ZResPUpg),
+   fD0ZResKUpg  (c.fD0ZResKUpg),
+   fD0ZResPiUpg (c.fD0ZResPiUpg),
+   fD0ZResEUpg  (c.fD0ZResEUpg),
+   fD0RPResPUpg (c.fD0RPResPUpg),
+   fD0RPResKUpg (c.fD0RPResKUpg),
+   fD0RPResPiUpg(c.fD0RPResPiUpg),
+   fD0RPResEUpg (c.fD0RPResEUpg),
+   fPt1ResPUpg  (c.fPt1ResPUpg),
+   fPt1ResKUpg  (c.fPt1ResKUpg),
+   fPt1ResPiUpg (c.fPt1ResPiUpg),
+   fPt1ResEUpg  (c.fPt1ResEUpg)
+{
+  //
+  // Copy constructor
+  // Performs a deep copy
+  //
+}
+
+//______________________________________________________
+AliHFEsmearDCA::~AliHFEsmearDCA(){
+  //
+  // Destructor
+  //
+}
+
+//______________________________________________________
+void AliHFEsmearDCA::SetRecEventInfo(const TObject *event){
+  //
+  // Set Virtual event an make a copy
+  //
+  if (!event) {
+    AliError("Pointer to AliVEvent !");
+    return;
+  }
+  TString className(event->ClassName());
+  if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
+    AliError("argument must point to an AliESDEvent or AliAODEvent !");
+    return ;
+  }
+  fEvent = (AliVEvent*) event;
+
+}
+
+//______________________________________________________
+void AliHFEsmearDCA::GetImproveITSImpactParameters(AliVTrack *track, Double_t &dcaxyn, Double_t &dcaxyo, Double_t &dcaxySign, Double_t &dcaxySigo, Double_t &dcazn, Double_t &dcazo, Double_t &dcazSign, Double_t &dcazSigo){
+  //
+  // Get HFE impact parameter (with recalculated primary vertex)
+  //
+  if(!fEvent){
+    AliDebug(1, "No Input event available\n");
+    return;
+  }
+  const Double_t kBeampiperadius=3.;
+  const Double_t kBeampiperadiusnew=1.9;
+  TString type = track->IsA()->GetName();
+  Double_t dcan[2]={-999.,-999.};
+  Double_t covn[3]={-999.,-999.,-999.};
+  Double_t dcao[2]={-999.,-999.};
+  Double_t covo[3]={-999.,-999.,-999.};
+
+  // Get reconstructed track parameters
+
+  Double_t xyz[3],pxpypz[3],cv[21];
+  track->GetXYZ(xyz);
+  pxpypz[0]=track->Px();
+  pxpypz[1]=track->Py(); 
+  pxpypz[2]=track->Pz(); 
+  track->GetCovarianceXYZPxPyPz(cv);
+  Short_t sign = (Short_t)track->Charge();
+
+  //AliExternalTrackParam et(track); //it doesn't work, i don't know why yet! so, I use the line below
+  AliExternalTrackParam et(xyz,pxpypz,cv,sign);
+  Double_t *param=const_cast<Double_t*>(et.GetParameter());
+
+  const AliMCParticle *mc=static_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
+  Double_t mcx[3];
+  Double_t mcp[3];
+  Double_t mccv[36]={0.};
+  Short_t  mcc;
+  mc->XvYvZv(mcx);
+  mc->PxPyPz(mcp);
+  mcc=mc->Charge();
+  AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
+  const Double_t *parammc=mct.GetParameter();
+  const AliVVertex *tmpvtx = fMCEvent->GetPrimaryVertex();
+
+  AliVertex vtx(mcx,1.,1);
+  AliVertex *tmpvtx2 = (AliVertex *)fMCEvent->GetPrimaryVertex();
+
+  mct.PropagateToDCA(tmpvtx2,track->GetBz(),10.); //mjtest
+  // Correct reference points and frames according to MC
+  // TODO: B-Field correct?
+  // TODO: failing propagation....
+  et.PropagateToDCA(tmpvtx2,track->GetBz(),10.);
+  et.Rotate(mct.GetAlpha());
+
+
+  // Select appropriate smearing functions
+  Double_t ptmc=TMath::Abs(mc->Pt());
+  Double_t sd0rpn=0.;
+  Double_t sd0zn =0.;
+  Double_t spt1n =0.;
+  Double_t sd0rpo=0.;
+  Double_t sd0zo =0.;
+  Double_t spt1o =0.;
+
+  switch (mc->Particle()->GetPdgCode()) {
+  case 2212: case -2212:
+    sd0rpo=EvalGraph(fD0RPResPCur,ptmc);
+    sd0zo =EvalGraph(fD0ZResPCur ,ptmc);
+    spt1o =EvalGraph(fPt1ResPCur ,ptmc);
+    sd0rpn=EvalGraph(fD0RPResPUpg,ptmc);
+    sd0zn =EvalGraph(fD0ZResPUpg ,ptmc);
+    spt1n =EvalGraph(fPt1ResPUpg ,ptmc);
+    break;
+  case 321: case -321:
+    sd0rpo=EvalGraph(fD0RPResKCur,ptmc); 
+    sd0zo =EvalGraph(fD0ZResKCur ,ptmc);
+    spt1o =EvalGraph(fPt1ResKCur ,ptmc);
+    sd0rpn=EvalGraph(fD0RPResKUpg,ptmc);
+    sd0zn =EvalGraph(fD0ZResKUpg ,ptmc);
+    spt1n =EvalGraph(fPt1ResKUpg ,ptmc);
+    break;
+  case 211: case -211:
+    sd0rpo=EvalGraph(fD0RPResPiCur,ptmc); 
+    sd0zo =EvalGraph(fD0ZResPiCur ,ptmc);
+    spt1o =EvalGraph(fPt1ResPiCur ,ptmc);
+    sd0rpn=EvalGraph(fD0RPResPiUpg,ptmc);
+    sd0zn =EvalGraph(fD0ZResPiUpg ,ptmc);
+    spt1n =EvalGraph(fPt1ResPiUpg ,ptmc);
+    break;
+  case 11: case -11:
+    sd0rpo=EvalGraph(fD0RPResECur,ptmc); 
+    sd0zo =EvalGraph(fD0ZResECur ,ptmc);
+    spt1o =EvalGraph(fPt1ResECur ,ptmc);
+    sd0rpn=EvalGraph(fD0RPResEUpg,ptmc);
+    sd0zn =EvalGraph(fD0ZResEUpg ,ptmc);
+    spt1n =EvalGraph(fPt1ResEUpg ,ptmc);
+    break;
+  default:
+    return;
+  }
+
+  //mj special to check resolution smeared by 10 %
+  sd0rpn = sd0rpo * 1.1;
+  sd0zn = sd0zo * 1.1;
+  spt1n = spt1o;
+
+  // Use the same units (i.e. cm and GeV/c)! TODO: pt!
+  sd0rpo*=1.e-4;
+  sd0zo *=1.e-4;
+  sd0rpn*=1.e-4;
+  sd0zn *=1.e-4;
+
+  // Apply the smearing
+  Double_t d0zo  =param  [1];
+  Double_t d0zmc =parammc[1];
+  Double_t d0rpo =param  [0];
+  Double_t d0rpmc=parammc[0];
+  Double_t pt1o  =param  [4];
+  Double_t pt1mc =parammc[4];
+  Double_t dd0zo =d0zo-d0zmc;
+  Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
+  Double_t d0zn  =d0zmc+dd0zn;
+  Double_t dd0rpo=d0rpo-d0rpmc;
+  Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
+  Double_t d0rpn =d0rpmc+dd0rpn;
+  Double_t dpt1o =pt1o-pt1mc;
+  Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
+  Double_t pt1n  =pt1mc+dpt1n;
+
+  param[0]=d0rpn;
+  param[1]=d0zn ;
+  param[4]=pt1n ;
+
+  //printf("d0rpo= %lf d0rpmc= %lf dd0rpo= %lf dd0rpn= %lf\n",d0rpo,d0rpmc,dd0rpo,dd0rpn);
+  // Copy the smeared parameters to the AOD track
+  Double_t x[3];
+  Double_t p[3];
+  et.GetXYZ(x);
+  et.GetPxPyPz(p);
+//  printf("x[0]= %lf x[1]= %lf x[2]= %lf \n",x[0],x[1],x[2]);
+//  track->SetPosition(x,kFALSE);
+//  track->SetP(p,kTRUE);
+
+  AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
+  const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+   // recalculate primary vertex for peri. and pp
+   AliVertexerTracks vertexer(fEvent->GetMagneticField());
+   vertexer.SetITSMode();
+   vertexer.SetMinClusters(4);
+   Int_t skipped[2];
+   skipped[0] = track->GetID();
+   vertexer.SetSkipTracks(1,skipped);
+   vertexer.SetConstraintOn();
+ //----diamond constraint-----------------------------
+   Float_t diamondcovxy[3];
+   esdevent->GetDiamondCovXY(diamondcovxy);
+   Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
+   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
+   AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
+   vertexer.SetVtxStart(diamond);
+   delete diamond; diamond=NULL;
+ //----------------------------------------------------
+   vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
+
+   // Getting the DCA
+   // Propagation always done on a working copy to not disturb the track params of the original track
+   AliESDtrack *esdtrack = NULL;
+   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
+     // Case ESD track: take copy constructor
+     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
+     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
+   } else {
+    // Case AOD track: take different constructor
+    esdtrack = new AliESDtrack(track);
+   }
+   if(esdtrack){
+    et.PropagateToDCA(tmpvtx, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
+    //et.PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadiusnew, dcan, covn);
+    esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcao, covo);
+    dcaxyn = dcan[0];
+    dcaxyo = dcao[0];
+    dcazn = dcan[1];
+    dcazo = dcao[1];
+    Double_t resfactor=1., resfactorz=1.;
+    if(sd0rpo) resfactor = sd0rpn/sd0rpo;
+    if(sd0zo) resfactorz = sd0zn/sd0zo;
+    if(covo[0]) dcaxySign = dcan[0]/(TMath::Sqrt(covo[0])*resfactor);
+    if(covo[0]) dcaxySigo = dcao[0]/TMath::Sqrt(covo[0]);
+    if(covo[2]) dcazSign = dcan[1]/(TMath::Sqrt(covo[2])*resfactorz);
+    if(covo[2]) dcazSigo = dcao[1]/TMath::Sqrt(covo[2]);
+    //if(dcaxyo) printf("pt= %lf resolusion xy= %lf dcaxyo= %lf dcaxyn=  %lf ratio= %lf\n",pt1mc,resfactor,dcaxyo,dcaxyn,dcaxyn/dcaxyo); 
+
+    if(!covo[0]){
+     dcaxySigo = -49.;
+     dcaxySign = -49.;
+    }
+    if(!covo[2]){
+     dcazSigo = -49.;
+     dcazSign = -49.;
+    }
+   }
+   delete esdtrack;
+   //delete vtxESDSkip;
+
+}
+
+//______________________________________________________
+Double_t AliHFEsmearDCA::EvalGraph(const TGraph *graph,Double_t x) const {
+  //
+  // Evaluates a TGraph without linear extrapolation. Instead the last
+  // valid point of the graph is used when out of range.
+  // The function assumes an ascending order of X.
+  //
+
+  // TODO: find a pretty solution for this:
+  Int_t    n   =graph->GetN();
+  Double_t xmin=graph->GetX()[0  ];
+  Double_t xmax=graph->GetX()[n-1];
+  if (x<xmin) return graph->Eval(xmin);
+  if (x>xmax) return graph->Eval(xmax);
+  return graph->Eval(x);
+}
diff --git a/PWGHF/hfe/AliHFEsmearDCA.h b/PWGHF/hfe/AliHFEsmearDCA.h
new file mode 100644 (file)
index 0000000..7b29779
--- /dev/null
@@ -0,0 +1,76 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Get improved dca info. by ITS upgrade implemented by the ALICE Heavy Flavour Electron Group
+//
+#ifndef ALIHFESMEARDCA_H
+#define ALIHFESMEARDCA_H
+
+class TList;
+class TGraph;
+
+class AliVEvent;
+class AliVParticle;
+class AliVTrack;
+class AliMCEvent;
+
+class AliHFEsmearDCA : public TObject{
+  public:
+  AliHFEsmearDCA(const Char_t */*name*/, const char *resfileCurURI, const char *resfileUpgURI, const Char_t */*title*/);
+    AliHFEsmearDCA(const AliHFEsmearDCA &c);
+    AliHFEsmearDCA &operator=(const AliHFEsmearDCA &);
+    virtual ~AliHFEsmearDCA();
+    
+    virtual void SetRecEventInfo(const TObject *event);
+
+    void SetMCEvent(AliMCEvent* const mcEvent){fMCEvent=mcEvent;};  // set stack pointer
+    void GetImproveITSImpactParameters(AliVTrack *track, Double_t &dcaxyn, Double_t &dcaxyo, Double_t &dcaxysign, Double_t &dcaxysigo, Double_t &dcazn, Double_t &dcazo, Double_t &dcazsign, Double_t &dcazsigo); // to check improvement by the ITS upgrade
+    Double_t EvalGraph(const TGraph *graph,Double_t x) const;
+    
+
+  private:
+    AliVEvent *fEvent;                //! working event
+    AliMCEvent *fMCEvent;             //! MCEvent pointer
+
+    TGraph *fD0ZResPCur  ; // old pt dep. d0 res. in z for protons
+    TGraph *fD0ZResKCur  ; // old pt dep. d0 res. in z for kaons
+    TGraph *fD0ZResPiCur ; // old pt dep. d0 res. in z for pions
+    TGraph *fD0ZResECur  ; // old pt dep. d0 res. in z for electrons 
+    TGraph *fD0RPResPCur ; // old pt dep. d0 res. in rphi for protons
+    TGraph *fD0RPResKCur ; // old pt dep. d0 res. in rphi for kaons
+    TGraph *fD0RPResPiCur; // old pt dep. d0 res. in rphi for pions
+    TGraph *fD0RPResECur ; // old pt dep. d0 res. in rphi for electrons
+    TGraph *fPt1ResPCur  ; // old pt dep. 1/pt res. for protons
+    TGraph *fPt1ResKCur  ; // old pt dep. 1/pt res. for kaons
+    TGraph *fPt1ResPiCur ; // old pt dep. 1/pt res. for pions
+    TGraph *fPt1ResECur  ; // old pt dep. 1/pt res. for electrons
+    TGraph *fD0ZResPUpg  ; // new pt dep. d0 res. in z for protons
+    TGraph *fD0ZResKUpg  ; // new pt dep. d0 res. in z for kaons
+    TGraph *fD0ZResPiUpg ; // new pt dep. d0 res. in z for pions
+    TGraph *fD0ZResEUpg  ; // new pt dep. d0 res. in z for electrons
+    TGraph *fD0RPResPUpg ; // new pt dep. d0 res. in rphi for protons
+    TGraph *fD0RPResKUpg ; // new pt dep. d0 res. in rphi for kaons
+    TGraph *fD0RPResPiUpg; // new pt dep. d0 res. in rphi for pions
+    TGraph *fD0RPResEUpg ; // new pt dep. d0 res. in rphi for electrons
+    TGraph *fPt1ResPUpg  ; // new pt dep. 1/pt res. for protons
+    TGraph *fPt1ResKUpg  ; // new pt dep. 1/pt res. for kaons
+    TGraph *fPt1ResPiUpg ; // new pt dep. 1/pt res. for pions
+    TGraph *fPt1ResEUpg  ; // new pt dep. 1/pt res. for electrons
+  
+  ClassDef(AliHFEsmearDCA, 1)      // Additional cuts implemented by the ALICE HFE group
+};
+
+#endif
index 4f032a4..3349817 100644 (file)
@@ -390,7 +390,10 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
         ptcorrelation->Draw("colz");
       }
     }
-    if(fWriteToFile) ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
+    if(fWriteToFile){ 
+      ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
+      ccorrelation->SaveAs("correlationMatrix.eps");
+    }
   }
 
   TFile *file = TFile::Open("tests.root","recreate");