added correction for events with vertex but 0 tracks
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Jan 2008 17:14:07 +0000 (17:14 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Jan 2008 17:14:07 +0000 (17:14 +0000)
added material budget study

13 files changed:
PWG0/AliCorrection.cxx
PWG0/AliCorrectionMatrix.cxx
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.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/esdTrackCuts/AliCutTask.cxx
PWG0/esdTrackCuts/draw.C
PWG0/esdTrackCuts/printStats.C

index bd551d5..3e909ed 100644 (file)
@@ -368,14 +368,31 @@ void AliCorrection::PrintInfo(Float_t ptCut)
   TH2* measuredEvents = GetEventCorrection()->GetMeasuredHistogram();
   TH2* generatedEvents = GetEventCorrection()->GetGeneratedHistogram();
 
+  Printf("AliCorrection::PrintInfo: Whole phasespace:");
+
   Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %1.f, events generated %1.f", measured->Integral(), generated->Integral(), measuredEvents->Integral(), generatedEvents->Integral());
 
-  // normalize to number of events;
-  measured->Scale(1.0 / measuredEvents->Integral());
-  generated->Scale(1.0 / generatedEvents->Integral());
+  Printf("AliCorrection::PrintInfo: Values in |eta| < 0.8, |vtx-z| < 10 and pt > %.2f:", ptCut);
+
+  Float_t tracksM = measured->Integral(measured->GetXaxis()->FindBin(-9.9), measured->GetXaxis()->FindBin(9.9), measured->GetYaxis()->FindBin(-0.79), measured->GetYaxis()->FindBin(0.79), measured->GetZaxis()->FindBin(ptCut), measured->GetZaxis()->GetNbins());
+  Float_t tracksG = generated->Integral(generated->GetXaxis()->FindBin(-9.9), generated->GetXaxis()->FindBin(9.9), generated->GetYaxis()->FindBin(-0.79), generated->GetYaxis()->FindBin(0.79), generated->GetZaxis()->FindBin(ptCut), generated->GetZaxis()->GetNbins());
+
+  Float_t eventsM = measuredEvents->Integral(measuredEvents->GetXaxis()->FindBin(-9.9), measuredEvents->GetXaxis()->FindBin(9.9), 1, measuredEvents->GetNbinsY());
+  Float_t eventsG = generatedEvents->Integral(generatedEvents->GetXaxis()->FindBin(-9.9), generatedEvents->GetXaxis()->FindBin(9.9), 1, generatedEvents->GetNbinsY());
+
+  Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %1.f, events generated %1.f", tracksM, tracksG, eventsM, eventsG);
 
-  Float_t nMeasured = measured->Integral(-1, -1, -1, -1, measured->GetZaxis()->FindBin(ptCut), measured->GetZaxis()->GetNbins());
-  Float_t nGenerated = generated->Integral(-1, -1, -1, -1, generated->GetZaxis()->FindBin(ptCut), generated->GetZaxis()->GetNbins());
+  if (tracksM > 0)
+    Printf("Effective average correction factor for TRACKS: %.3f", tracksG / tracksM);
+  if (eventsM > 0)
+  Printf("Effective average correction factor for EVENTS: %.3f", eventsG / eventsM);
 
-  Printf("%.2f tracks/event measured, %.2f tracks after correction --> effective average correction factor is %.2f (pt cut %.2f GeV/c)", nMeasured, nGenerated, nGenerated / nMeasured, ptCut);
+  if (eventsM > 0 && eventsG > 0)
+  {
+    // normalize to number of events;
+    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);
+  }
 }
index 990a1d4..66b7048 100644 (file)
@@ -362,7 +362,10 @@ void AliCorrectionMatrix::SetCorrectionToUnity()
   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)
+      {
         fhCorr->SetBinContent(x, y, z, 1);
+        fhCorr->SetBinError(x, y, z, 0);
+      }
 }
 
 //____________________________________________________________________
index 261d339..0e54435 100644 (file)
@@ -354,6 +354,7 @@ Float_t AlidNdEtaCorrection::GetMeasuredFraction(CorrectionType correctionType,
     etaEnd = etaBegin;
   }
 
+  //Printf("AlidNdEtaCorrection::GetMeasuredFraction: Using vtx range of +- 10 cm");
   Int_t vertexBegin = generated->GetXaxis()->FindBin(-9.99);
   Int_t vertexEnd = generated->GetXaxis()->FindBin(9.99);
 
@@ -387,6 +388,26 @@ Float_t AlidNdEtaCorrection::GetMeasuredFraction(CorrectionType correctionType,
   return fraction;
 }
 
