updated dndeta analysis for common plots for MB&UE working group (Chiara)
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Oct 2010 18:27:12 +0000 (18:27 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Oct 2010 18:27:12 +0000 (18:27 +0000)
12 files changed:
PWG0/AliCorrection.cxx
PWG0/AliCorrectionMatrix3D.cxx
PWG0/AliCorrectionMatrix3D.h
PWG0/AliPWG0Helper.cxx
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
PWG0/dNdEta/AlidNdEtaCorrectionTask.h
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/AlidNdEtaTask.h
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/dNdEtaAnalysis.h

index bca4767..f09ec2c 100644 (file)
@@ -103,21 +103,30 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Hel
   //Float_t binLimitsEta[] = {-2.8,-2.4,-2,-1.6,-1.2,-0.8,-0.5,0,0.5,0.8,1.2,1.6,2.0,2.4,2.8}; 
 
   //Float_t binLimitsEta[] = {-3,-2.6,-2.2,-1.8,-1.4,-1,-0.6,-0.2,0.2,0.6,1,1.4,1.8,2.2,2.6,3.0};
-  Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
+  //1. standard
+  //  Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
+  // 2. MB Working Group definition
+  Float_t binLimitsEta[] = {-2, -1.8, -1.6, -1.4, -1.2, -1., -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0};
 //  Float_t binLimitsEta[] = { -7.0, -6.9, -6.8, -6.7, -6.6, -6.5, -6.4, -6.3, -6.2, -6.1, -6.0, -5.9, -5.8, -5.7, -5.6, -5.5, -5.4, -5.3, -5.2, -5.1, -5.0, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4.0, -3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3.0, -2.9, -2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, -0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0 };
 
   fEventCorr = new AliCorrectionMatrix2D("EventCorrection", Form("%s EventCorrection", fTitle.Data()), 22, binLimitsVtx, nBinsN, binLimitsN);
 
   if (analysisMode & AliPWG0Helper::kSPD)
   {
-    TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 40, binLimitsEta, nBinsN2, binLimitsN2);
+    // 1. standard
+    //TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 40, binLimitsEta, nBinsN2, binLimitsN2);
+    // 2. MB Working Group definition
+    TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 20, binLimitsEta, nBinsN2, binLimitsN2);
     fTrackCorr = new AliCorrectionMatrix3D("TrackCorrection", Form("%s TrackCorrection", fTitle.Data()), dummyBinning);
     fTrackCorr->SetAxisTitles("vtx-z (cm)", "#eta", title3);
     delete dummyBinning;
   }
   else
   {
-    TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 40, binLimitsEta , nBinsPt, binLimitsPt);
+    // 1. standard
+    //TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 40, binLimitsEta , nBinsPt, binLimitsPt);
+    // 2. MB Working Group definition
+    TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",22, binLimitsVtx, 20, binLimitsEta , nBinsPt, binLimitsPt);
     fTrackCorr = new AliCorrectionMatrix3D("TrackCorrection", Form("%s TrackCorrection", fTitle.Data()), dummyBinning);
     fTrackCorr->SetAxisTitles("vtx-z (cm)", "#eta", "p_{T} (GeV/c)");
     delete dummyBinning;
index a54424e..f4ee9fc 100644 (file)
@@ -282,17 +282,17 @@ TH1* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Option_t* opt, Float_t aMin
 }
 
 //____________________________________________________________________
-void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az)
+void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az, Double_t weight)
 {
   // add value to measured histogram
-  GetMeasuredHistogram()->Fill(ax, ay, az);
+  GetMeasuredHistogram()->Fill(ax, ay, az, weight);
 }
 
 //____________________________________________________________________
-void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az)
+void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az, Double_t weight)
 {
   // add value to generated histogram
-  GetGeneratedHistogram()->Fill(ax, ay, az);
+  GetGeneratedHistogram()->Fill(ax, ay, az, weight);
 }
 
 //____________________________________________________________________
index e27167c..4e54a21 100644 (file)
@@ -50,8 +50,8 @@ public:
   TH2* Get2DCorrectionHistogram(Option_t* opt, Float_t aMin, Float_t aMax)     {return Get2DCorrection(opt,aMin,aMax)->GetCorrectionHistogram();}
   TH1* Get1DCorrectionHistogram(Option_t* opt, Float_t aMins1=0, Float_t aMax1=0, Float_t aMins2=0, Float_t aMax2=0);
 
-  void FillMeas(Float_t ax, Float_t ay, Float_t az);
-  void FillGene(Float_t ax, Float_t ay, Float_t az);
+  void FillMeas(Float_t ax, Float_t ay, Float_t az, Double_t weight = 1.);
+  void FillGene(Float_t ax, Float_t ay, Float_t az, Double_t weight = 1.);
 
   Float_t GetCorrection(Float_t ax, Float_t ay, Float_t az) const;
 
index a629c07..7ca5f73 100644 (file)
@@ -90,7 +90,7 @@ const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMo
     vertex = aEsd->GetPrimaryVertexTracks();
     if (debug)
       Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks");
