removing checking of vertex resolution in getvertex
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 10 Sep 2008 12:57:42 +0000 (12:57 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 10 Sep 2008 12:57:42 +0000 (12:57 +0000)
--> new function AliPWG0Helper::TestVertex
adding optional multiplicity axis for track level corrections for SPD

18 files changed:
PWG0/AliCorrection.cxx
PWG0/AliCorrectionMatrix.cxx
PWG0/AliCorrectionMatrix2D.cxx
PWG0/AliPWG0Helper.cxx
PWG0/AliPWG0Helper.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
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/drawSystematics.C
PWG0/dNdEta/drawSystematicsNew.C
PWG0/dNdEta/run.C
PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx
PWG0/multiplicity/correct.C
PWG0/multiplicity/plots.C

index 44198fa1fe5dcd4e668b827357d701e44d7c043b..89843e75ba54d0fd3611d4a8b8878fa90b5aee7a 100644 (file)
@@ -59,6 +59,10 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Hel
     return;
   }
 
+  Int_t nBinsN2 = 1;
+  Float_t binLimitsN2[]   = {-0.5, 1000};
+  //Int_t nBinsN2 = 21;
+  //Float_t binLimitsN2[]   = {-0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 1000.5};
   Int_t nBinsN = 22;
   Float_t binLimitsN[]   = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
   //Int_t nBinsN = 52;
@@ -77,15 +81,25 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title, AliPWG0Hel
 //  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};
 //  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 };
 
-  TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",18, binLimitsVtx, 30, binLimitsEta , nBinsPt, binLimitsPt);
-
   fEventCorr = new AliCorrectionMatrix2D("EventCorrection", Form("%s EventCorrection", fTitle.Data()), 18, binLimitsVtx, nBinsN, binLimitsN);
-  fTrackCorr = new AliCorrectionMatrix3D("TrackCorrection", Form("%s TrackCorrection", fTitle.Data()), dummyBinning);
 
-  delete dummyBinning;
+  if (analysisMode == AliPWG0Helper::kSPD)
+  {
+    TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",18, binLimitsVtx, 30, binLimitsEta, nBinsN2, binLimitsN2);
+    fTrackCorr = new AliCorrectionMatrix3D("TrackCorrection", Form("%s TrackCorrection", fTitle.Data()), dummyBinning);
+    fTrackCorr->SetAxisTitles("vtx z [cm]", "#eta", "Ntracks");
+    delete dummyBinning;
+  }
+  else
+  {
+    TH3F* dummyBinning = new TH3F("dummyBinning","dummyBinning",18, binLimitsVtx, 30, 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;
+  }
+
 
   fEventCorr->SetAxisTitles("vtx z [cm]", "Ntracks");
-  fTrackCorr->SetAxisTitles("vtx z [cm]", "#eta", "p_{T} [GeV/c]");
 }
 
 //____________________________________________________________________
@@ -404,7 +418,7 @@ void AliCorrection::PrintStats(Float_t zRange, Float_t etaRange, Float_t ptCut)
     tracksM /= eventsM;
     tracksG /= eventsG;
 
-    Printf("%.2f tracks/event measured, %.2f tracks/event after correction --> effective average correction factor is %.3f (pt cut %.2f GeV/c)", tracksM, tracksG, tracksG / tracksM, ptCut);
+    Printf("%.2f tracks/event measured, %.2f tracks/event after correction --> effective average correction factor is %.3f (pt cut %.2f GeV/c)", tracksM, tracksG, (tracksM > 0) ? (tracksG / tracksM) : -1, ptCut);
   }
 }
 
@@ -423,6 +437,8 @@ void AliCorrection::PrintInfo(Float_t ptCut)
 
   Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f", measured->Integral(), generated->Integral(), measuredEvents->Integral(), generatedEvents->Integral());
 
+  Printf("Example centered bin: tracks measured: %.1f tracks generated: %.1f, events measured: %.1f, events generated %.1f", measured->GetBinContent(measured->GetNbinsX() / 2, measured->GetNbinsY() / 2, measured->GetNbinsZ() / 2), generated->GetBinContent(measured->GetNbinsX() / 2, measured->GetNbinsY() / 2, measured->GetNbinsZ() / 2), measuredEvents->GetBinContent(measuredEvents->GetNbinsX() / 2, measuredEvents->GetNbinsY() / 2), generatedEvents->GetBinContent(measuredEvents->GetNbinsX() / 2, measuredEvents->GetNbinsY() / 2));
+
   PrintStats(10, 1.0, ptCut);
   PrintStats(10, 1.5, ptCut);
 }
index 82fac7236e0832df5fa3e49fb96093e55722df4c..6f283d7c0ee871c4c931247ff6cc8323e8a5d4e7 100644 (file)
@@ -353,9 +353,9 @@ void AliCorrectionMatrix::SetCorrectionToUnity()
   if (!fhCorr)
     return;
 
-  for (Int_t x=1; x<=fhCorr->GetNbinsX(); ++x)
-    for (Int_t y=1; y<=fhCorr->GetNbinsY(); ++y)
-      for (Int_t z=1; z<=fhCorr->GetNbinsZ(); ++z)
+  for (Int_t x=0; x<=fhCorr->GetNbinsX()+1; ++x)
+    for (Int_t y=0; y<=fhCorr->GetNbinsY()+1; ++y)
+      for (Int_t z=0; z<=fhCorr->GetNbinsZ()+1; ++z)
       {
         fhCorr->SetBinContent(x, y, z, 1);
         fhCorr->SetBinError(x, y, z, 0);
index 891de0b512d9fec6dde391f5046b8865fe8af8ed..6878e719ca8c71fc93ae41165fe1ccded342b5c7 100644 (file)
@@ -162,10 +162,15 @@ TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(Char_t* opt, Float_t min, F
 
   gene1D->SetName(Form("corr_1D_%s",fName.Data()));
   gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
+  TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
+  divided->Reset();
+  divided->Divide(gene1D, meas1D, 1, 1, "B");
 
-  gene1D->Divide(gene1D, meas1D, 1, 1, "B");
+  Printf("%p %p", gene1D, meas1D);
   
-  return (TH1F*)gene1D;   
+  return (TH1F*)divided;   
 }
 
 //____________________________________________________________________
index 9593f8e67dbb8034e3ef502bd0885fb5ee6b8cb6..160a9466ebd84f69a4dbd165e09c3e8ca74ed0cf 100644 (file)
@@ -78,6 +78,31 @@ Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
   return kFALSE;
 }
 