+//____________________________________________________________________
+TH1* AlidNdEtaCorrection::GetMeasuredEventFraction(CorrectionType correctionType, Int_t multCut)
+{
+  // calculates the fraction of events above multCut (but including it)
+  //
+  // uses the generated event histogram from the correction passed, e.g. pass GetTrack2ParticleCorrection()
+
+  if (!GetCorrection(correctionType))
+    return 0;
+
+  const TH2F* generated = GetCorrection(correctionType)->GetEventCorrection()->GetGeneratedHistogram();
+
+  TH1* allEvents = generated->ProjectionX(Form("%s_all", generated->GetName()), 1, generated->GetNbinsY());
+  TH1* aboveEvents = generated->ProjectionX(Form("%s_above", generated->GetName()), generated->GetYaxis()->FindBin(multCut), generated->GetNbinsY());
+  
+  aboveEvents->Divide(aboveEvents, allEvents, 1, 1, "B");
+
+  return aboveEvents;  
+}
+
 //____________________________________________________________________
 void AlidNdEtaCorrection::ReduceInformation()
 {
index 6a9a704..c56c95f 100644 (file)
@@ -23,6 +23,7 @@
 #include "AliPWG0Helper.h"
 
 class AliCorrection;
+class TH1;
 
 class AlidNdEtaCorrection : public TNamed
 {
@@ -65,6 +66,7 @@ public:
   void    DrawOverview(const char* canvasName = 0);
 
   Float_t GetMeasuredFraction(CorrectionType correctionType, Float_t ptCutOff, Float_t eta = -100, Bool_t debug = kFALSE);
+  TH1*    GetMeasuredEventFraction(CorrectionType correctionType, Int_t multCut);
 
   void ReduceInformation();
 
index b693bfb..6c32b96 100644 (file)
@@ -33,6 +33,7 @@ dNdEtaAnalysis::dNdEtaAnalysis() :
 
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
+    fdNdEtaNotEventCorrected[i] = 0;
     fdNdEta[i] = 0;
     fdNdEtaPtCutOffCorrected[i] = 0;
   }
@@ -57,11 +58,13 @@ dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title, AliPWG0Helper::Analy
 
   fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", 40, -2, 2);
 
+  fdNdEtaNotEventCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_noteventcorrected", fdNdEta[0]->GetName())));
   fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
 
   for (Int_t i=1; i<kVertexBinning; ++i)
   {
     fdNdEta[i]    = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
+    fdNdEtaNotEventCorrected[i]    = dynamic_cast<TH1F*> (fdNdEtaNotEventCorrected[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
     fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
   }
 
@@ -83,6 +86,11 @@ dNdEtaAnalysis::~dNdEtaAnalysis()
 
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
+    if (fdNdEtaNotEventCorrected[i])
+    {
+      delete fdNdEtaNotEventCorrected[i];
+      fdNdEtaNotEventCorrected[i] = 0;
+    }
     if (fdNdEta[i])
     {
       delete fdNdEta[i];
@@ -139,6 +147,7 @@ void dNdEtaAnalysis::Copy(TObject &c) const
 
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
+    target.fdNdEtaNotEventCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaNotEventCorrected[i]->Clone());
     target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
     target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
   }
@@ -165,15 +174,22 @@ void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
 }
 
 //____________________________________________________________________
-void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType)
+void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, Int_t multCut)
 {
   //
   // correct with the given correction values and calculate dNdEta and pT distribution
   // the corrections that are applied can be steered by the flag correctionType
+  // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
   //
 
   // TODO put tag somewhere which corrections have been applied
 
+  if (multCut > 1)
+  {
+    Printf("ERROR: A bigger multiplicity cut than 1 is not possible in the current implementation");
+    return;
+  }
+
   // set corrections to 1
   fData->SetCorrectionToUnity();
 
@@ -218,12 +234,26 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     printf("INFO: No correction applied\n");
 
   fData->Multiply();
+  fData->PrintInfo(ptCut);
 
   TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
 
   // integrate multiplicity axis out (include under/overflow bins!!!)
   TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
-  TH1D* vertexHist = tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
+
+  // multiplicity cut correction
+  // correct for not-measured events with too small multiplicity
+  // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
+  TH1D* vertexHistUncorrected = tmp->ProjectionX("_px", tmp->GetYaxis()->FindBin(multCut), tmp->GetNbinsY() + 1, "e");
+  TH1D* vertexHist = (TH1D*) vertexHistUncorrected->Clone("vertexHist");
+
+  Printf("Using %d events (above a cut of %d)", (Int_t) vertexHistUncorrected->Integral(), multCut);
+  if (correction && correctionType >= AlidNdEtaCorrection::kVertexReco) 
+  {
+    TH1* eventFraction = correction->GetMeasuredEventFraction(correctionType, multCut);
+    vertexHist->Divide(eventFraction);
+    Printf("Multiplicity cut correction: Corrected from %d events to %d events", (Int_t) vertexHistUncorrected->Integral(), (Int_t) vertexHist->Integral());
+  }
 
   // create pt hist
   {
@@ -309,6 +339,13 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
 
       //printf("vertexBinBegin = %d, vertexBinEnd = %d\n", vertexBinBegin, vertexBinEnd);
 
+      Float_t totalEventsUncorrected = vertexHistUncorrected->Integral(vertexBinBegin, vertexBinEnd - 1);
+      if (totalEventsUncorrected == 0)
+      {
+        printf("WARNING: No events (uncorrected) for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
+        continue;
+      }
+
       Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd - 1);
       if (totalEvents == 0)
       {
@@ -346,8 +383,17 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
       if (bin > 0 && bin < fdNdEta[vertexPos]->GetNbinsX())
       {
-        Float_t dndeta = sum / totalEvents;
-        Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
+        Float_t dndeta = sum / totalEventsUncorrected;
+        Float_t error  = TMath::Sqrt(sumError2) / totalEventsUncorrected;
+
+        dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
+        error  = error / fdNdEta[vertexPos]->GetBinWidth(bin);
+
+        fdNdEtaNotEventCorrected[vertexPos]->SetBinContent(bin, dndeta);
+        fdNdEtaNotEventCorrected[vertexPos]->SetBinError(bin, error);
+  
+        dndeta = sum / totalEvents;
+        error  = TMath::Sqrt(sumError2) / totalEvents;
 
         dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
         error  = error / fdNdEta[vertexPos]->GetBinWidth(bin);
@@ -391,6 +437,11 @@ void dNdEtaAnalysis::SaveHistograms()
 
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
+    if (fdNdEtaNotEventCorrected[i])
+      fdNdEtaNotEventCorrected[i]->Write();
+    else
+      printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaNotEventCorrected[%d] is 0\n", i);
+
     if (fdNdEta[i])
       fdNdEta[i]->Write();
     else
@@ -418,6 +469,9 @@ void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
 
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
+    if (dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaNotEventCorrected[i]->GetName())))
+      fdNdEtaNotEventCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaNotEventCorrected[i]->GetName()));
+  
     fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
   }
@@ -540,7 +594,7 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list)
   TObject* obj;
 
   // sub collections
-  const Int_t nCollections = 2 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
+  const Int_t nCollections = 3 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
   TList* collections[nCollections];
   for (Int_t i=0; i<nCollections; ++i)
     collections[i] = new TList;
@@ -559,6 +613,7 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list)
     {
       collections[2+i]->Add(entry->fdNdEta[i]);
       collections[2+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
+      collections[2+2*kVertexBinning+i]->Add(entry->fdNdEtaNotEventCorrected[i]);
     }
 
     ++count;
@@ -569,7 +624,8 @@ Long64_t dNdEtaAnalysis::Merge(TCollection* list)
   for (Int_t i=0; i<kVertexBinning; ++i)
   {
     fdNdEta[i]->Merge(collections[2+i]);
-    fdNdEta[i]->Merge(collections[2+kVertexBinning+i]);
+    fdNdEtaPtCutOffCorrected[i]->Merge(collections[2+kVertexBinning+i]);
+    fdNdEtaNotEventCorrected[i]->Merge(collections[2+kVertexBinning+i]);
   }
 
   for (Int_t i=0; i<nCollections; ++i)
index 20d0e3d..1ea6462 100644 (file)
@@ -42,7 +42,7 @@ public:
   void FillTrack(Float_t vtx, Float_t eta, Float_t pt);
   void FillEvent(Float_t vtx, Float_t n);
 
-  void Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType);
+  void Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, Int_t multCut = 0);
 
   void DrawHistograms(Bool_t simple = kFALSE);
   void LoadHistograms(const Char_t* dir = 0);
@@ -62,8 +62,9 @@ protected:
 
   TH1F* fPtDist; // pt distribution
 
-  TH1F* fdNdEta[kVertexBinning]; // dndeta results for different vertex bins (0 = full range)
-  TH1F* fdNdEtaPtCutOffCorrected[kVertexBinning];  // dndeta results for different vertex bins (0 = full range), pt cut off corrected
+  TH1F* fdNdEtaNotEventCorrected[kVertexBinning];  // dndeta results for different vertex bins (0 = full range)
+  TH1F* fdNdEta[kVertexBinning];                   // dndeta results for different vertex bins (0 = full range), mult cut off corrected
+  TH1F* fdNdEtaPtCutOffCorrected[kVertexBinning];  // dndeta results for different vertex bins (0 = full range), mult + pt cut off corrected 
 
   ClassDef(dNdEtaAnalysis, 1)
 };