-    if (vertex && vertex->GetNContributors() <= 0)
+    if (!vertex || vertex->GetNContributors() <= 0)
     {
       if (debug)
         Printf("AliPWG0Helper::GetVertex: Vertex from tracks has no contributors. Falling back to SPD vertex.");
index 65c6c88..2f52464 100644 (file)
@@ -330,11 +330,11 @@ void AlidNdEtaCorrection::FillMCParticle(Float_t vtx, Float_t eta, Float_t pt, B
 }
 
 //____________________________________________________________________
-void AlidNdEtaCorrection::FillTrackedParticle(Float_t vtx, Float_t eta, Float_t pt)
+void AlidNdEtaCorrection::FillTrackedParticle(Float_t vtx, Float_t eta, Float_t pt, Double_t weight)
 {
   // fills a tracked particle in the corrections
 
-  fTrack2ParticleCorrection->GetTrackCorrection()->FillMeas(vtx, eta, pt);
+  fTrack2ParticleCorrection->GetTrackCorrection()->FillMeas(vtx, eta, pt, weight);
 }
 
 //____________________________________________________________________
index c94978b..980a5a7 100644 (file)
@@ -43,7 +43,7 @@ public:
   ~AlidNdEtaCorrection();
 
   void FillMCParticle(Float_t vtx, Float_t eta, Float_t pt, Bool_t trigger, Bool_t vertex, Int_t processType);
-  void FillTrackedParticle(Float_t vtx, Float_t eta, Float_t pt);
+  void FillTrackedParticle(Float_t vtx, Float_t eta, Float_t pt, Double_t weight=1.);
   void FillEvent(Float_t vtx, Float_t n, Bool_t trigger, Bool_t vertex, Int_t processType);
 
   void Finish();
index 0593bec..300eb8d 100644 (file)
@@ -76,7 +76,23 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
   fMultAll(0),
   fMultTr(0),
   fMultVtx(0),
-  fEventStats(0)
+  fEventStats(0),
+  fEsdTrackCutsCheck(0),
+  fEtaCorrelationAllESD(0),
+  fpTCorrelation(0),
+  fpTCorrelationShift(0),
+  fpTCorrelationAllESD(0),
+  fpTCorrelationShiftAllESD(0),
+  fPtMin(0.15),
+  fPtMC(0),
+  fEtaMC(0),
+  fPtESD(0),
+  fEtaESD(0),
+  fVtxMC(0),
+  fNumberEventMC(0),
+  fNumberEvent(0),
+  fEventNumber(-1),
+  fWeightSecondaries(kFALSE)
 {
   //
   // Constructor. Initialization of pointers
@@ -131,7 +147,23 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   fMultAll(0),
   fMultTr(0),
   fMultVtx(0),
-  fEventStats(0)
+  fEventStats(0),
+  fEsdTrackCutsCheck(0),
+  fEtaCorrelationAllESD(0),
+  fpTCorrelation(0),
+  fpTCorrelationShift(0),
+  fpTCorrelationAllESD(0),
+  fpTCorrelationShiftAllESD(0),
+  fPtMin(0.15),
+  fPtMC(0),
+  fEtaMC(0),
+  fPtESD(0),
+  fEtaESD(0),
+  fVtxMC(0),
+  fNumberEventMC(0),
+  fNumberEvent(0),
+  fEventNumber(-1),
+  fWeightSecondaries(kFALSE)
 {
   //
   // Constructor. Initialization of pointers
@@ -199,6 +231,7 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
 {
   // create result objects and add to output list
 
+  AliDebug(2,Form("*********************************** fOption = %s", fOption.Data()));
   if (fOption.Contains("only-positive"))
   {
     Printf("INFO: Processing only positive particles.");
@@ -243,6 +276,9 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   {
     fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim"));
     fEsdTrackCutsSec  = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec"));
+    fEsdTrackCutsCheck  = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsCheck"));
+    fEsdTrackCutsCheck->SetPtRange(0.15);
+    fEsdTrackCutsCheck->SetEtaRange(-1.,1.);
     fOutput->Add(fEsdTrackCutsPrim);
     fOutput->Add(fEsdTrackCutsSec);
   }
@@ -288,6 +324,8 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
 
   fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
   fOutput->Add(fEtaCorrelation);
+  fEtaCorrelationAllESD = new TH2F("fEtaCorrelationAllESD", "fEtaCorrelationAllESD;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
+  fOutput->Add(fEtaCorrelationAllESD);
   fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
   fOutput->Add(fEtaCorrelationShift);
   fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
@@ -298,6 +336,11 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
   fpTResolution = new TH2F("fpTResolution", ";MC p_{T} (GeV/c);(MC p_{T} - ESD p_{T}) / MC p_{T}", 160, 0, 20, 201, -0.2, 0.2);
   fOutput->Add(fpTResolution);
 
+  fpTCorrelationAllESD = new TH2F("fpTCorrelationAllESD", "fpTCorrelationAllESD;MC p_{T} (GeV/c);ESD p_{T}", 160, 0, 20, 160, 0, 20);
+  fOutput->Add(fpTCorrelationAllESD);
+  fpTCorrelationShiftAllESD = new TH2F("fpTCorrelationShiftAllESD", "fpTCorrelationShiftAllESD;MC p_{T} (GeV/c);MC p_{T} - ESD p_{T}", 160, 0, 20, 100, -1, 1);
+  fOutput->Add(fpTCorrelationShiftAllESD);
+
   fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
   fOutput->Add(fMultAll);
   fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
@@ -336,12 +379,30 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
     // TODO like this we send an empty object back...
     fOutput->Add(fEsdTrackCuts->Clone());
   }
+  fPtMC = new TH1F("fPtMC", "Pt from MC for selected tracks;MC p_{T} (GeV/c)", 160, 0, 20);
+  fOutput->Add(fPtMC); 
+  fEtaMC = new TH1F("fEtaMC", "Eta from MC for selected tracks;MC #eta", 120, -3, 3);
+  fOutput->Add(fEtaMC);
+  fPtESD = new TH1F("fPtESD", "Pt from ESD for selected tracks;ESD p_{T} (GeV/c)", 160, 0, 20);
+  fOutput->Add(fPtESD);
+  fEtaESD = new TH1F("fEtaESD", "Eta from ESD for selected tracks;ESD #eta", 120, -3, 3);
+  fOutput->Add(fEtaESD);
+
+  fVtxMC = new TH1F("fVtxMC", "Vtx,z from MC for all events;MC vtx_z (cm)", 100, -30, 30);
+  fOutput->Add(fVtxMC); 
+
+  fNumberEventMC = new TH1F("fNumberEventMC","Number of event accepted at MC level",600000,-0.5,600000-0.5);
+  fOutput->Add(fNumberEventMC);
+
+  fNumberEvent = new TH1F("fNumberEvent","Number of event accepted at Reco level",600000,-0.5,600000-0.5);
+  fOutput->Add(fNumberEvent);
 }
 
 void AlidNdEtaCorrectionTask::Exec(Option_t*)
 {
   // process the event
 
+  fEventNumber++;
   // Check prerequisites
   if (!fESD)
   {
@@ -423,6 +484,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   AliGenEventHeader* genHeader = header->GenEventHeader();
   TArrayF vtxMC(3);
   genHeader->PrimaryVertex(vtxMC);
+  fVtxMC->Fill(vtxMC[2]);
+  AliDebug(2,Form("MC vtx: x = %.6f, y = %.6f, z = %.6f",vtxMC[0],vtxMC[1],vtxMC[2]));
 
   // get the ESD vertex
   Bool_t eventVertex = kFALSE;
@@ -541,6 +604,44 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     
     if (vtxESD)
     {
+      // control histograms on pT
+      Int_t nfromstack = stack->GetNtrack();
+      AliDebug(3,Form(" from stack we have %d tracks\n",nfromstack));
+      for (Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
+       AliESDtrack* esdTrackcheck = dynamic_cast<AliESDtrack*> (fESD->GetTrack(itrack));
+       if (!esdTrackcheck){
+         AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", itrack));
+         continue;
+       }
+       if (fOnlyPrimaries){
+         Int_t label = TMath::Abs(esdTrackcheck->GetLabel());
+         AliDebug(4,Form("label = %d\n",label));
+         if (label == 0 || label > nfromstack) continue;
+         if (stack->IsPhysicalPrimary(label) == kFALSE) continue;
+       }
+        
+       Int_t label = TMath::Abs(esdTrackcheck->GetLabel());
+       if (label == 0 || label > nfromstack) continue;
+       if (!stack->Particle(label)){
+         AliDebug(4,Form("WARNING: No particle for %d", label));
+       }
+       else{   
+         if (!fEsdTrackCuts->AcceptTrack(esdTrackcheck)){
+           TParticle* particle = stack->Particle(label);
+           if (fEsdTrackCutsCheck->AcceptTrack(esdTrackcheck)){
+             //if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > fPtMin){
+             Float_t ptMC = particle->Pt();
+             Float_t etaMC = particle->Eta();
+             Float_t ptESD = esdTrackcheck->Pt();
+             Float_t etaESD = esdTrackcheck->Eta();
+             fEtaCorrelationAllESD->Fill(etaMC,etaESD);
+             fpTCorrelationAllESD->Fill(ptMC,ptESD);
+             fpTCorrelationShiftAllESD->Fill(ptMC,ptMC-ptESD);
+           }
+         }
+       }
+      } // end loop over all ESDs
+
       // get multiplicity from ESD tracks
       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
       Int_t nGoodTracks = list->GetEntries();
@@ -563,7 +664,8 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           continue;
         }
         
-        if (esdTrack->Pt() < 0.15)
+       AliDebug(3,Form("particle %d: pt = %.6f, eta = %.6f",i,esdTrack->Pt(), esdTrack->Eta())); 
+        if (esdTrack->Pt() < fPtMin)
           continue;
         
         if (fOnlyPrimaries)
@@ -576,10 +678,28 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
             continue;
         }
         
-        // INEL>0 trigger
-        if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
+       Int_t label = TMath::Abs(esdTrack->GetLabel());
+       if (!stack->Particle(label)){
+         AliDebug(3,Form("WARNING: No particle for %d", label));
+       }
+       else{
+         TParticle* particle = stack->Particle(label);
+         Float_t ptMC = particle->Pt();
+         Float_t etaMC = particle->Eta();
+         fPtMC->Fill(ptMC);
+         fEtaMC->Fill(etaMC);
+         fPtESD->Fill(esdTrack->Pt());
+         fEtaESD->Fill(esdTrack->Eta());
+       }
+
+        // 2 Options for INEL>0 trigger - choose one
+       // 1. HL
+        //if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
+       // foundInEta10 = kTRUE;
+       // 2. MB Working Group definition
+       if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin)
           foundInEta10 = kTRUE;
-  
+
         etaArr[inputCount] = esdTrack->Eta();
         if (fSymmetrize)
           etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
@@ -656,6 +776,10 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     
     if (!foundInEta10)
       eventTriggered = kFALSE;
+    else{
+     //Printf("The event triggered. Its number in file is %d",fESD->GetEventNumberInFile());
+      fNumberEvent->Fill(fESD->GetEventNumberInFile());
+    }
   }
   else
     return;
@@ -724,7 +848,11 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
     if (SignOK(particle->GetPDG()) == kFALSE)
       continue;
     
-    if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
+    // for INEL > 0, MB Working Group definition use the second option
+    // 1. standard
+    //if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
+    // 2. MB Working Group definition
+    if (fPIDParticles && TMath::Abs(particle->Eta()) < 0.8 && particle->Pt() > fPtMin)
       fPIDParticles->Fill(particle->GetPdgCode());
 
     Float_t eta = particle->Eta();
@@ -839,7 +967,11 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       continue;
     }
 
-    if (TMath::Abs(particle->Eta()) < 1.0)
+    // for INEL > 0, MB Working Group definition use the second option
+    // 1. standard
+    //if (TMath::Abs(particle->Eta()) < 1.0)
+    // 2. INEL > 0 MB Working Group definition
+    if (TMath::Abs(particle->Eta()) < 0.8 && particle->Pt()>fPtMin)
       fPIDTracks->Fill(particle->GetPdgCode());
     
     // find particle that is filled in the correction map
@@ -878,9 +1010,17 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       else
         fEtaResolution->Fill(particle->Eta() - etaArr[i]);
 
-      if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
-        if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
+      if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS){
+       // for INEL > 0, MB Working Group definition use the second option
+       // 1. standard
+        //if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
+       // 2. INEL > 0 MB WOrking Group definition
+       if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > 0){
           fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
+         fpTCorrelation->Fill(particle->Pt(),thirdDimArr[i]);
+         fpTCorrelationShift->Fill(particle->Pt(),particle->Pt()-thirdDimArr[i]);
+       }
+      }
 
       Float_t eta = -999;
       Float_t thirdDim = -1;