+//____________________________________________________________________
+const Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug)
+{
+    // Checks if a vertex meets the needed quality criteria
+
+  Float_t requiredZResolution = -1;
+  if (analysisMode == kSPD || analysisMode == kTPCITS)
+  {
+    requiredZResolution = 0.1;
+  }
+  else if (analysisMode == kTPC)
+    requiredZResolution = 10.;
+
+  // check resolution
+  Double_t zRes = vertex->GetZRes();
+
+  if (zRes > requiredZResolution) {
+    if (debug)
+      Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
+    return kFALSE;
+  }
+
+  return kTRUE;
+}
+
 //____________________________________________________________________
 const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug,Bool_t bRedoTPC)
 {
@@ -87,31 +112,28 @@ const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode ana
   // also the quality criteria that are applied)
 
   const AliESDVertex* vertex = 0;
-  Float_t requiredZResolution = -1;
   if (analysisMode == kSPD || analysisMode == kTPCITS)
   {
     vertex = aEsd->GetPrimaryVertexSPD();
-    requiredZResolution = 0.1;
     if (debug)
       Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
   }
-  else if (analysisMode == kTPC) 
+  else if (analysisMode == kTPC)
   {
     if(bRedoTPC){
-      Double_t kBz = aEsd->GetMagneticField(); 
+      Double_t kBz = aEsd->GetMagneticField();
       AliVertexerTracks vertexer(kBz);
-      vertexer.SetTPCMode(); 
+      vertexer.SetTPCMode();
       AliESDVertex *vTPC = vertexer.FindPrimaryVertex(aEsd);
       aEsd->SetPrimaryVertexTPC(vTPC);
       for (Int_t i=0; i<aEsd->GetNumberOfTracks(); i++) {
        AliESDtrack *t = aEsd->GetTrack(i);
        t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
-      } 
+      }
       delete vTPC;
     }
 
     vertex = aEsd->GetPrimaryVertexTPC();
-    requiredZResolution = 10.;
     if (debug)
       Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks");
   }
@@ -129,17 +151,15 @@ const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode ana
     if (debug){
       Printf("AliPWG0Helper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors());
       Printf("AliPWG0Helper::GetVertex: NIndices(): %d",vertex->GetNIndices());
-
+      vertex->Print();
     }
     return 0;
   }
 
   // check resolution
   Double_t zRes = vertex->GetZRes();
-
-  if (zRes == 0 || zRes > requiredZResolution) {
-    if (debug)
-      Printf("AliPWG0Helper::GetVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
+  if (zRes == 0) {
+    Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0.");
     return 0;
   }
 
index 55132c370fe0c36eef1db3bac3c89cb7ac67675a..022b2339e5239f74b0f22bd0722b294bb8574eb1 100644 (file)
@@ -30,6 +30,7 @@ class AliPWG0Helper : public TObject
     static Bool_t IsEventTriggered(const AliESD* aEsd, Trigger trigger = kMB2);
     static Bool_t IsEventTriggered(ULong64_t triggerMask, Trigger trigger = kMB2);
     static const AliESDVertex* GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMethod, Bool_t debug = kFALSE,Bool_t bRedoTPC = false);
+    static const Bool_t TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug = kFALSE);
 
     static Bool_t IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug = kFALSE);
 
index 457ef11b3d4e782f3fe828819d139b63c3a99e0d..a167463fe4e7b95cff032d48f648c3561fc3f008 100644 (file)
@@ -51,10 +51,12 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
   fPIDParticles(0),
   fPIDTracks(0),
   fVertexCorrelation(0),
+  fVertexCorrelationShift(0),
   fVertexProfile(0),
   fVertexShift(0),
   fVertexShiftNorm(0),
   fEtaCorrelation(0),
+  fEtaCorrelationShift(0),
   fEtaProfile(0),
   fEtaResolution(0),
   fpTResolution(0),
@@ -94,10 +96,12 @@ AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
   fPIDParticles(0),
   fPIDTracks(0),
   fVertexCorrelation(0),
+  fVertexCorrelationShift(0),
   fVertexProfile(0),
   fVertexShift(0),
   fVertexShiftNorm(0),
   fEtaCorrelation(0),
+  fEtaCorrelationShift(0),
   fEtaProfile(0),
   fEtaResolution(0),
   fpTResolution(0),
@@ -230,6 +234,8 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
 
   fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
   fOutput->Add(fVertexCorrelation);
+  fVertexCorrelationShift = new TH2F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx", 120, -30, 30, 100, -1, 1);
+  fOutput->Add(fVertexCorrelationShift);
   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
   fOutput->Add(fVertexProfile);
   fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
@@ -239,6 +245,8 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
 
   fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
   fOutput->Add(fEtaCorrelation);
+  fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta - ESD #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);
   fOutput->Add(fEtaProfile);
   fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2);
@@ -349,21 +357,28 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   {
     Double_t vtx[3];
     vtxESD->GetXYZ(vtx);
-    eventVertex = kTRUE;
 
     Double_t diff = vtxMC[2] - vtx[2];
     fVertexShift->Fill(diff);
     if (vtxESD->GetZRes() > 0)
-      fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
+        fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
 
-    if (eventTriggered)
+    if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
     {
-      fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
-      fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
+        vtxESD = 0;
+    }
+    else
+    {
+      eventVertex = kTRUE;
+
+      if (eventTriggered)
+      {
+        fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
+        fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
+        fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
+      }
     }
   }
-  else
-    Printf("No vertex found");
 
   // fill process type
   Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
@@ -564,28 +579,28 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       fPIDParticles->Fill(particle->GetPdgCode());
 
     Float_t eta = particle->Eta();
-    Float_t pt = particle->Pt();
+    Float_t thirdDim = (fAnalysisMode == AliPWG0Helper::kSPD) ? inputCount : particle->Pt();
 