index 5c435d7..4e1f19f 100644 (file)
@@ -102,7 +102,7 @@ void InitPadCOLZ()
 
 // --- end of helpers --- begin functions ---
 
-void DrawOverview(const char* fileName, const char* dirName)
+void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
 {
   loadlibs();
   if (!TFile::Open(fileName))
@@ -247,7 +247,7 @@ Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc =
   return diffFullRange;
 }
 
-void dNdEta(Bool_t onlyESD = kFALSE)
+void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
 {
   TFile* file = TFile::Open("analysis_esd.root");
   TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
@@ -258,6 +258,8 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
   TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
 
+  TH1* histESDNoEvCorr = (TH1*) file->Get("dndeta/dNdEta_noteventcorrected");
+
   Prepare1DPlot(histESD);
   Prepare1DPlot(histESDMB);
   Prepare1DPlot(histESDMBVtx);
@@ -267,6 +269,8 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   Prepare1DPlot(histESDMBVtxNoPt);
   Prepare1DPlot(histESDMBTracksNoPt);
 
+  Prepare1DPlot(histESDNoEvCorr);
+  
   histESD->SetLineWidth(0);
   histESDMB->SetLineWidth(0);
   histESDMBVtx->SetLineWidth(0);
@@ -279,11 +283,17 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMB->SetMarkerColor(2);
   histESDMBVtx->SetMarkerColor(3);
 
+  histESD->SetLineColor(1);
+  histESDMB->SetLineColor(2);
+  histESDMBVtx->SetLineColor(3);
+
   histESDNoPt->SetMarkerColor(1);
   histESDMBNoPt->SetMarkerColor(2);
   histESDMBVtxNoPt->SetMarkerColor(3);
   histESDMBTracksNoPt->SetMarkerColor(4);
 
+  histESDNoEvCorr->SetMarkerColor(6);
+
   histESD->SetMarkerStyle(20);
   histESDMB->SetMarkerStyle(21);
   histESDMBVtx->SetMarkerStyle(22);
@@ -292,6 +302,8 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMBNoPt->SetMarkerStyle(21);
   histESDMBVtxNoPt->SetMarkerStyle(22);
   histESDMBTracksNoPt->SetMarkerStyle(23);
+  
+  histESDNoEvCorr->SetMarkerStyle(29);
 
   TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 0, histESDMBVtx->GetMaximum() * 1.1);
   Prepare1DPlot(dummy);
@@ -300,7 +312,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   dummy->SetYTitle("dN_{ch}/d#eta");
   dummy->GetYaxis()->SetTitleOffset(1);
 
-  Float_t etaLimit = 0.7999;
+  Float_t etaLimit = 1.2999;
 
   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
@@ -318,8 +330,11 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMB->Draw("SAME");
   histESD->Draw("SAME");
 
-  canvas->SaveAs("dNdEta1.gif");
-  canvas->SaveAs("dNdEta1.eps");
+  if (save)
+  {
+    canvas->SaveAs("dNdEta1.gif");
+    canvas->SaveAs("dNdEta1.eps");
+  }
 
   if (onlyESD)
     return;
@@ -333,8 +348,13 @@ void dNdEta(Bool_t onlyESD = kFALSE)
 
   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
   fdNdEtaAnalysis->LoadHistograms();
+  dNdEtaAnalysis* fdNdEtaAnalysis2 = (dNdEtaAnalysis*) fdNdEtaAnalysis->Clone();
+
   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
-  TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
+  TH1* histMCPtCut = (TH1*) fdNdEtaAnalysis->GetdNdEtaHistogram(0)->Clone("histMCPtCut");
+
+  fdNdEtaAnalysis2->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
+  TH1* histMCPtCutNoEvCorr = (TH1*) fdNdEtaAnalysis2->GetdNdEtaHistogram(0)->Clone("histMCPtCutNoEvCorr");
 
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
   fdNdEtaAnalysis->LoadHistograms();
@@ -343,12 +363,12 @@ void dNdEta(Bool_t onlyESD = kFALSE)
 
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
   fdNdEtaAnalysis->LoadHistograms();
-  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
   TH1* histMCTrVtxPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
   fdNdEtaAnalysis->LoadHistograms();
-  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
   TH1* histMCTracksPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
 
   TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
@@ -358,6 +378,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   Prepare1DPlot(histMCTrVtx);
 
   Prepare1DPlot(histMCPtCut);
+  Prepare1DPlot(histMCPtCutNoEvCorr);
   Prepare1DPlot(histMCTrPtCut);
   Prepare1DPlot(histMCTrVtxPtCut);
   if (histMCTracksPtCut)
@@ -372,6 +393,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histMCTrVtxPtCut->SetLineColor(3);
   if (histMCTracksPtCut)
     histMCTracksPtCut->SetLineColor(4);
+  histMCPtCutNoEvCorr->SetLineColor(6);
 
   TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
   Prepare1DPlot(dummy2);
@@ -388,14 +410,19 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMBNoPt->Draw("SAME");
   histESDMBVtxNoPt->Draw("SAME");
   histESDMBTracksNoPt->Draw("SAME");
+  histESDNoEvCorr->Draw("SAME");
   histMCPtCut->Draw("SAME");
+  histMCPtCutNoEvCorr->Draw("SAME");
   histMCTrPtCut->Draw("SAME");
   histMCTrVtxPtCut->Draw("SAME");
   if (histMCTracksPtCut)
     histMCTracksPtCut->Draw("SAME");
 
-  canvas2->SaveAs("dNdEta2.gif");
-  canvas2->SaveAs("dNdEta2.eps");
+  if (save)
+  {
+    canvas2->SaveAs("dNdEta2.gif");
+    canvas2->SaveAs("dNdEta2.eps");
+  }
 
   TH1* ratio = (TH1*) histMC->Clone("ratio");
   TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
@@ -403,7 +430,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   ratio->Divide(histESD);
   ratioNoPt->Divide(histESDNoPt);
 
-  ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
+  ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 
   ratio->SetLineColor(1);
   ratioNoPt->SetLineColor(2);
@@ -486,8 +513,11 @@ void dNdEta(Bool_t onlyESD = kFALSE)
 
   canvas3->Modified();
 
-  canvas3->SaveAs("dNdEta.gif");
-  canvas3->SaveAs("dNdEta.eps");
+  if (save)
+  {
+    canvas3->SaveAs("dNdEta.gif");
+    canvas3->SaveAs("dNdEta.eps");
+  }
 
   TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 
@@ -1444,15 +1474,85 @@ void CompareCorrection2Measured(const char* dataInput = "analysis_esd_raw.root",
   
   TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
   hist1->SetTitle("mc");
+  Printf("mc contains %f entries", hist1->Integral());
+  Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), 1, hist1->GetNbinsY(), hist1->GetZaxis()->FindBin(0.301), hist1->GetNbinsZ()));
 
   TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
   hist2->SetTitle("esd");