@@ -942,7 +1082,13 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
 
       if (fillTrack)
       {
-        fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim);
+        Double_t weight = 1.;
+       if (fWeightSecondaries){
+         if (!firstIsPrim){
+           weight = GetSecondaryCorrection(thirdDim);
+         }
+       }        
+       fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim, weight);
         fTemp2->Fill(vtxMC[2], eta);
       }
       
@@ -954,7 +1100,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
           eta2 = TMath::Abs(eta2);
         
         fEtaProfile->Fill(eta2, eta2 - etaArr[i]);
-        fEtaCorrelation->Fill(etaArr[i], eta2);
+        fEtaCorrelation->Fill(eta2, etaArr[i]);
         fEtaCorrelationShift->Fill(eta2, eta2 - etaArr[i]);
       }
 
@@ -1213,6 +1359,9 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
   if (fEtaCorrelation)
     fEtaCorrelation->Write();
+  fEtaCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationAllESD"));
+  if (fEtaCorrelationAllESD)
+    fEtaCorrelationAllESD->Write();
   fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
   if (fEtaCorrelationShift)
   {
@@ -1225,6 +1374,15 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution"));
   if (fEtaResolution)
     fEtaResolution->Write();
+  fpTCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelation"));
+  if (fpTCorrelation)
+    fpTCorrelation->Write();
+  fpTCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShift"));
+  if (fpTCorrelationShift)
+  {
+    fpTCorrelationShift->FitSlicesY();
+    fpTCorrelationShift->Write();
+  }
   fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution"));
   if (fpTResolution)
   {
@@ -1232,6 +1390,15 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
     fpTResolution->Write();
   }
 
+  fpTCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationAllESD"));
+  if (fpTCorrelationAllESD)
+    fpTCorrelationAllESD->Write();
+  fpTCorrelationShiftAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShiftAllESD"));
+  if (fpTCorrelationShiftAllESD)
+  {
+    fpTCorrelationShiftAllESD->FitSlicesY();
+    fpTCorrelationShiftAllESD->Write();
+  }
   fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
   if (fMultAll)
     fMultAll->Write();
@@ -1263,9 +1430,32 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   if (fPIDTracks)
     fPIDTracks->Write();
 
+  fPtMC = dynamic_cast<TH1F*> (fOutput->FindObject("fPtMC"));
+  if (fPtMC)
+    fPtMC->Write();
+  fEtaMC = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaMC"));
+  if (fEtaMC)
+    fEtaMC->Write();
+  fPtESD = dynamic_cast<TH1F*> (fOutput->FindObject("fPtESD"));
+  if (fPtESD)
+    fPtESD->Write();
+  fEtaESD = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaESD"));
+  if (fEtaESD)
+    fEtaESD->Write();
+  fVtxMC = dynamic_cast<TH1F*> (fOutput->FindObject("fVtxMC"));
+  if (fVtxMC)
+    fVtxMC->Write();
+
+  fNumberEventMC = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEventMC"));
+  if (fNumberEventMC)
+    fNumberEventMC->Write();
+  fNumberEvent = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEvent"));
+  if (fNumberEvent)
+    fNumberEvent->Write();
+
  //fdNdEtaCorrection->DrawHistograms();
 
-  Printf("Writting result to %s", fileName.Data());
+  Printf("Writing result to %s", fileName.Data());
 
   if (fPIDParticles && fPIDTracks)
   {
@@ -1315,3 +1505,27 @@ Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
   return kTRUE;
 }
 
+Double_t AlidNdEtaCorrectionTask::GetSecondaryCorrection(Double_t pt){
+
+       // getting the data driven correction factor to correct for 
+       // the underestimate of secondaries in Pythia
+       // (A. Dainese + J. Otwinowski
+
+       if (pt <= 0.17) return 1.0;
+       if (pt <= 0.4) return GetLinearInterpolationValue(0.17,1.0,0.4,1.07, pt);
+       if (pt <= 0.6) return GetLinearInterpolationValue(0.4,1.07,0.6,1.25, pt);
+       if (pt <= 1.2) return GetLinearInterpolationValue(0.6,1.25,1.2,1.5,  pt);
+       return 1.5;
+
+}
+
+Double_t AlidNdEtaCorrectionTask::GetLinearInterpolationValue(Double_t const x1, Double_t const y1, Double_t const x2, Double_t const y2, const Double_t pt)
+{
+
+       //
+       // linear interpolation to be used to get the secondary correction (see AlidNdEtaCorrectionTask::GetSecondaryCorrection)
+       //
+
+       return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
+}
+
index 8a4946b..cdd1a94 100644 (file)
@@ -6,8 +6,8 @@
 #include "AliAnalysisTask.h"
 #include <TString.h>
 #include "AliPWG0Helper.h"
+#include "AliESDtrackCuts.h"
 
-class AliESDtrackCuts;
 class dNdEtaAnalysis;
 class AlidNdEtaCorrection;
 class TH1;
@@ -41,6 +41,10 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     void SetSkipParticles(Bool_t flag = kTRUE) { fSystSkipParticles = flag; }
 
     void SetOption(const char* opt) { fOption = opt; }
+    void SetPtMin(Float_t ptMin) {fPtMin = ptMin;}
+    void SetWeightSecondaries(Bool_t flag) { fWeightSecondaries = flag;}
+    Double_t GetSecondaryCorrection(Double_t pt);
+    Double_t GetLinearInterpolationValue(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t pt);
 
  protected:
     Bool_t SignOK(TParticlePDG* particle);
@@ -106,12 +110,28 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     AlidNdEtaCorrection* fdNdEtaCorrectionSpecial[4];   //! correction maps used for systematic studies, may contain:
                                                         // for specific process type (ND, SD, DD), enable with option: process-types
                                                         // for particle species (pi, K, p, rest), enable with: particle-species
+    AliESDtrackCuts*  fEsdTrackCutsCheck;        //! Object containing the parameters of the esd track cuts
+    TH2F* fEtaCorrelationAllESD;                       //! ESD eta vs MC eta
+    TH2F* fpTCorrelation;                         //! ESD pT vs MC pT in |eta| < 0.9
+    TH2F* fpTCorrelationShift;                    //! (MC pT - ESD pT) vs MC pT in |eta| < 0.9
+    TH2F* fpTCorrelationAllESD;                         //! ESD pT vs MC pT in |eta| < 0.9
+    TH2F* fpTCorrelationShiftAllESD;                    //! (MC pT - ESD pT) vs MC pT in |eta| < 0.9
+    Float_t fPtMin;    // ptMin for kOneTrack
+    TH1F* fPtMC;       //! pT histogram for MC information for selected tracks
+    TH1F* fEtaMC;      //! eta histogram for MC information for selected tracks
+    TH1F* fPtESD;      //! pT histogram for ESD information for selected tracks
+    TH1F* fEtaESD;     //! eta histogram for ESD information for selected tracks
+    TH1F* fVtxMC;      //! vtx_z histogram for MC information for all events
+    TH1F* fNumberEventMC;      //! number of accepted event histogram for MC information for all events
+    TH1F* fNumberEvent;      //! number of accepted event histogram for reco information for all events
+    Int_t fEventNumber;      // number of the event - useful when running on one file, on one worker
+    Bool_t fWeightSecondaries; // is true if calculating corrections to be applied to real data (secondaries correction should be on)
 
  private:
     AlidNdEtaCorrectionTask(const AlidNdEtaCorrectionTask&);
     AlidNdEtaCorrectionTask& operator=(const AlidNdEtaCorrectionTask&);
 
-  ClassDef(AlidNdEtaCorrectionTask, 1);
+  ClassDef(AlidNdEtaCorrectionTask, 2);
 };
 
 #endif