-    fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
+    fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 
     if (fdNdEtaCorrectionProcessType[0])
     {
       // non diffractive
       if (processType==AliPWG0Helper::kND)
-        fdNdEtaCorrectionProcessType[0]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
+        fdNdEtaCorrectionProcessType[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 
       // single diffractive
       if (processType==AliPWG0Helper::kSD)
-        fdNdEtaCorrectionProcessType[1]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
+        fdNdEtaCorrectionProcessType[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
 
       // double diffractive
       if (processType==AliPWG0Helper::kDD)
-        fdNdEtaCorrectionProcessType[2]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
+        fdNdEtaCorrectionProcessType[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
     }
 
     if (eventTriggered)
       if (eventVertex)
-        fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
+        fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, thirdDim);
 
     // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
     if (TMath::Abs(eta) < 1.5) // && pt > 0.2)
@@ -658,21 +673,36 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       fEtaResolution->Fill(particle->Eta() - etaArr[i]);
       fpTResolution->Fill(particle->Pt() - ptArr[i]);
 
+      Float_t eta = -999;
+      Float_t thirdDim = -1;
+
       Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
       // in case of primary the MC values are filled, otherwise (background) the reconstructed values
       if (label == label2 && firstIsPrim)
       {
-        fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+        thirdDim = particle->Pt();
+        eta = particle->Eta();
       }
       else
       {
-        fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], etaArr[i], ptArr[i]);
+        thirdDim = ptArr[i];
+        eta = etaArr[i];
       }
 
-      fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
-      fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+      if (fAnalysisMode == AliPWG0Helper::kSPD)
+        thirdDim = inputCount;
+
+      fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim);
+
+      // eta comparison for tracklets with the same label (others are background)
+      if (label == label2)
+      {
+        fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
+        fEtaCorrelation->Fill(etaArr[i], particle->Eta());
+        fEtaCorrelationShift->Fill(etaArr[i], particle->Eta() - etaArr[i]);
+      }
 
-      fEtaCorrelation->Fill(etaArr[i], particle->Eta());
+      fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
 
       if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
       {
@@ -683,15 +713,15 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
       {
         // non diffractive
         if (processType == AliPWG0Helper::kND)
-          fdNdEtaCorrectionProcessType[0]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+          fdNdEtaCorrectionProcessType[0]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
 
         // single diffractive
         if (processType == AliPWG0Helper::kSD)
-          fdNdEtaCorrectionProcessType[1]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+          fdNdEtaCorrectionProcessType[1]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
 
         // double diffractive
         if (processType == AliPWG0Helper::kDD)
-          fdNdEtaCorrectionProcessType[2]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+          fdNdEtaCorrectionProcessType[2]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
       }
 
       // control histograms
@@ -845,6 +875,9 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
   if (fVertexCorrelation)
     fVertexCorrelation->Write();
+  fVertexCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelationShift"));
+  if (fVertexCorrelationShift)
+    fVertexCorrelationShift->Write();
   fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
   if (fVertexProfile)
     fVertexProfile->Write();
@@ -858,6 +891,9 @@ void AlidNdEtaCorrectionTask::Terminate(Option_t *)
   fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
   if (fEtaCorrelation)
     fEtaCorrelation->Write();
+  fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
+  if (fEtaCorrelationShift)
+    fEtaCorrelationShift->Write();
   fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
   if (fEtaProfile)
     fEtaProfile->Write();
index fe6a0e357fd3ea782e1164476d41c34dcbfdd4a3..ddc15a865804dbc7693ad69e024922a071233f35 100644 (file)
@@ -59,11 +59,13 @@ class AlidNdEtaCorrectionTask : public AliAnalysisTask {
     TH1F* fPIDTracks;                            //! pid of reconstructed tracks
 
     TH2F* fVertexCorrelation;                    //! ESD z-vtx vs MC z-vtx
+    TH2F* fVertexCorrelationShift;               //! (MC z-vtx - ESD z-vtx) vs MC z-vtx
     TProfile* fVertexProfile;                    //! Profile of MC z-vtx - ESD z-vtx vs. MC z-vtx
     TH1F* fVertexShift;                          //! (MC z-vtx - ESD z-vtx) in +- 10 cm
     TH1F* fVertexShiftNorm;                      //! (MC z-vtx - ESD z-vtx) / (sigma_ESD-z-vtx) histogrammed
 
     TH2F* fEtaCorrelation;                       //! ESD eta vs MC eta
+    TH2F* fEtaCorrelationShift;                  //! (MC eta - ESD eta) vs MC eta
     TProfile* fEtaProfile;                       //! Profile of MC eta - ESD eta vs. MC eta
     TH1F* fEtaResolution;                        //! MC eta - ESD eta in |eta| < 1
 
index c5f46bac77e2cd1e0ef967a859317e7a12e8a086..0cf621a3b2718360c2bd219ec10a50302bc2b28e 100644 (file)
@@ -59,8 +59,8 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fdNdEtaAnalysisTr(0),
   fdNdEtaAnalysisTrVtx(0),
   fdNdEtaAnalysisTracks(0),
-  fVertex(0),
   fPartPt(0),
+  fVertex(0),
   fPhi(0),
   fEtaPhi(0),
   fDeltaPhi(0)
@@ -163,7 +163,7 @@ void AlidNdEtaTask::CreateOutputObjects()
     fOutput->Add(fPartEta[i]);
   }
 
-  fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 160, -40, 40);
+  fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
   fOutput->Add(fEvents);
 
   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
@@ -178,6 +178,9 @@ void AlidNdEtaTask::CreateOutputObjects()
   if (fAnalysisMode == AliPWG0Helper::kSPD)
     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 18*50, -3.14, 3.14);
 