+  Printf("esd contains %f entries", hist2->Integral());
+  Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), 1, hist2->GetNbinsY(), hist2->GetZaxis()->FindBin(0.301), hist2->GetNbinsZ()));
 
   AliPWG0Helper::CreateDividedProjections(hist1, hist2);
 
   new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
   new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
   new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
+}
+
+void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
+{
+  loadlibs();
+
+  TFile* file = TFile::Open(dataInput);
+
+  if (!file)
+  {
+    cout << "Error. File not found" << endl;
+    return;
+  }
+
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+
+  TFile* file = TFile::Open(dataInput2);
+
+  if (!file)
+  {
+    cout << "Error. File not found" << endl;
+    return;
+  }
+
+  dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
+  fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
+
+  gROOT->cd();
+
+  TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
+  hist1->SetTitle("esd1");
+  Printf("esd1 contains %f entries", hist1->GetEntries());
+  Printf("esd1 contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), 1, hist1->GetNbinsY(), hist1->GetZaxis()->FindBin(0.301), hist1->GetNbinsZ()));
+
+  TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
+  hist2->SetTitle("esd2");
+  Printf("esd2 contains %f entries", hist2->GetEntries());
+  Printf("esd2 contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), 1, hist2->GetNbinsY(), hist2->GetZaxis()->FindBin(0.301), hist2->GetNbinsZ()));
+
+  AliPWG0Helper::CreateDividedProjections(hist1, hist2);
+
+  new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
+  new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
+  new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
+
+  TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
+  TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
+
+  Printf("event1 contains %f entries", event1->GetEntries());
+  Printf("event2 contains %f entries", event2->GetEntries());
+  Printf("event1 integral is %f", event1->Integral());
+  Printf("event2 integral is %f", event2->Integral());
+  Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
+  Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
+
+  projx1 = event1->ProjectionX();
+  projx2 = event2->ProjectionX();
+
+  new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
+
+  projx1->Divide(projx2);
+  new TCanvas; projx1->Draw();
+
+  event1->Divide(event2);
+  new TCanvas; event1->Draw("COLZ");
 
 }
-  
+
index 691fbb7..97dd60e 100644 (file)
@@ -676,187 +676,6 @@ void Composition()
   Track2Particle1DComposition(fileNames, nCompositions, folderNames);
 }
 
-Double_t ConvSigma1To2D(Double_t sigma)
-{
-  return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
-}
-
-Double_t ConvDistance1DTo2D(Double_t distance)
-{
-  return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
-}
-
-Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
-{
-  Double_t count = 0;
-
-  //nSigma = ConvSigma1To2D(nSigma);
-
-  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
-    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
-    {
-      Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
-      Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
-
-      Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
-
-      d = ConvDistance1DTo2D(d);
-
-      if (d < nSigma)
-        count += tracks->GetBinContent(x, y);
-    }
-
-  return count;
-}
-
-TH2F* Sigma2VertexGaussianTracksHist()
-{
-  TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
-
-  TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
-  gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
-
-  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
-    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
-      tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
-
-  //normalize
-  tracks->Scale(1.0 / tracks->Integral());
-
-  return tracks;
-}
-
-TH1F* Sigma2VertexGaussian()
-{
-  TH2F* tracks = Sigma2VertexGaussianTracksHist();
-
-  TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
-  canvas->Divide(2, 2);
-
-  canvas->cd(1);
-  tracks->Draw("COLZ");
-
-  TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
-  for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
-    ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
-  ratio->SetMarkerStyle(21);
-
-  canvas->cd(2);
-  ratio->DrawCopy("P");
-
-  TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 4 sigma / % included n sigma", 50, 0.05, 5.05);
-  Double_t sigma3 = Sigma2VertexCount(tracks, 4);
-  for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
-    ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
-  ratio2->SetMarkerStyle(21);
-
-  canvas->cd(3);
-  ratio2->DrawCopy("P");
-
-  canvas->SaveAs("Sigma2Vertex.eps");
-
-  return ratio2;
-}
-
-TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
-{
-  TFile* file = TFile::Open(fileName);
-
-  TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
-  TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
-  if (!sigmavertex || !sigmavertexPrim)
-  {
-    printf("Could not read histogram(s)\n");
-    return;
-  }
-
-  // calculate ratio
-  TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 4 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
-
-  // calculate contamination
-  TH1F* contamination = ratio->Clone("sigmavertexsimulation_contamination");
-  contamination->SetTitle("sigmavertexsimulation_contamination;N#sigma;1 + N_{secondaries} / N_{all}");
-
-  for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
-  {
-    ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(4)) / sigmavertex->GetBinContent(i));
-    contamination->SetBinContent(i, 1 + (sigmavertex->GetBinContent(i) - sigmavertexPrim->GetBinContent(i)) / sigmavertex->GetBinContent(i));
-  }
-
-  // print stats
-  for (Float_t sigma = 2.0; sigma < 5.25; sigma += 0.5)
-  {
-    Float_t error1 = 1 - ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma - 0.5));
-    Float_t error2 = -1 + ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma + 0.5));
-    Float_t cont = -1 + contamination->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma));
-    Printf("%.2f sigma --> syst. error = - %.2f %% + %.2f %%, cont. = %.2f %%", sigma, error1 * 100, error2 * 100, cont * 100);
-  }
-
-  TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
-  canvas->Divide(2, 1);
-
-  canvas->cd(1);
-  sigmavertex->SetMarkerStyle(21);
-  sigmavertex->Draw("P");
-
-  canvas->cd(2);
-  ratio->SetMarkerStyle(21);
-  ratio->DrawCopy("P");
-
-  contamination->DrawCopy("SAME");
-
-  TH1F** returnContainer = new TH1F*[2];
-  returnContainer[0] = ratio;
-  returnContainer[1] = contamination;
-
-  return returnContainer;
-}
-
-void Sigma2VertexCompare(const char* fileName = "systematics.root")
-{
-  TH1F* ratio1 = Sigma2VertexGaussian();
-
-  TH1F** hists = Sigma2VertexSimulation(fileName);
-  TH1F* ratio2 = hists[0];
-  TH1F* contamination = hists[1];
-
-  ratio1->SetStats(kFALSE);
-  ratio2->SetStats(kFALSE);
-
-  ratio1->SetMarkerStyle(0);
-  ratio2->SetMarkerStyle(0);
-
-  ratio1->SetLineWidth(2);
-  ratio2->SetLineWidth(2);
-
-  TLegend* legend = new TLegend(0.6, 0.8, 0.95, 0.95);
-  legend->SetFillColor(0);
-  legend->AddEntry(ratio1, "Gaussian");
-  legend->AddEntry(ratio2, "Simulation");
-  legend->AddEntry(contamination, "1 + Contamination");
-
-  ratio2->SetTitle("");
-  ratio2->GetYaxis()->SetTitleOffset(1.5);
-  ratio2->GetXaxis()->SetRangeUser(2, 5);
-
-  TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
-  InitPad();
-
-  ratio2->SetMarkerStyle(21);
-  ratio1->SetMarkerStyle(22);
-
-  ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
-  ratio2->SetLineColor(kRed);
-  ratio2->SetMarkerColor(kRed);
-  ratio2->Draw("PL");
-  ratio1->Draw("SAMEPL");
-
-  contamination->Draw("SAME");
-
-  legend->Draw();
-
-  canvas->SaveAs("Sigma2VertexCompare.eps");
-}
 
 void drawSystematics()
 {
index de420b0..0e5cb81 100644 (file)
@@ -12,6 +12,21 @@ void loadlibs()
   gSystem->Load("libPWG0base");
 }
 