index fe54c9e..ca0419b 100644 (file)
@@ -9,6 +9,7 @@
 #include <TChain.h>
 #include <TFile.h>
 #include <TH1F.h>
+#include <TH1D.h>
 #include <TH2F.h>
 #include <TH3F.h>
 #include <TParticle.h>
@@ -94,7 +95,18 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fTrackletsVsClusters(0),
   fTrackletsVsUnassigned(0),
   fStats(0),
-  fStats2(0)
+  fStats2(0),
+  fPtMin(0.15),
+  fEta(0x0),
+  fEtaMC(0x0),
+  fHistEvents(0),
+  fHistEventsMC(0),
+  fTrigEffNum(0),
+  fTrigEffDen(0),
+  fVtxEffNum(0),
+  fVtxEffDen(0),
+  fVtxTrigEffNum(0),
+  fVtxTrigEffDen(0)
 {
   //
   // Constructor. Initialization of pointers
@@ -327,6 +339,31 @@ void AlidNdEtaTask::UserCreateOutputObjects()
     fEsdTrackCuts->SetName("fEsdTrackCuts");
     fOutput->Add(fEsdTrackCuts);
   }
+
+  fEta = new TH1D("fEta", "Eta;#eta;count", 80, -4, 4);
+  fOutput->Add(fEta);
+  
+  fEtaMC = new TH1D("fEtaMC", "Eta, MC;#eta;count", 80, -4, 4);
+  fOutput->Add(fEtaMC);
+
+  fHistEvents = new TH1D("fHistEvents", "N. of Events;accepted;count", 2, 0, 2);
+  fOutput->Add(fHistEvents);
+  
+  fHistEventsMC = new TH1D("fHistEventsMC", "N. of MC Events;accepted;count", 2, 0, 2);
+  fOutput->Add(fHistEventsMC);
+
+  fTrigEffNum = new TH1D("fTrigEffNum", "N. of triggered events", 100,0,100); 
+  fOutput->Add(fTrigEffNum);
+  fTrigEffDen = new TH1D("fTrigEffDen", "N. of true events", 100,0,100); 
+  fOutput->Add(fTrigEffDen);
+  fVtxEffNum = new TH1D("fVtxEffNum", "N. of events with vtx", 100,0,100); 
+  fOutput->Add(fVtxEffNum);
+  fVtxEffDen = new TH1D("fVtxEffDen", "N. of true events", 100,0,100); 
+  fOutput->Add(fVtxEffDen);
+  fVtxTrigEffNum = new TH1D("fVtxTrigEffNum", "N. of triggered events with vtx", 100,0,100); 
+  fOutput->Add(fVtxTrigEffNum);
+  fVtxTrigEffDen = new TH1D("fVtxTrigEffDen", "N. of triggered true events", 100,0,100); 
+  fOutput->Add(fVtxTrigEffDen);
   
   //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
 }