+  fVertex = new TH3F("vertex_check", "vertex_check", 100, -1, 1, 100, -1, 1, 100, -30, 30);
+  fOutput->Add(fVertex);
+
   if (fReadMC)
   {
     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
@@ -195,9 +198,6 @@ void AlidNdEtaTask::CreateOutputObjects()
     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
     fOutput->Add(fdNdEtaAnalysisTracks);
 
-    fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
-    fOutput->Add(fVertex);
-
     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
     fPartPt->Sumw2();
     fOutput->Add(fPartPt);
@@ -218,6 +218,9 @@ void AlidNdEtaTask::Exec(Option_t*)
   Bool_t eventTriggered = kFALSE;
   const AliESDVertex* vtxESD = 0;
 
+  // post the data already here
+  PostData(0, fOutput);
+
   // ESD analysis
   if (fESD)
   {
@@ -228,18 +231,20 @@ void AlidNdEtaTask::Exec(Option_t*)
     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 
     Double_t vtx[3];
-    if (vtxESD) {
-      vtxESD->GetXYZ(vtx);
-    }
-    else
-      Printf("No vertex");
 
     // fill z vertex resolution
     if (vtxESD)
+    {
       fVertexResolution->Fill(vtxESD->GetZRes());
+      fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
 
-    // post the data already here
-    PostData(0, fOutput);
+      if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
+      {
+          vtxESD->GetXYZ(vtx);
+      }
+      else
+        vtxESD = 0;
+    }
 
     // needed for syst. studies
     AliStack* stack = 0;
@@ -336,14 +341,30 @@ void AlidNdEtaTask::Exec(Option_t*)
 
         fDeltaPhi->Fill(deltaPhi);
 
-        // TODO implement fUseMCKine here
+        Int_t label = mult->GetLabel(i, 0);
+        Float_t eta = mult->GetEta(i);
+
+        if (fUseMCKine)
+        {
+          if (label > 0)
+          {
+            TParticle* particle = stack->Particle(label);
+            eta = particle->Eta();
+          }
+          else
+            Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
+        }
 
-        etaArr[inputCount] = mult->GetEta(i);
-        labelArr[inputCount] = mult->GetLabel(i, 0);
+        etaArr[inputCount] = eta;
+        labelArr[inputCount] = label;
         ptArr[inputCount] = 0; // no pt for tracklets
         ++inputCount;
       }
 
+      // fill multiplicity in pt bin
+      for (Int_t i=0; i<inputCount; ++i)
+        ptArr[i] = inputCount;
+
       Printf("Accepted %d tracklets", inputCount);
     }
     else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
@@ -465,7 +486,7 @@ void AlidNdEtaTask::Exec(Option_t*)
               continue;
             }
 
-            fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+            fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), (fAnalysisMode == AliPWG0Helper::kSPD) ? inputCount : particle->Pt());
           } // end of track loop
 
           // for event count per vertex
@@ -530,6 +551,7 @@ void AlidNdEtaTask::Exec(Option_t*)
 
     Int_t nAcceptedParticles = 0;
 
+    // count particles first, then fill
     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
     {
       if (!stack->IsPhysicalPrimary(iMc))
@@ -546,27 +568,43 @@ void AlidNdEtaTask::Exec(Option_t*)
 
       //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
       Float_t eta = particle->Eta();
-      Float_t pt = particle->Pt();
 
        // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
       if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
         nAcceptedParticles++;
+    }
+
+    for (Int_t iMc = 0; iMc < nPrim; ++iMc)
+    {
+      if (!stack->IsPhysicalPrimary(iMc))
+        continue;
 
-      fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
-      fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
+      //Printf("Getting particle %d", iMc);
+      TParticle* particle = stack->Particle(iMc);
+
+      if (!particle)
+        continue;
+
+      if (particle->GetPDG()->Charge() == 0)
+        continue;
+
+      Float_t eta = particle->Eta();
+      Float_t thirdDim = (fAnalysisMode == AliPWG0Helper::kSPD) ? nAcceptedParticles : particle->Pt();
+
+      fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
 
       if (processType != AliPWG0Helper::kSD )
-        fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, pt);
+        fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
 
       if (eventTriggered)
       {
-        fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
+        fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
         if (vtxESD)
-          fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
+          fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
       }
 
       if (TMath::Abs(eta) < 0.8)
-        fPartPt->Fill(pt);
+        fPartPt->Fill(particle->Pt());
     }
 
     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
@@ -612,85 +650,90 @@ void AlidNdEtaTask::Terminate(Option_t *)
   if (!fdNdEtaAnalysisESD)
   {
     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
-    return;
   }
-
-  if (fMult && fMultVtx)
+  else
   {
-    new TCanvas;
-    fMult->Draw();
-    fMultVtx->SetLineColor(2);
-    fMultVtx->Draw("SAME");
-  }
+    if (fMult && fMultVtx)
+    {
+      new TCanvas;
+      fMult->Draw();
+      fMultVtx->SetLineColor(2);
+      fMultVtx->Draw("SAME");
+    }
 