+void InitPad()
+{
+  if (!gPad)
+    return;
+
+  gPad->Range(0, 0, 1, 1);
+  gPad->SetLeftMargin(0.15);
+  //gPad->SetRightMargin(0.05);
+  //gPad->SetTopMargin(0.13);
+  //gPad->SetBottomMargin(0.1);
+
+  gPad->SetGridx();
+  gPad->SetGridy();
+}
+
 void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
 {
   //gROOT->ProcessLine(".L drawPlots.C");
@@ -22,8 +37,6 @@ void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
   const char* legendNames[] = { "#pi", "K", "p", "standard" };
   Int_t folderCount = 3;
 
-
-
   TH2F* h2DCorrections[4];
   TH1F* h1DCorrections[4];
   for (Int_t i=0; i<4; i++) {
@@ -275,10 +288,13 @@ void MisalignmentShowRawTrackPlots(const char* dirName = "fdNdEtaAnalysisESD")
 
 void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files, const char** dirs, const char** names, Int_t* histID)
 {
-  gSystem->Load("libPWG0base");
+  loadlibs();
+  gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C");
 
   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
   canvas->Divide(2, 1);
+  canvas->cd(2)->SetGridx();
+  canvas->cd(2)->SetGridy();
 
   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
   legend->SetFillColor(0);
@@ -307,9 +323,11 @@ void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files, const
 
     if (i != 0)
     {
+      PrintIntegratedDeviation(hist, base, names[i]);
+
       canvas->cd(2);
-      hist->Divide(hist, base, 1, 1, "B");
-      hist->GetYaxis()->SetRangeUser(0.98, 1.02);
+      hist->Divide(hist, base, 1, 1);
+      hist->GetYaxis()->SetRangeUser(0.9, 1.1);
       hist->DrawCopy((i == 1) ? "" : "SAME");
     }
   }
@@ -370,6 +388,47 @@ void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files1, cons
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 }
 
+void MaterialBudgetChange()
+{
+  const char* files[] =
+    { "Material-normal-mcvtx/analysis_esd.root",
+      "Material-increased-mcvtx/analysis_esd.root",
+      "Material-decreased-mcvtx/analysis_esd.root" };
+
+  const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
+  const char* names[] = { "no change", "+ 10 % material", "- 10 % material" };
+  Int_t hist[] = { 0, 0, 0 };
+
+  drawdNdEtaRatios("MaterialBudgetChange", 3, files, dirs, names, hist);
+}
+
+void MisalignmentChange()
+{
+  const char* files[] =
+    { "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-sim-mcvtx/analysis_esd.root",
+      "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-simrec-mcvtx/analysis_esd.root" };
+
+  const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
+  const char* names[] = { "no change", "increased in sim/rec", "increased in sim" };
+  Int_t hist[] = { 0, 0, 0 };
+
+  drawdNdEtaRatios("MisalignmentChange", 2, files, dirs, names, hist);
+}
+
+void dNdEtaVertexRanges()
+{
+  const char* files[] =
+    { "analysis_esd.root",
+      "analysis_esd.root",
+      "analysis_esd.root" };
+
+  const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
+  const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
+  Int_t hist[] = { 0, 1, 2 };
+
+  drawdNdEtaRatios("dNdEtaVertexRanges", 3, files, dirs, names, hist);
+}
+
 void vertexShiftStudy(Int_t histID)
 {
   const char* files[] = { "maps/idealA/mc-vertex/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.05/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.1/analysis_esd.root", "maps/idealA/mc-vertex-shift-dep/analysis_esd.root" };
@@ -411,18 +470,17 @@ void vertexShift()
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 }
 
-void CompareRawTrackPlots(Float_t ptCut = 0.0)
+void CompareRawTrackPlots(const char* fileName1, const char* fileName2, Float_t ptCut = 0.0)
 {
   loadlibs();
 
   const char* dirName = "fdNdEtaAnalysisESD";
 
-  TFile* file = TFile::Open("fullA-simrec/MB2/analysis_esd_raw.root");
+  TFile* file = TFile::Open(fileName1);
   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
   fdNdEtaAnalysis->LoadHistograms();
 
-  //TFile* file2 = TFile::Open("fullA-sim/analysis_esd_raw.root");
-  TFile* file2 = TFile::Open("fullA-simrecexcepttpc/analysis_esd_raw.root");
+  TFile* file2 = TFile::Open(fileName2);
   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
   fdNdEtaAnalysis2->LoadHistograms();
 
@@ -437,8 +495,8 @@ void CompareRawTrackPlots(Float_t ptCut = 0.0)
   track1->Scale(1.0 / event1Count);
   track2->Scale(1.0 / event2Count);
 
-  Float_t nTrack1 = track1->Integral(-1, -1, -1, -1, track1->GetZaxis()->FindBin(ptCut), track1->GetZaxis()->GetNbins());
-  Float_t nTrack2 = track2->Integral(-1, -1, -1, -1, track2->GetZaxis()->FindBin(ptCut), track2->GetZaxis()->GetNbins());
+  Float_t nTrack1 = track1->Integral(track1->GetXaxis()->FindBin(-9.9), track1->GetXaxis()->FindBin(9.9), track1->GetYaxis()->FindBin(-0.79), track1->GetYaxis()->FindBin(0.79), track1->GetZaxis()->FindBin(ptCut), track1->GetZaxis()->GetNbins());
+  Float_t nTrack2 = track2->Integral(track2->GetXaxis()->FindBin(-9.9), track2->GetXaxis()->FindBin(9.9), track2->GetYaxis()->FindBin(-0.79), track2->GetYaxis()->FindBin(0.79), track2->GetZaxis()->FindBin(ptCut), track2->GetZaxis()->GetNbins());
 
   Printf("There are %.2f tracks/event in the first sample and %.2f tracks/event in the second sample. %.2f %% difference (with pt cut at %.2f GeV/c)", nTrack1, nTrack2, 100.0 * (nTrack1 - nTrack2) / nTrack1, ptCut);
 
@@ -450,20 +508,26 @@ void CompareRawTrackPlots(Float_t ptCut = 0.0)
   new TCanvas; gROOT->FindObject("track1_zx_div_track2_zx")->Draw("COLZ");
   new TCanvas; gROOT->FindObject("track1_zy_div_track2_zy")->Draw("COLZ");
 
-  new TCanvas;
-  proj1 = track1->Project3D("ze2");
-  proj2 = track2->Project3D("ze2");
-  AliPWG0Helper::NormalizeToBinWidth(proj1);
-  AliPWG0Helper::NormalizeToBinWidth(proj2);
-
-  proj1->DrawCopy();
-  proj2->SetLineColor(2);
-  proj2->SetMarkerColor(2);
-  proj2->DrawCopy("SAME");
-
-  AliPWG0Helper::CreateDividedProjections(track1, track2, "ze2");
-  TH1* pt = gROOT->FindObject("track1_ze2_div_track2_ze2");
-  new TCanvas; pt->Draw();
+  for (Int_t i=0; i<3; i++)
+  {
+    char c = 'x' + (char) i;
+
+    /*proj1 = track1->Project3D(Form("%ce2", c));
+    proj2 = track2->Project3D(Form("%ce2", c));
+    AliPWG0Helper::NormalizeToBinWidth(proj1);
+    AliPWG0Helper::NormalizeToBinWidth(proj2);
+
+    new TCanvas;
+    proj1->DrawCopy();
+    proj2->SetLineColor(2);
+    proj2->SetMarkerColor(2);
+    proj2->DrawCopy("SAME");*/
+
+    AliPWG0Helper::CreateDividedProjections(track1, track2, Form("%ce2", c));
+    TH1* pt = gROOT->FindObject(Form("track1_%ce2_div_track2_%ce2", c, c));
+    new TCanvas; pt->DrawCopy();
+    gPad->SetGridx(); gPad->SetGridy();
+  }
 }
 
 void MagnitudeOfCorrection(const char* fileName, const char* dirName = "dndeta", Float_t ptCut = 0.3)