@@ -335,6 +372,7 @@ Bool_t AlidNdEtaTask::IsEventInBinZero()
 {
   // checks if the event goes to the 0 bin
   
+  Bool_t isZeroBin = kTRUE;
   fESD = (AliESDEvent*) fInputEvent;
   
   AliInputEventHandler* inputHandler = static_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
@@ -361,9 +399,29 @@ Bool_t AlidNdEtaTask::IsEventInBinZero()
   if (!triggerAnalysis->IsTriggerFired(fESD, fTrigger))
     return kFALSE;
   
-  // TODO 0 bin check
+  // 0 bin check - from Michele
   
-  return kTRUE;
+  const AliMultiplicity* mult = fESD->GetMultiplicity();
+  if (!mult){
+    Printf("AlidNdEtaTask::IsBinZero: Can't get multiplicity object from ESDs");
+    return kFALSE;
+  }
+  Int_t ntracklet = mult->GetNumberOfTracklets();
+  const AliESDVertex * vtxESD = fESD->GetPrimaryVertexSPD();
+  if(vtxESD) {
+         // If there is a vertex from vertexer z with delta phi > 0.02 we
+         // don't consider it rec (we keep the event in bin0). If quality
+         // is good eneough we check the number of tracklets
+         // if the vertex is more than 15 cm away, this is autamatically bin0
+         if( TMath::Abs(vtxESD->GetZ()) <= 15 ) {
+                 if (vtxESD->IsFromVertexerZ()) {
+                         if (vtxESD->GetDispersion()<=0.02 ) {
+                                 if(ntracklet>0) isZeroBin = kFALSE;
+                         }
+                 } else if(ntracklet>0) isZeroBin = kFALSE; // if the event is not from Vz we chek the n of tracklets
+         } 
+  }
+  return isZeroBin;
 }
 
 void AlidNdEtaTask::UserExec(Option_t*)