-  if (fPartEta[0])
-  {
-    Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
-    Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
+    if (fPartEta[0])
+    {
+      Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
+      Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
 
-    Printf("%d events with vertex used", events1 + events2);
+      Printf("%d events with vertex used", events1 + events2);
 
-    if (events1 > 0 && events2 > 0)
-    {
-    fPartEta[0]->Scale(1.0 / (events1 + events2));
-    fPartEta[1]->Scale(1.0 / events1);
-    fPartEta[2]->Scale(1.0 / events2);
+      if (events1 > 0 && events2 > 0)
+      {
+        fPartEta[0]->Scale(1.0 / (events1 + events2));
+        fPartEta[1]->Scale(1.0 / events1);
+        fPartEta[2]->Scale(1.0 / events2);
 
-    for (Int_t i=0; i<3; ++i)
-      fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
+        for (Int_t i=0; i<3; ++i)
+          fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
 
-    new TCanvas("control", "control", 500, 500);
-    for (Int_t i=0; i<3; ++i)
-    {
-      fPartEta[i]->SetLineColor(i+1);
-      fPartEta[i]->Draw((i==0) ? "" : "SAME");
+        new TCanvas("control", "control", 500, 500);
+        for (Int_t i=0; i<3; ++i)
+        {
+          fPartEta[i]->SetLineColor(i+1);
+          fPartEta[i]->Draw((i==0) ? "" : "SAME");
+        }
+      }
     }
+
+    if (fEvents)
+    {
+        new TCanvas("control3", "control3", 500, 500);
+        fEvents->Draw();
     }
-  }
 
-  if (fEvents)
-  {
-    new TCanvas("control3", "control3", 500, 500);
-    fEvents->Draw();
-  }
+    TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
 
-  TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
+    if (fdNdEtaAnalysisESD)
+      fdNdEtaAnalysisESD->SaveHistograms();
 
-  if (fdNdEtaAnalysisESD)
-    fdNdEtaAnalysisESD->SaveHistograms();
+    if (fEsdTrackCuts)
+      fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
 
-  if (fEsdTrackCuts)
-    fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
+    if (fMult)
+      fMult->Write();
 
-  if (fMult)
-    fMult->Write();
+    if (fMultVtx)
+      fMultVtx->Write();
 
-  if (fMultVtx)
-    fMultVtx->Write();
+    for (Int_t i=0; i<3; ++i)
+      if (fPartEta[i])
+        fPartEta[i]->Write();
 
-  for (Int_t i=0; i<3; ++i)
-    if (fPartEta[i])
-      fPartEta[i]->Write();
+    if (fEvents)
+      fEvents->Write();
 
-  if (fEvents)
-    fEvents->Write();
+    if (fVertexResolution)
+      fVertexResolution->Write();
 
-  if (fVertexResolution)
-    fVertexResolution->Write();
+    if (fDeltaPhi)
+      fDeltaPhi->Write();
 
-  if (fDeltaPhi)
-    fDeltaPhi->Write();
+    if (fPhi)
+      fPhi->Write();
 
-  if (fPhi)
-    fPhi->Write();
+    if (fEtaPhi)
+      fEtaPhi->Write();
 
-  if (fEtaPhi)
-    fEtaPhi->Write();
+    fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
+    if (fVertex)
+      fVertex->Write();
 
-  fout->Write();
-  fout->Close();
+    fout->Write();
+    fout->Close();
 
-  Printf("Writting result to analysis_esd_raw.root");
+    Printf("Writting result to analysis_esd_raw.root");
+  }
 
   if (fReadMC)
   {
@@ -702,7 +745,6 @@ void AlidNdEtaTask::Terminate(Option_t *)
       fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
       fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
       fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
-      fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
     }
 
     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
@@ -721,7 +763,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
     fPartPt->Scale(1.0/events);
     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
 
-    fout = new TFile("analysis_mc.root","RECREATE");
+    TFile* fout = new TFile("analysis_mc.root","RECREATE");
 
     fdNdEtaAnalysis->SaveHistograms();
     fdNdEtaAnalysisNSD->SaveHistograms();
@@ -732,9 +774,6 @@ void AlidNdEtaTask::Terminate(Option_t *)
     if (fPartPt)
       fPartPt->Write();
 
-    if (fVertex)
-      fVertex->Write();
-
     fout->Write();
     fout->Close();
 
index ee612ec329e88a153c0d91aa17c0449428650254..f61fa830a84f1a7c4d198349c7b1d8f05770fe80 100644 (file)
@@ -67,10 +67,10 @@ class AlidNdEtaTask : public AliAnalysisTask {
     dNdEtaAnalysis* fdNdEtaAnalysisTracks;  //! contains the dndeta from the triggered events with vertex counted from the mc particles associated to the tracks (comparing this to the raw values from the esd shows the effect of the detector resolution)
 
     // control histograms (MC)
-    TH3F* fVertex;                //! vertex of counted particles
     TH1F* fPartPt;                //! counted particles as function of pt
 
     // control histograms (ESD)
+    TH3F* fVertex;                //! 3d vertex distribution
     TH1F* fPhi;                   //! raw phi distribution
     TH2F* fEtaPhi;                //! raw eta - phi distribution
     TH1F* fDeltaPhi;              //! histogram of delta_phi values for tracklets (only for SPD analysis)
index f1537d385f7315c4b7eff3606e6849088be6463d..fe2ebe92e5a3b156984a70c28a1bcc085b4d3be6 100644 (file)
@@ -201,7 +201,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
   // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
   //
 
-  fTag.Form("Correcting dN/deta spectrum >>> %s <<<. Correction type: %d, pt cut: %.2f.", tag, (Int_t) correctionType, ptCut);
+  fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut);
   Printf("\n\n%s", fTag.Data());
 
   // set corrections to 1
@@ -268,6 +268,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
     TH2* eTrig =    correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
     TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
+    //TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1);
     TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1);
 
     //new TCanvas; correctedEvents->DrawCopy("TEXT");
@@ -281,6 +282,9 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
     Printf("%d triggered events with 0 mult. -- %d triggered events with vertex", triggeredEventsWith0Mult, allEventsWithVertex);
 
+    TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
+    kineBias->Reset();
+
     for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
     {
       Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
@@ -291,6 +295,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
       Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
         eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
+      kineBias->SetBinContent(i, fZ);
 
       events *= fZ;
 
@@ -303,6 +308,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     }
 
     //new TCanvas; correctedEvents->DrawCopy("TEXT");
+    //new TCanvas; kineBias->Draw();
   }
 
   fData->PrintInfo(ptCut);
@@ -315,6 +321,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
   TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
 
   // create pt hist
+  if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
   {
     // reset all ranges
     dataHist->GetXaxis()->SetRange(0, 0);
@@ -353,13 +360,13 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
   dataHist->GetYaxis()->SetRange(0, 0);
   dataHist->GetZaxis()->SetRange(0, 0);
 
-  // integrate over pt (with pt cut)
+  // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD)
   Int_t ptLowBin = 1;
-  if (ptCut > 0)
+  if (ptCut > 0 && fAnalysisMode != AliPWG0Helper::kSPD)
     ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
 
   dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
-  printf("pt range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
+  printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
   TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
 
   dataHist->GetZaxis()->SetRange(0, 0);
@@ -368,10 +375,12 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
   if (vtxVsEta == 0)
   {
-    printf("ERROR: pt integration failed\n");
+    printf("ERROR: pt/multiplicity integration failed\n");
     return;
   }
 
+  //new TCanvas(tag, tag, 800, 600);  vtxVsEta->DrawCopy("COLZ");
+
   // clear result histograms
   for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
   {
@@ -415,8 +424,16 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         // TODO create map
       }
 
-      Int_t vtxBegin = binBegin[iEta - 1];
-      Int_t vtxEnd   = 18 + 1 - binBegin[30 - iEta];
+      Int_t vtxBegin = 1;
+      Int_t vtxEnd   = 30;
+
+      if (binBegin)
+      {
+        vtxBegin = binBegin[iEta - 1];
+        vtxEnd = 18 + 1 - binBegin[30 - iEta];
+      }
+      else
+        Printf("WARNING: No acceptance applied!");
       
       // eta range not accessible
       if (vtxBegin == -1)
@@ -427,10 +444,10 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta));
 
       Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
-      printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
+      //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
       vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
       vertexBinEnd =   TMath::Min(vertexBinEnd, vtxEnd);