@@ -476,3 +540,184 @@ void MagnitudeOfCorrection(const char* fileName, const char* dirName = "dndeta",
   fdNdEtaAnalysis->GetData()->PrintInfo(ptCut);
 }
 
+Double_t ConvSigma1To2D(Double_t sigma)
+{
+  return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
+}
+
+Double_t ConvDistance1DTo2D(Double_t distance)
+{
+  return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
+}
+
+Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
+{
+  Double_t count = 0;
+
+  //nSigma = ConvSigma1To2D(nSigma);
+
+  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
+    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
+    {
+      Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
+      Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
+
+      Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
+
+      d = ConvDistance1DTo2D(d);
+
+      if (d < nSigma)
+        count += tracks->GetBinContent(x, y);
+    }
+
+  return count;
+}
+
+TH2F* Sigma2VertexGaussianTracksHist()
+{
+  TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
+
+  TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
+  gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
+
+  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
+    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
+      tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
+
+  //normalize
+  tracks->Scale(1.0 / tracks->Integral());
+
+  return tracks;
+}
+
+TH1F* Sigma2VertexGaussian()
+{
+  TH2F* tracks = Sigma2VertexGaussianTracksHist();
+
+  TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
+  canvas->Divide(2, 2);
+
+  canvas->cd(1);
+  tracks->Draw("COLZ");
+
+  TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
+  for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
+    ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
+  ratio->SetMarkerStyle(21);
+
+  canvas->cd(2);
+  ratio->DrawCopy("P");
+
+  TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 4 sigma / % included n sigma", 50, 0.05, 5.05);
+  Double_t sigma3 = Sigma2VertexCount(tracks, 4);
+  for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
+    ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
+  ratio2->SetMarkerStyle(21);
+
+  canvas->cd(3);
+  ratio2->DrawCopy("P");
+
+  canvas->SaveAs("Sigma2Vertex.eps");
+
+  return ratio2;
+}
+
+TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
+{
+  TFile* file = TFile::Open(fileName);
+
+  TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
+  TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
+  if (!sigmavertex || !sigmavertexPrim)
+  {
+    printf("Could not read histogram(s)\n");
+    return;
+  }
+
+  // calculate ratio
+  TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 4 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
+
+  // calculate contamination
+  TH1F* contamination = ratio->Clone("sigmavertexsimulation_contamination");
+  contamination->SetTitle("sigmavertexsimulation_contamination;N#sigma;1 + N_{secondaries} / N_{all}");
+
+  for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
+  {
+    ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(4)) / sigmavertex->GetBinContent(i));
+    contamination->SetBinContent(i, 1 + (sigmavertex->GetBinContent(i) - sigmavertexPrim->GetBinContent(i)) / sigmavertex->GetBinContent(i));
+  }
+
+  // print stats
+  for (Float_t sigma = 2.0; sigma < 5.25; sigma += 0.5)
+  {
+    Float_t error1 = 1 - ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma - 0.5));
+    Float_t error2 = -1 + ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma + 0.5));
+    Float_t cont = -1 + contamination->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma));
+    Printf("%.2f sigma --> syst. error = - %.2f %% + %.2f %%, cont. = %.2f %%", sigma, error1 * 100, error2 * 100, cont * 100);
+  }
+
+  TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
+  canvas->Divide(2, 1);
+
+  canvas->cd(1);
+  sigmavertex->SetMarkerStyle(21);
+  sigmavertex->Draw("P");
+
+  canvas->cd(2);
+  ratio->SetMarkerStyle(21);
+  ratio->DrawCopy("P");
+
+  contamination->DrawCopy("SAME");
+
+  TH1F** returnContainer = new TH1F*[2];
+  returnContainer[0] = ratio;
+  returnContainer[1] = contamination;
+
+  return returnContainer;
+}
+
+void Sigma2VertexCompare(const char* fileName = "systematics.root")
+{
+  TH1F* ratio1 = Sigma2VertexGaussian();
+
+  TH1F** hists = Sigma2VertexSimulation(fileName);
+  TH1F* ratio2 = hists[0];
+  TH1F* contamination = hists[1];
+
+  ratio1->SetStats(kFALSE);
+  ratio2->SetStats(kFALSE);
+
+  ratio1->SetMarkerStyle(0);
+  ratio2->SetMarkerStyle(0);
+
+  ratio1->SetLineWidth(2);
+  ratio2->SetLineWidth(2);
+
+  TLegend* legend = new TLegend(0.6, 0.8, 0.95, 0.95);
+  legend->SetFillColor(0);
+  legend->AddEntry(ratio1, "Gaussian");
+  legend->AddEntry(ratio2, "Simulation");
+  legend->AddEntry(contamination, "1 + Contamination");
+
+  ratio2->SetTitle("");
+  ratio2->GetYaxis()->SetTitleOffset(1.5);
+  ratio2->GetXaxis()->SetRangeUser(2, 5);
+
+  TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
+  InitPad();
+
+  ratio2->SetMarkerStyle(21);
+  ratio1->SetMarkerStyle(22);
+
+  ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
+  ratio2->SetLineColor(kRed);
+  ratio2->SetMarkerColor(kRed);
+  ratio2->Draw("PL");
+  ratio1->Draw("SAMEPL");
+
+  contamination->Draw("SAME");
+
+  legend->Draw();
+
+  canvas->SaveAs("Sigma2VertexCompare.eps");
+}
index 76db7e4..66754d5 100644 (file)
@@ -125,7 +125,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
 
   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL, 1);
   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
   TFile* file2 = TFile::Open(dataOutput, "RECREATE");
   fdNdEtaAnalysis->SaveHistograms();