@@ -546,6 +604,7 @@ isManager()->GetInputEventHandler());
         
       const AliMultiplicity* mult = fESD->GetMultiplicity();
       if (!mult)
+       Printf("Returning, no Multiplicity found");
         return;
       
       if (mult->GetNumberOfTracklets() == 0)
@@ -821,33 +880,45 @@ isManager()->GetInputEventHandler());
           if (eventTriggered && vtxESD)
             fRawPt->Fill(pT);
   
-          if (esdTrack->Pt() < 0.15)
+          if (esdTrack->Pt() < fPtMin) // even if the pt cut is already in the esdTrackCuts....
             continue;
           
-          // INEL>0 trigger
-          if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
-            foundInEta10 = kTRUE;
-  
           if (fOnlyPrimaries)
           {
-            if (label == 0)
-              continue;
-            
-            if (stack->IsPhysicalPrimary(label) == kFALSE)
-              continue;
+          if (label == 0)
+            continue;
+          
+          if (stack->IsPhysicalPrimary(label) == kFALSE)
+            continue;
           }
+
+         // 2 types of INEL>0 trigger - choose one
+
+         // 1. HL 
+          // INEL>0 trigger
+          // if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
+          // foundInEta10 = kTRUE;
+
+          // 2. MB working group
+         if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin){  // this should be in the trigger selection as well, so one particle should always be found
+           foundInEta10 = kTRUE;
+         }
   
           if (fUseMCKine)
           {
-            if (label > 0)
+           if (label > 0)
             {
               TParticle* particle = stack->Particle(label);
               eta = particle->Eta();
               pT = particle->Pt();
+             // check when using INEL>0, MB Working Group definition
+             if (TMath::Abs(eta) >=0.8 || pT <= fPtMin){
+               AliDebug(2,Form("WARNING ************* USING KINE: eta = %f, pt = %f",eta,pT));
+             }
             }
             else
               Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
-          }
+         }
   
           if (fSymmetrize)
             eta = TMath::Abs(eta);
@@ -865,8 +936,20 @@ isManager()->GetInputEventHandler());
         delete list;
       }
       
-      if (!foundInEta10)
+      if (!foundInEta10){
         eventTriggered = kFALSE;
+       fHistEvents->Fill(0.1);
+       AliDebug(3,"Event rejected");
+      }
+      else{
+       if (eventTriggered){
+         fHistEvents->Fill(1.1);
+         AliDebug(3,Form("Event Accepted, with inputcount = %d",inputCount));
+       }
+       else{
+         AliDebug(3,"Event has foundInEta10 but was not triggered");
+       }
+      }
     }
     else
       return;
@@ -906,7 +989,7 @@ isManager()->GetInputEventHandler());
           
           if (TMath::Abs(eta) < 0.5)
             part05++;
-          if (TMath::Abs(eta) < 1.0)
+          if (TMath::Abs(eta) < 1.0) // in the INEL>0, MB WG definition, this is in any case equivalent to <0.8, due to the EsdTrackCuts
             part10++;
         }
         
@@ -1028,6 +1111,8 @@ isManager()->GetInputEventHandler());
     Int_t nAcceptedParticles = 0;
     Bool_t oneParticleEvent = kFALSE;
 
+    Int_t nMCparticlesinRange = 0; // number of true particles in my range of interest
+
     // count particles first, then fill
     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
     {
@@ -1046,14 +1131,48 @@ isManager()->GetInputEventHandler());
       //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
       Float_t eta = particle->Eta();
 
-      if (TMath::Abs(eta) < 1.0)
-        oneParticleEvent = kTRUE;
+      AliDebug(2,Form("particle %d: eta = %f, pt = %f",iMc,particle->Eta(),particle->Pt()));
+
+      // INEL>0: choose one
+      // 1.HL
+      // if (TMath::Abs(eta) < 1.0){
+      //   oneParticleEvent = kTRUE;
+      //   nMCparticlesinRange++;
+      //}
+      // 2.MB Working Group definition
+      if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
+       oneParticleEvent = kTRUE;
+       nMCparticlesinRange++;
+      }
 
       // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: these histograms are only gathered for comparison))
       if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
         nAcceptedParticles++;
     }
 
+    if (TMath::Abs(vtxMC[2]) < 10.){
+      // if MC vertex is inside 10
+      if (vtxESD){
+       fVtxEffNum->Fill(nMCparticlesinRange);
+      }
+      fVtxEffDen->Fill(nMCparticlesinRange);
+      if (eventTriggered){
+       if (vtxESD){
+          fVtxTrigEffNum->Fill(nMCparticlesinRange);
+       }
+       fVtxTrigEffDen->Fill(nMCparticlesinRange);
+       fTrigEffNum->Fill(nMCparticlesinRange);
+      }
+      fTrigEffDen->Fill(nMCparticlesinRange);
+    }
+
+    if (oneParticleEvent){
+           fHistEventsMC->Fill(1.1);
+    }
+    else{
+           fHistEventsMC->Fill(0.1);
+    }  
+
     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
     {
       if (!stack->IsPhysicalPrimary(iMc))
@@ -1095,6 +1214,7 @@ isManager()->GetInputEventHandler());
         fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
       
       if (oneParticleEvent)