-      Printf(" after:  %d %d", vertexBinBegin, vertexBinEnd);
+      //Printf(" after:  %d %d", vertexBinBegin, vertexBinEnd);
 
       // no data for this bin
       if (vertexBinBegin > vertexBinEnd)
@@ -473,13 +490,18 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       }
 
       Float_t ptCutOffCorrection = 1;
-      if (correction && ptCut > 0)
-        ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd);
 
-      if (ptCutOffCorrection <= 0)
+      // find pt cut off correction factor
+      if (fAnalysisMode != AliPWG0Helper::kSPD)
       {
-        printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
-        continue;
+        if (correction && ptCut > 0)
+            ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd);
+
+        if (ptCutOffCorrection <= 0)
+        {
+            printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
+            continue;
+        }
       }
 
       //printf("Eta: %d Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f\n", iEta, vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection);
@@ -502,7 +524,7 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
 
-        Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents);
+        //Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents);
       }
     }
   }
index bebc840d1d1c2a35a190e2dc795b2c29e90a6f01..38475da8d11d0afbdd15b8d5da2a4936fd19658f 100644 (file)
@@ -59,6 +59,8 @@ public:
   TH1F* GetdNdEtaHistogram(Int_t i = 0) { return fdNdEta[i]; }
   TH1F* GetdNdEtaPtCutOffCorrectedHistogram(Int_t i = 0) { return fdNdEtaPtCutOffCorrected[i]; }
 
+  void SetAnalysisMode(AliPWG0Helper::AnalysisMode mode) { fAnalysisMode = mode; }
+
 protected:
   Float_t GetVtxMin(Float_t eta);
 
index 1590528538f587e155855bdd7d0a184ffd6a1680..436a9238b296c917f86685e8b4e2cb4709f446a1 100644 (file)
@@ -118,6 +118,48 @@ void DrawOverview(const char* fileName = "correction_map.root", const char* dirN
   dNdEtaCorrection->DrawOverview();
 }
 
+void PrintInfo(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
+{
+  loadlibs();
+  if (!TFile::Open(fileName))
+    return;
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
+  if (!dNdEtaCorrection->LoadHistograms())
+    return;
+
+  dNdEtaCorrection->Finish();
+
+  for (Int_t i=AlidNdEtaCorrection::kTrack2Particle; i<=AlidNdEtaCorrection::kND; i++)
+  {
+    Printf("Correction %d", i);
+    dNdEtaCorrection->GetCorrection(i)->PrintInfo(0.3);
+  }
+}
+
+void PrintAllInfos()
+{
+  PrintInfo();
+
+  Printf("RAW ESD");
+  TFile::Open("analysis_esd_raw.root");
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+  fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
+
+  const Int_t num = 3;
+  const char* results[] = { "dndeta", "dndetaTr", "dndetaTrVtx" };
+
+  TFile::Open("analysis_esd.root");
+  for (Int_t i=0; i<num; i++)
+  {
+    Printf("CORRECTED %s", results[i]);
+    dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+    fdNdEtaAnalysis->LoadHistograms(results[i]);
+    fdNdEtaAnalysis->GetData()->PrintInfo(0.3);
+  }
+}  
+
 void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
 {
   gSystem->Load("libPWG0base");
@@ -455,6 +497,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   DrawdNdEtaRatio(histESD, histMC, "full_inelastic", etaPlotLimit);
   DrawdNdEtaRatio(histESDMB, histMCTr, "triggered", etaPlotLimit);
   DrawdNdEtaRatio(histESDMBVtx, histMCTrVtx, "triggered_vertex", etaPlotLimit);
+  DrawdNdEtaRatio(histESDnsd, histMCnsd, "NSD", etaPlotLimit);
 
   new TCanvas;
   dummy2->DrawCopy();
@@ -480,6 +523,7 @@ void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
   Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
 
   PrintIntegratedDeviation(histMC, histESD, "all events");
+  PrintIntegratedDeviation(histMCnsd, histESDnsd, "all events (NSD)");
   PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
   PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
   PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
@@ -1184,6 +1228,44 @@ void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
   }
 }
 
+void CompareTrack2Particle1D(const char* file1, const char* file2, Float_t upperPtLimit = 9.9)
+{
+  loadlibs();
+
+  const char* folderName = "dndeta_correction";
+
+  c = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
+  c->Divide(3, 1);
+
+  for (Int_t fileId = 0; fileId < 2; fileId++)
+  {
+    const char* file = ((fileId == 0) ? file1 : file2);
+    Correction1DCreatePlots(file, folderName, upperPtLimit, 1);
+
+    TH1* corr[3];
+    corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
+    corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
+    corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
+    /*corr[0] = dynamic_cast<TH1*> (gROOT->FindObject("generated_x"))->Clone(Form("hist_x_%d", fileId));
+    corr[1] = dynamic_cast<TH1*> (gROOT->FindObject("generated_y"))->Clone(Form("hist_y_%d", fileId));
+    corr[2] = dynamic_cast<TH1*> (gROOT->FindObject("generated_z"))->Clone(Form("hist_z_%d", fileId));*/
+
+    for (Int_t i=0; i<3; i++)
+    {
+      c->cd(i+1);
+      InitPad();
+      corr[i]->GetYaxis()->SetRangeUser(0.8, 2);
+      corr[i]->SetLineColor(fileId+1);
+      corr[i]->DrawCopy((fileId == 0) ? "" : "SAME");
+    }
+  }
+
+  return;
+
+  c->SaveAs(Form("%s.gif", canvas->GetName()));
+  c->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
 {
   TFile::Open(fileName);
@@ -1246,8 +1328,8 @@ TCanvas* Track2Particle2D(const char* fileName = "correction_map.root", const ch
   InitPadCOLZ();
   corrZY->Draw("COLZ");
 
-  canvas->SaveAs(Form("spd_corr_track2particle_%d.gif", gMax));
-  canvas->SaveAs(Form("spd_corr_track2particle_%d.eps", gMax));
+  canvas->SaveAs(Form("corr_track2particle_%d.gif", gMax));
+  canvas->SaveAs(Form("corr_track2particle_%d.eps", gMax));
 
   return canvas;
 }
@@ -1976,3 +2058,45 @@ void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root",
 
   input->Multiply(proj);
 }
+
+void MakeGaussianProfile(const char* histName = "fVertexCorrelation", Bool_t subtractMean = kFALSE)
+{
+    TFile::Open("correction_map.root");
+
+    TH2* hist2d = (TH2*) gFile->Get(histName);
+    hist2d->Sumw2();
+
+    TH1* result = hist2d->ProjectionX("result");
+    result->GetYaxis()->SetTitle(hist2d->GetYaxis()->GetTitle());
+    result->Reset();
+
+    for (Int_t x=1; x<hist2d->GetNbinsX(); ++x)
+    {
+        hist = hist2d->ProjectionY(Form("temp_%d", x), x, x);
+        if (hist->GetEntries() == 0)
+            continue;
+        if (hist->Fit("gaus") == 0)
+        {
+            func = hist->GetFunction("gaus");
+            mean = func->GetParameter(1);
+            error = func->GetParError(1);
+
+            if (subtractMean)
+                mean = hist2d->GetXaxis()->GetBinCenter(x) - mean;
+
+            result->SetBinContent(x, mean);
+            result->SetBinError(x, error);
+
+            if (x % 10 == 0)
+            {
+                new TCanvas;
+                ((TH1*) hist->Clone())->DrawCopy();
+            }
+        }
+        //break;
+    }
+
+    new TCanvas;
+    result->GetYaxis()->SetRangeUser(-0.2, 0.2);
+    result->Draw();
+}
index 0d55f9afe5d54ebe23ffe930bd9e76be3f239987..843128161a20af2122f8045d0a95ed59916ef8dd 100644 (file)
@@ -728,7 +728,7 @@ void DrawdNdEtaDifferences()
   canvas2->SaveAs("particlecomposition_result.eps");
 }
 