@@ -133,7 +133,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   file->cd();
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco);
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco, 1);
   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
   file2->cd();
   fdNdEtaAnalysis->SaveHistograms();
@@ -141,7 +141,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   file->cd();
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle);
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, 1);
   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
   file2->cd();
   fdNdEtaAnalysis->SaveHistograms();
@@ -149,7 +149,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   file->cd();
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
   fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
   //fdNdEtaAnalysis->DrawHistograms(kTRUE);
   file2->cd();
   fdNdEtaAnalysis->SaveHistograms();
@@ -178,7 +178,7 @@ void* FinishAnalysis(const char* analysisFile = "analysis_esd_raw.root", const c
     //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
   }
   else
-    fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone);
+    fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
 
   fdNdEtaAnalysis->DrawHistograms(simple);
 
index b408fa4..26b0794 100644 (file)
@@ -212,18 +212,18 @@ void AliCutTask::Exec(Option_t *)
       Int_t label = TMath::Abs(track->GetLabel());
       if (stack->IsPhysicalPrimary(label) == kTRUE)
       {
-       if (label >= max)
-       {
-         Printf("Warning label %d is higher than number of primaries %d", label, max);
-         continue;
-       }
-           
-        if (fTrackCutsPrimaries->AcceptTrack(track)) 
-       {
-         primAcc[label]++;
-       } 
-       else
-         primRej[label]++;
+        if (label >= max)
+        {
+          Printf("Warning label %d is higher than number of primaries %d", label, max);
+          continue;
+        }
+
+        if (fTrackCutsPrimaries->AcceptTrack(track))
+        {
+          primAcc[label]++;
+        }
+        else
+          primRej[label]++;
 
       }
       else
@@ -253,20 +253,20 @@ void AliCutTask::Exec(Option_t *)
 
     for (Int_t i=0; i<max; i++) {
       if (primAcc[i] == 1) {
-       fPrimStats->Fill(1);
+        fPrimStats->Fill(1);
       } else if (primAcc[i] > 1) {
-       fPrimStats->Fill(2);
-       fPrimStats->Fill(3, primAcc[i]);
+        fPrimStats->Fill(2);
+        fPrimStats->Fill(3, primAcc[i]);
       }
 
       if (primRej[i] > 0) {
-       if (primAcc[i] == 0) {
-         fPrimStats->Fill(4);
-         fPrimStats->Fill(5, primRej[i]);
-       } else if (primAcc[i] > 0) {
-         fPrimStats->Fill(6);
-         fPrimStats->Fill(7, primRej[i]);
-       }
+        if (primAcc[i] == 0) {
+          fPrimStats->Fill(4);
+          fPrimStats->Fill(5, primRej[i]);
+        } else if (primAcc[i] > 0) {
+          fPrimStats->Fill(6);
+          fPrimStats->Fill(7, primRej[i]);
+        }
       }
     }
 
index d458a72..5f150b8 100644 (file)
@@ -1,22 +1,33 @@
 void draw()
 {
-  const char* files[] = { "trackCuts_normal.root", "trackCuts_increased.root", "trackCuts_decreased.root" };
+  //const char* files[] = { "trackCuts_normal.root", "trackCuts_increased.root", "trackCuts_decreased.root" };
+  const char* files[] =
+    { "Material-normal/trackCuts.root",
+      "Material-increased-mcvtx/trackCuts.root",
+      "Material-decreased-mcvtx/trackCuts.root" };
   const char* titles[] = { "default geometry", "+ 10% material", "- 10% material" };
   Int_t colors[] = { 1, 2, 4 };
+  Int_t markers[] = {24, 25, 26, 27 };
 
   TCanvas* c = new TCanvas;
+  TCanvas* c2 = new TCanvas;
 
-  TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
+  TLegend* legend = new TLegend(0.6, 0.5, 0.9, 0.9);
+  TLegend* legend2 = new TLegend(0.7, 0.7, 0.9, 0.9);
 
   for (Int_t i=0; i<3; i++) {
-    TFile::Open(files[i]);
-    
+    if (!TFile::Open(files[i]))
+      return;
+
     TH1* ptPrim = gFile->Get("fTrackCutsPrimaries/before_cuts/pt");
     TH1* ptPrimCut = gFile->Get("fTrackCutsPrimaries/after_cuts/pt_cut");
 
     TH1* ptSec = gFile->Get("fTrackCutsSecondaries/before_cuts/pt");
     TH1* ptSecCut = gFile->Get("fTrackCutsSecondaries/after_cuts/pt_cut");
 
+    TH1* vertex = gFile->Get("fVertex");
+    Int_t nEvents = vertex->GetEntries();
+
     ptPrim->Add(ptSec);
     ptPrimCut->Add(ptSecCut);
 
@@ -24,34 +35,83 @@ void draw()
     ptPrimCut->Sumw2();
     ptSec->Sumw2();
     ptSecCut->Sumw2();
-    
+
     Printf("%s", titles[i]);
+    Printf("Total particles: %d", (Int_t) (ptPrim->GetEntries() + ptSec->GetEntries()));
+    Printf("Total particles/event: %.2f", (ptPrim->GetEntries() + ptSec->GetEntries()) / nEvents);
+    Printf("Primaries/event: %.2f", ptPrim->GetEntries() / nEvents);
+    Printf("Secondaries/event: %.2f", ptSec->GetEntries() / nEvents);
+    Printf("Primaries > 0.2 GeV/c/event: %.2f", ptPrim->Integral(ptPrim->GetXaxis()->FindBin(0.21), ptPrim->GetNbinsX()) / nEvents);
+    Printf("Secondaries > 0.2 GeV/c/event: %.2f", ptSec->Integral(ptSec->GetXaxis()->FindBin(0.21), ptSec->GetNbinsX()) / nEvents);
+    Printf("Primaries after cuts > 0.2 GeV/c/event: %.2f +- %.2f", ptPrimCut->Integral(ptPrimCut->GetXaxis()->FindBin(0.21), ptPrimCut->GetNbinsX()) / nEvents, TMath::Sqrt(ptPrimCut->Integral(ptPrimCut->GetXaxis()->FindBin(0.21), ptPrimCut->GetNbinsX())) / nEvents);
+    Printf("Secondaries after cuts > 0.2 GeV/c/event: %.2f +- %.2f", ptSecCut->Integral(ptSecCut->GetXaxis()->FindBin(0.21), ptSecCut->GetNbinsX()) / nEvents, TMath::Sqrt(ptSecCut->Integral(ptSecCut->GetXaxis()->FindBin(0.21), ptSecCut->GetNbinsX())) / nEvents);
     Printf("%.2f %% secondaries before cuts", 100.0 * ptSec->GetEntries() / ptPrim->GetEntries());
     Printf("%.2f %% secondaries after cuts", 100.0 * ptSecCut->GetEntries() / ptPrimCut->GetEntries());
     Printf("");
 
-    ptSec->Divide(ptSec, ptPrim, 1, 1, "B");
-    ptSecCut->Divide(ptSecCut, ptPrimCut, 1, 1, "B");
-
+    ptPrim->SetLineColor(colors[i]);
+    ptPrimCut->SetLineColor(colors[i]);
     ptSec->SetLineColor(colors[i]);
     ptSecCut->SetLineColor(colors[i]);
+
+    ptPrim->SetMarkerColor(colors[i]);
+    ptPrimCut->SetMarkerColor(colors[i]);
+    ptSec->SetMarkerColor(colors[i]);
+    ptSecCut->SetMarkerColor(colors[i]);
+
+    ptPrim->SetStats(kFALSE);
+    ptPrim->SetTitle("");
+    ptPrim->GetYaxis()->SetTitle("N");
     ptSec->SetStats(kFALSE);
+    ptSec->SetTitle("");
+
+    ptPrim->SetMarkerStyle(markers[0]);
+    ptPrimCut->SetMarkerStyle(markers[1]);
+    ptSec->SetMarkerStyle(markers[2]);
+    ptSecCut->SetMarkerStyle(markers[3]);
 
+    if (i == 0) {
+      legend->AddEntry(ptPrim->Clone(), "Primaries");
+      legend->AddEntry(ptPrimCut->Clone(), "Primaries after cuts");
+      legend->AddEntry(ptSec->Clone(), "Secondaries");
+      legend->AddEntry(ptSecCut->Clone(), "Secondaries after cuts");
+    }
+
+    ptPrim->GetXaxis()->SetRangeUser(0, 2);
     ptSec->GetXaxis()->SetRangeUser(0, 2);
+    //ptPrim->GetYaxis()->SetRangeUser(1e-5, ptSec->GetMaximum() * 1.1);
+
+    c->cd();
+    ptPrim->DrawCopy((i > 0) ? "SAME" : "");
+    ptPrimCut->DrawCopy("SAME");
+    ptSec->DrawCopy("SAME");
+    ptSecCut->DrawCopy("SAME");
+
+    ptSec->Divide(ptSec, ptPrim, 1, 1, "B");
+    ptSecCut->Divide(ptSecCut, ptPrimCut, 1, 1, "B");
+
+    ptSec->SetMarkerStyle(1);
+    ptSecCut->SetMarkerStyle(1);
+
     ptSec->GetYaxis()->SetRangeUser(0, 1);
 
-    ptSec->SetTitle("");
     ptSec->GetYaxis()->SetTitle("N_{Secondaries} / N_{All}");
-    
+
+    c2->cd();
     ptSec->DrawCopy((i > 0) ? "SAME" : "");
     ptSecCut->DrawCopy("SAME");
 
     legend->AddEntry(ptSec, titles[i]);
+    legend2->AddEntry(ptSec, titles[i]);
   }
 
+  c->cd();
   legend->Draw();
+  c->SaveAs("changedmaterial_absolute.gif");
 
-  c->SaveAs("secondaries_changedmaterial.gif");
+  c2->cd();
+  legend2->Draw();
+  c2->SaveAs("changedmaterial_relative.gif");
 }
 
 void drawStats(const char* fileName = "trackCuts.root")