+       AliDebug(3,Form("filling dNdEtaAnalysis object:: vtx = %f, eta = %f, pt = %f",vtxMC[2], eta, thirdDim));
         fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
 
       if (eventTriggered)
@@ -1109,6 +1229,9 @@ isManager()->GetInputEventHandler());
         //Float_t value = 1. / TMath::TwoPi() / particle->Pt() * particle->Energy() / particle->P();
         Float_t value = 1;
         fPartPt->Fill(particle->Pt(), value);
+       if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
+         fEtaMC->Fill(eta);
+       }
       }
     }
 
@@ -1139,6 +1262,8 @@ void AlidNdEtaTask::Terminate(Option_t *)
   if (!fOutput)
     Printf("ERROR: fOutput not available");
 
+  TH1D* dNdEta = new TH1D("dNdEta","#eta;counts",80, -4, 4);
+  TH1D* dNdEtaMC = new TH1D("dNdEtaMC","#eta,MC;counts",80, -4, 4);
   if (fOutput)
   {
     fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
@@ -1166,6 +1291,52 @@ void AlidNdEtaTask::Terminate(Option_t *)
     fStats2 = dynamic_cast<TH2F*> (fOutput->FindObject("fStats2"));
 
     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
+
+    fEta = dynamic_cast<TH1D*>(fOutput->FindObject("fEta"));
+    fEtaMC = dynamic_cast<TH1D*>(fOutput->FindObject("fEtaMC"));
+    fHistEvents = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEvents"));
+    fHistEventsMC = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEventsMC"));
+    fVtxEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffDen"));
+    fVtxEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffNum"));
+    fTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffDen"));
+    fTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffNum"));
+    fVtxTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffDen"));
+    fVtxTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffNum"));
+    Printf("Events selected = %f", fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1)));
+    Printf("Events selected frm MC = %f", fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1)));
+    Float_t nevents = fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1));
+    Float_t neventsMC = fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1));
+    for (Int_t ibin = 1; ibin <= fEta->GetNbinsX(); ibin++){
+      Float_t eta = fEta->GetBinContent(ibin)/nevents/fEta->GetBinWidth(ibin);
+      Float_t etaerr = fEta->GetBinError(ibin)/nevents/fEta->GetBinWidth(ibin);
+      Float_t etaMC = fEtaMC->GetBinContent(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
+      Float_t etaerrMC = fEtaMC->GetBinError(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
+      dNdEta->SetBinContent(ibin,eta);
+      dNdEta->SetBinError(ibin,etaerr);
+      dNdEtaMC->SetBinContent(ibin,etaMC);
+      dNdEtaMC->SetBinError(ibin,etaerrMC);
+    }
+    new TCanvas("eta", " eta ",50, 50, 550, 550) ;
+    fEta->Draw();
+    new TCanvas("etaMC", " etaMC ",50, 50, 550, 550) ;
+    fEtaMC->Draw();
+    new TCanvas("dNdEta", "#eta;dNdEta ",50, 50, 550, 550) ;
+    dNdEta->Draw();
+    new TCanvas("dNdEtaMC", "#eta,MC;dNdEta ",50, 50, 550, 550) ;
+    dNdEtaMC->Draw();
+    new TCanvas("Events", "Events;Events ",50, 50, 550, 550) ;
+    fHistEvents->Draw();
+    new TCanvas("Events, MC", "Events, MC;Events ",50, 50, 550, 550) ;
+    fHistEventsMC->Draw();
+    
+    TFile* outputFileCheck = new TFile("histogramsCheck.root", "RECREATE");
+    fEta->Write();
+    fEtaMC->Write();
+    dNdEta->Write();
+    dNdEtaMC->Write();
+    outputFileCheck->Write();
+    outputFileCheck->Close();
+
   }
 
   if (!fdNdEtaAnalysisESD)
@@ -1297,6 +1468,24 @@ void AlidNdEtaTask::Terminate(Option_t *)
     if (fVertexVsMult)
       fVertexVsMult->Write();
     
+    if (fVtxEffDen) fVtxEffDen->Write();
+    else Printf("fVtxEffDen not found");
+
+    if (fVtxEffNum) fVtxEffNum->Write();
+    else Printf("fVtxEffNum not found");
+
+    if (fVtxTrigEffNum) fVtxTrigEffNum->Write();
+    else Printf("fVtxTrigEffNum not found");
+
+    if (fVtxTrigEffDen) fVtxTrigEffDen->Write();
+    else Printf("fVtxTrigEffDen not found");
+
+    if (fTrigEffNum) fTrigEffNum->Write();
+    else Printf("fTrigEffNum not found");
+
+    if (fTrigEffDen) fTrigEffDen->Write();
+    else Printf("fTrigEffDen not found");
+
     fout->Write();
     fout->Close();
 
index b9585aa..9615bef 100644 (file)
@@ -15,6 +15,7 @@ class TH2F;
 class TH3F;
 class AliESDEvent;
 class AliTriggerAnalysis;
+class TH1D;
 
 class AlidNdEtaTask : public AliAnalysisTaskSE {
   public:
@@ -43,6 +44,7 @@ class AlidNdEtaTask : public AliAnalysisTaskSE {
     void SetDiffTreatment(AliPWG0Helper::DiffTreatment diffTreatment) { fDiffTreatment = diffTreatment; }
     
     void SetOption(const char* opt) { fOption = opt; }
+    void SetPtMin(Float_t ptMin) { fPtMin = ptMin;}
 
  protected:
     AliESDEvent *fESD;                         //! ESD object
@@ -100,12 +102,22 @@ class AlidNdEtaTask : public AliAnalysisTaskSE {
     TH2F* fTrackletsVsUnassigned; //! number of tracklets vs. number of unassigned clusters in L1 (only for SPD analysis)
     TH1F* fStats;                 //! further statistics : bin 1 = vertexer 3d, bin 2 = vertexer z, etc (see CreateOutputObjects)
     TH2F* fStats2;                //! V0 vs SPD statistics
-
+    Float_t fPtMin;               // pt min, to be used in kOneTrack case
+    TH1D* fEta;                   //! eta distribution from ESD
+    TH1D* fEtaMC;                 //! eta distribution from MC
+    TH1D* fHistEvents;            //! histo for n. of selected ESD events
+    TH1D* fHistEventsMC;          //! histo for n. of selected MC events
+    TH1D* fTrigEffNum;            //!
+    TH1D* fTrigEffDen;            //!
+    TH1D* fVtxEffNum;             //!
+    TH1D* fVtxEffDen;             //!
+    TH1D* fVtxTrigEffNum;         //!
+    TH1D* fVtxTrigEffDen;         //!
  private:
     AlidNdEtaTask(const AlidNdEtaTask&);
     AlidNdEtaTask& operator=(const AlidNdEtaTask&);
 
-  ClassDef(AlidNdEtaTask, 1);
+  ClassDef(AlidNdEtaTask, 2);
 };
 
 #endif
index ae390f2..5e16d08 100644 (file)
@@ -32,7 +32,9 @@ dNdEtaAnalysis::dNdEtaAnalysis() :
   fMult(0),
   fPtDist(0),
   fAnalysisMode(AliPWG0Helper::kInvalid),
-  fTag()
+  fTag(),
+  fvtxMin(-9.99),
+  fvtxMax(9.99)
 {
   // default constructor
 
@@ -50,7 +52,9 @@ dNdEtaAnalysis::dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0H
   fMult(0),
   fPtDist(0),
   fAnalysisMode(analysisMode),
-  fTag()
+  fTag(),
+  fvtxMin(-9.99),
+  fvtxMax(9.99)
 {
   // constructor
 
@@ -123,7 +127,9 @@ dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
   fMult(0),
   fPtDist(0),
   fAnalysisMode(AliPWG0Helper::kInvalid),
-  fTag()
+  fTag(),
+  fvtxMin(-9.99),
+  fvtxMax(9.99)
 {
   //
   // dNdEtaAnalysis copy constructor
@@ -433,8 +439,10 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     fdNdEtaPtCutOffCorrected[vertexPos]->Reset();
   }
 
-  const Float_t vertexRangeBegin[kVertexBinning] = { -9.99,  -9.99,  0.01 };
-  const Float_t vertexRangeEnd[kVertexBinning]   = {  9.99,  -0.01,  9.99 };
+  //const Float_t vertexRangeBegin[kVertexBinning] = { -9.99,  -9.99,  0.01 };
+  //const Float_t vertexRangeEnd[kVertexBinning]   = {  9.99,  -0.01,  9.99 };
+  const Float_t vertexRangeBegin[kVertexBinning] = { fvtxMin,  fvtxMin,  0.01 };
+  const Float_t vertexRangeEnd[kVertexBinning]   = { fvtxMax,  -0.01,  fvtxMax };
 
   for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
   {
@@ -445,31 +453,32 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       Int_t vertexBinEnd   = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
 
       const Int_t *binBegin = 0;
-      const Int_t maxBins = 40;
+      const Int_t maxBins = 20;
+      const Int_t maxBins1 = 40;
       
       if (dataHist->GetNbinsY() != maxBins)
         AliFatal(Form("Binning of acceptance is different from data histogram: data=%d, acceptance=%d", dataHist->GetNbinsY(), maxBins));
 
       // adjust acceptance range
       // produce with drawPlots.C: DetermineAcceptance(...)
-      const Int_t binBeginSPD[maxBins] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
-      const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1};
-      const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, -1, 14, 10, 8, 7, 6, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1};
+
+      //const Int_t binBeginSPD[maxBins] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
+      //const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1};
+      //const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, -1, 14, 10, 8, 7, 6, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1};
 
       if (fAnalysisMode & AliPWG0Helper::kSPD)
       {
-        //const Int_t binBeginSPD[maxBins] = {15, 14, 13, 12, 11, 10, 9, 9, 8, 7, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
-        
+       const Int_t binBeginSPD[maxBins1] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
         binBegin = binBeginSPD;
       }
       else if (fAnalysisMode & AliPWG0Helper::kTPC)
       {
-
+        const Int_t binBeginTPC[maxBins1] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1};
         binBegin = binBeginTPC;
       }
       else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
       {
-
+        const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1};
         binBegin = binBeginTPCITS;
       }
 
@@ -479,7 +488,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       if (binBegin)
       {
         vtxBegin = binBegin[iEta - 1];
-        vtxEnd = 18 + 1 - binBegin[maxBins - iEta];
+        vtxEnd = (Int_t) vtxVsEta->GetNbinsX() + 1 - binBegin[maxBins - iEta];
       }
       else
         Printf("WARNING: No acceptance applied!");
index 0d317ec..07760f2 100644 (file)
@@ -62,6 +62,9 @@ public:
   void SetAnalysisMode(AliPWG0Helper::AnalysisMode mode) { fAnalysisMode = mode; }
   AliPWG0Helper::AnalysisMode GetAnalysisMode() { return fAnalysisMode; }
 
+  void SetVtxMin(Float_t vtxMin) {fvtxMin = vtxMin;}
+  void SetVtxMax(Float_t vtxMax) {fvtxMax = vtxMax;}
+
 protected:
   Float_t GetVtxMin(Float_t eta);
 
@@ -76,7 +79,10 @@ protected:
   AliPWG0Helper::AnalysisMode fAnalysisMode;       // detector (combination) used for the analysis
   TString fTag;                                   // tag saved that describes the applied correction
 
-  ClassDef(dNdEtaAnalysis, 2)
+  Float_t fvtxMin;  // mininum vtx on which to integrate
+  Float_t fvtxMax;  // maximum vtx on which to intgrate
+
+  ClassDef(dNdEtaAnalysis, 3)
 };
 
 #endif