-mergeCorrectionsWithDifferentCrosssections(Char_t* correctionFileName="correction_map.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
+void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
   //
   // Function used to merge standard corrections with vertex
   // reconstruction corrections obtained by a certain mix of ND, DD
@@ -738,26 +738,32 @@ mergeCorrectionsWithDifferentCrosssections(Char_t* correctionFileName="correctio
   // (standard to changed x-section) of the different dN/deta
   // distributions are saved to a file.
   //
+  // correctionTarget is of type AlidNdEtaCorrection::CorrectionType
+  //    kINEL = 3
+  //    kNSD = 4
 
   loadlibs();
 
   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
 
-  /*
+
+
   const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25"};
-  // add scalesND!!!
+  Float_t scalesND[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0 };
   Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
   Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 0.5, 1.5, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
   Int_t nChanges = 17;
-  */
 
+  /*
   const Char_t* changes[]  = { "pythia", "qgsm", "phojet"};
   Float_t scalesND[] = {1.0, 1.10, 1.11};
   Float_t scalesSD[] = {1.0, 0.69, 0.86};
   Float_t scalesDD[] = {1.0, 0.98, 0.61};
   Int_t nChanges = 3;
-
+  */
+  
   // cross section from Pythia
+  // 14 TeV!
   Float_t sigmaND = 55.2;
   Float_t sigmaDD = 9.78;
   Float_t sigmaSD = 14.30;
@@ -852,8 +858,7 @@ mergeCorrectionsWithDifferentCrosssections(Char_t* correctionFileName="correctio
       dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
       fdNdEtaAnalysis->LoadHistograms();
 
-      //fdNdEtaAnalysis->Finish(current, 0.3, AlidNdEtaCorrection::kINEL, Form("%d %d", j, i));
-      fdNdEtaAnalysis->Finish(current, 0.3, AlidNdEtaCorrection::kNSD);
+      fdNdEtaAnalysis->Finish(current, 0.3, correctionTarget, Form("%d %d", j, i));
 
       name = "ratio";
       if (j==0) name.Append("_vetexReco_");
index 0f9aede4ff28741ba7a404e6753a669e61e76b58..3fc4c245e61e6be62a760c93a4c228c82e5d6f09 100644 (file)
@@ -65,18 +65,19 @@ void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
 }
 
 
-void DrawEffectOfChangeInCrossSection(const char* fileName) {
+void DrawEffectOfChangeInCrossSection(const char* fileName = "systematics_vtxtrigger_compositions.root") {
   
   //get the data
   TFile* fin = TFile::Open(fileName);
 
-  //const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore" };
-  const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
-  const Int_t nChanges = 3;
+  const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore" };
   //const Char_t* changes[]  = {"pythia","ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
+  const Int_t nChanges = 9;
+  //const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
+  //const Int_t nChanges = 3;
   Int_t colors[] = {1,1,kRed,kBlue,kGreen,kPink,1,kRed,kBlue};
 
-  TH1F* hRatios[3];
+  TH1F* hRatios[9];
   for(Int_t i=0; i<nChanges; i++) {
     hRatios[i] = (TH1F*)fin->Get(Form("ratio_vertexReco_triggerBias_%s",changes[i]));
     hRatios[i]->SetLineWidth(2);
@@ -97,7 +98,7 @@ void DrawEffectOfChangeInCrossSection(const char* fileName) {
   null->GetYaxis()->SetTitle("Ratio pythia/modified cross-sections");
   null->Draw();
 
-  TLatex* text[3];
+  TLatex* text[9];
 
   for(Int_t i=1; i<nChanges; i++) {
     hRatios[i]->Draw("same");
index f96e56162025e9c276d6e0f15907080d79eee473..0a8a8621c39a2beae9bfc37d88366381cbbbba24 100644 (file)
@@ -1,9 +1,9 @@
 void Load(const char* taskName, Bool_t debug)
 {
   TString compileTaskName;
-  compileTaskName.Form("%s.cxx+", taskName);
+  compileTaskName.Form("%s.cxx++", taskName);
   if (debug)
-    compileTaskName += "+g";
+    compileTaskName += "g";
 
   if (gProof) {
     gProof->Load(compileTaskName);
@@ -120,16 +120,17 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     if (mc)
       task->SetReadMC();
 
+    // syst. error flags
     //task->SetUseMCVertex();
     //task->SetUseMCKine();
     //task->SetOnlyPrimaries();
-  
+
     task->SetTrigger(trigger);
     task->SetAnalysisMode(analysisMode);
     task->SetTrackCuts(esdTrackCuts);
 
     mgr->AddTask(task);
+
     // Attach input
     mgr->ConnectInput(task, 0, cInput);
 
@@ -141,6 +142,8 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
   {
     Load("AlidNdEtaCorrectionTask", aDebug);
     task2 = new AlidNdEtaCorrectionTask(option);
+
+    // syst. error flags
     //task2->SetOnlyPrimaries();
 
     task2->SetTrigger(trigger);
@@ -148,10 +151,10 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
     task2->SetTrackCuts(esdTrackCuts);
 
     mgr->AddTask(task2);
-    
+
     // Attach input
     mgr->ConnectInput(task2, 0, cInput);
-    
+
     // Attach output
     cOutput = mgr->CreateContainer("cOutput2", TList::Class(), AliAnalysisManager::kOutputContainer);
     mgr->ConnectOutput(task2, 0, cOutput);
@@ -183,10 +186,18 @@ void run(Int_t runWhat, const Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool
 
     mgr->StartAnalysis("proof", data, nRuns, offset);
   }
+  else if (aProof == 3)
+  {
+    ROOT->ProcessLine(".L CreateChainFromDataSet.C");
+    ds = gProof->GetDataSet(data);
+    chain = CreateChainFromDataSet(ds);
+    mgr->StartAnalysis("local", chain, nRuns, offset);
+  }
   else
   {
     // Create chain of input files
     gROOT->LoadMacro("../CreateESDChain.C");
+
     chain = CreateESDChain(data, nRuns, offset);
     //chain = CreateChain("TE", data, nRuns, offset);
 
index c778f410ffceb194cfc4c565bc02426ae0629053..c9f278511fec785fd0934d6f2988c9f959a611b5 100644 (file)
@@ -888,13 +888,13 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
   xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
   xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
 
-  // 10^28 lum --> 1.2 kHz
-  // 10^31 lum --> 1200 kHz
-  //Float_t rate = 1200e3;
-  Float_t rate = 250e3;
+  // 10^28 lum --> 1.4 kHz
+  // 10^31 lum --> 1400 kHz
+  //Float_t rate = 1400e3;
+  Float_t rate = 1.4e3;
 
   // time in s
-  Float_t lengthRun = 1e6;
+  Float_t lengthRun = 1e5;
 
   Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
   Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
@@ -918,11 +918,14 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
     //Float_t ratePerTrigger[] = { 60, 13.3, 13.3, 13.3 };
 
     Int_t nCuts = 3;
-    Int_t cuts[] = { 0, 126, 162 };
-    Float_t ratePerTrigger[] = { 60, 20.0, 20.0 };
+    Int_t cuts[] = { 0, 114, 148 };
+
+    //Int_t nCuts = 3;
+    //Int_t cuts[] = { 0, 126, 162 };
+    //Float_t ratePerTrigger[] = { 60, 20.0, 20.0 };
 
     // desired trigger rate in Hz
-    //Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
+    Float_t ratePerTrigger[] = { 100, 1, 1, 1, 1, 1 };
 
     xSection->SetStats(kFALSE);
     xSection->SetTitle(""); //(i == 0) ? "SPD Layer 1" : "SPD Layer 2");
@@ -979,6 +982,9 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
       nTrigger = TMath::Nint(((Float_t) nTrigger) / 1000) * 1000;
 
       Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
+      if (nTrigger == 0)
+        continue;
+
       proj2->FillRandom(proj, nTrigger);
 
       Int_t removed = 0;
index 268a26d0710a9dc1ad96666a6f7a94c902f5420d..d4a75185a39366e1914f98d5787ade02d7db60c1 100644 (file)
@@ -133,8 +133,8 @@ void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", c
   histB->SetLineColor(2);
   histRAW->SetLineColor(3);
 
-  histC->DrawCopy();
-  histB->DrawCopy("SAME");
+  histC->DrawCopy("HISTE");
+  histB->DrawCopy("HISTE SAME");
   histRAW->DrawCopy("SAME");
 
   legend = new TLegend(0.2, 0.2, 0.4, 0.4);
@@ -144,6 +144,7 @@ void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", c
   legend->AddEntry(histB, label2);
   legend->AddEntry(histRAW, "raw ESD");
 
+  histGraph = 0;
   if (simpleCorrect > 0)
   {
     graph = new TGraph;
@@ -189,7 +190,7 @@ void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", c
     histMC->Sumw2();
     histMC->Scale(1.0 / histMC->Integral());
 
-    histMC->Draw("SAME");
+    histMC->Draw("HISTE SAME");
     histMC->SetLineColor(4);
     legend->AddEntry(histMC, "MC");
   }
@@ -216,12 +217,18 @@ void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", c
     {
       histC->SetBinContent(bin, histC->GetBinContent(bin) / histMC->GetBinContent(bin));
       histB->SetBinContent(bin, histB->GetBinContent(bin) / histMC->GetBinContent(bin));
-      histGraph->SetBinContent(bin, histGraph->GetBinContent(bin) / histMC->GetBinContent(bin));
+
+      /*
+      if (histGraph)
+      {
+        histGraph->SetBinContent(bin, histGraph->GetBinContent(bin) / histMC->GetBinContent(bin));
+        histGraph->SetBinError(bin, 0);
+      }
+      */
 
       // TODO errors?
       histC->SetBinError(bin, 0);
       histB->SetBinError(bin, 0);
-      histGraph->SetBinError(bin, 0);
     }
   }
 
@@ -230,7 +237,10 @@ void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", c
 
   histC->Draw("HIST");
   histB->Draw("HIST SAME");
-  histGraph->Draw("HIST SAME");
+  /*
+  if (histGraph)
+    histGraph->Draw("HIST SAME");
+  */
 
   /*
   if (simpleCorrect > 0)
index 1894423ef98a58af6476cd03ffbf097eb2d52047..24f6e249a9d3420a5bf522c4fa70e4ca26ca616b 100644 (file)
@@ -30,9 +30,9 @@
 
 #endif
 
-const char* correctionFile = "multiplicityMC_2M.root";
+const char* correctionFile = "multiplicity.root";
 const char* measuredFile   = "multiplicityMC_1M_3.root";
-Int_t etaRange = 3;
+Int_t etaRange = 2;
 Int_t displayRange = 200; // axis range
 Int_t ratioRange = 151;   // range to calculate difference
 Int_t longDisplayRange = 200;