index 685370f..c1d4bb4 100644 (file)
@@ -33,25 +33,30 @@ void printStats(const char* fileName = "trackCuts.root")
     Printf("UNEXPECTED: %f != %f", fPrimStats->GetBinContent(5) + fPrimStats->GetBinContent(7), fStatsPrim->GetBinContent(2));
 
   Float_t notLostPrimaries = 100.0 * fPrimStats->GetBinContent(7) / (fPrimStats->GetBinContent(5) + fPrimStats->GetBinContent(7));
+
   Printf("Accepted %d out of %d primary tracks (%.2f %%)", tracksPrimAc, tracksPrim, 100.0 * tracksPrimAc / tracksPrim);
+  Printf("Per Event: %.2f out of %.2f", (Float_t) tracksPrimAc / eventsVertex, (Float_t) tracksPrim / eventsVertex);
 
   Printf("Among the non accepted ones %.2f %% are from primaries that have been found with other tracks", notLostPrimaries);
 
   Printf("Accepted %d out of %d secondary tracks (%.2f %%)", tracksSecAc, tracksSec, 100.0 * tracksSecAc / tracksSec);
+  Printf("Per Event: %.2f out of %.2f", (Float_t) tracksSecAc / eventsVertex, (Float_t) tracksSec / eventsVertex);
 
   Printf("Before cuts: %.2f %% of the tracks are secondaries", 100.0 * tracksSec / (tracksPrim + tracksSec));
 
   Printf("After cuts: %.2f %% of the tracks are secondaries", 100.0 * tracksSecAc / (tracksPrimAc + tracksSecAc));
 }
 
-void ComparePlots(const char* fileName1, const char* fileName2, const char* dir = "AliESDtrackCuts/before_cuts")
+void ComparePlots(const char* fileName1, const char* fileName2, const char* dirName1 = "AliESDtrackCuts/before_cuts", const char* dirName2 = 0)
 {
   file1 = TFile::Open(fileName1);
   file2 = TFile::Open(fileName2);
 
-  dir1 = file1->GetDirectory(dir);
-  dir2 = file2->GetDirectory(dir);
+  if (!dirName2)
+    dirName1 = dirName2;
+
+  dir1 = file1->GetDirectory(dirName1);
+  dir2 = file2->GetDirectory(dirName2);
 
   keyList = dir1->GetListOfKeys();
   TIter iter(keyList);
@@ -68,11 +73,19 @@ void ComparePlots(const char* fileName1, const char* fileName2, const char* dir
 
     TString name(key->GetName());
 
-    c = new TCanvas(key->GetName(), key->GetName(), 600, 400);
+    c = new TCanvas(Form("%s_canvas", key->GetName()), key->GetName(), 600, 400);
+
     hist1->Draw();
     hist2->SetLineColor(2);
     hist2->Draw("SAME");
 
+    /* this does not work, although given here: http://pcroot.cern.ch/root/html/THistPainter.html ;-)
+    TPaveStats *st = (TPaveStats*) hist1->FindObject("stats");
+    st->SetTextColor(2);
+    st->SetY1NDC(st->GetY1NDC() - 0.2);
+    st->SetY2NDC(st->GetY2NDC() - 0.2);
+    */
+
     Float_t min = 0.9 * TMath::Min(hist1->GetMinimum(), hist2->GetMinimum());
 
     if (name.Contains("cov") || name.Contains("dXY") || name.Contains("dZ")) {
@@ -85,7 +98,7 @@ void ComparePlots(const char* fileName1, const char* fileName2, const char* dir
     hist1->GetYaxis()->SetRangeUser(min, 1.1 * TMath::Max(hist1->GetMaximum(), hist2->GetMaximum()));
 
     c->SaveAs(Form("%s.gif", key->GetName()));
-    
+
     //break;
   }
 }