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 44198fa..89843e7 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 82fac72..6f283d7 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 891de0b..6878e71 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 9593f8e..160a946 100644 (file)
@@ -79,6 +79,31 @@ Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
 }
 
 //____________________________________________________________________
+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)
 {
   // Get the vertex from the ESD and returns it if the vertex is valid
@@ -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 55132c3..022b233 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 457ef11..a167463 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 fe6a0e3..ddc15a8 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 c5f46ba..0cf621a 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 ee612ec..f61fa83 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 f1537d3..fe2ebe9 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 bebc840..38475da 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 1590528..436a923 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 0d55f9a..8431281 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 0f9aede..3fc4c24 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 f96e561..0a8a862 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 c778f41..c9f2785 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 268a26d..d4a7518 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 1894423..24f6e24 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;