merging post-CVS changes into SVN:
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Jan 2008 08:51:07 +0000 (08:51 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 21 Jan 2008 08:51:07 +0000 (08:51 +0000)
added function to scale the corrections
draw overview of corrections
change trigger definitions to new CTP conf. file (p-p.cfg)
updated merging of corr.maps of different x-sections
moved normalizetobinwidth to alipwg0helper

24 files changed:
PWG0/AliCorrection.cxx
PWG0/AliCorrection.h
PWG0/AliCorrectionMatrix.cxx
PWG0/AliCorrectionMatrix.h
PWG0/AliCorrectionMatrix2D.cxx
PWG0/AliCorrectionMatrix3D.cxx
PWG0/AliCorrectionMatrix3D.h
PWG0/AliPWG0Helper.cxx
PWG0/AliPWG0Helper.h
PWG0/dNdEta/AliMultiplicityCorrection.cxx
PWG0/dNdEta/AliMultiplicityCorrection.h
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/CreateCuts.C
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/drawSystematics.C
PWG0/dNdEta/drawSystematicsNew.C
PWG0/dNdEta/plotsMultiplicity.C
PWG0/dNdEta/run.C
PWG0/dNdEta/runCorrection.C
PWG0/dNdEta/rundNdEtaAnalysis.C
PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx

index 46b76c3..22044fa 100644 (file)
@@ -184,7 +184,7 @@ void AliCorrection::Divide()
   fEventCorr->Divide();
   fTrackCorr->Divide();
 
-  Int_t emptyBins = fTrackCorr->CheckEmptyBins(-9.99, 9.99, -0.79, 0.79, 0.3, 9.9);
+  Int_t emptyBins = fTrackCorr->CheckEmptyBins(-9.99, 9.99, -0.79, 0.79, 0.3, 4.9);
   printf("INFO: In the central region the track correction of %s has %d empty bins\n", GetName(), emptyBins);
 }
 
@@ -284,6 +284,39 @@ void AliCorrection::DrawHistograms(const Char_t* name)
     fTrackCorr->DrawHistograms(Form("%s track", name));
 }
 
+void AliCorrection::DrawOverview(const char* canvasName)
+{
+  // draw projection of the corrections
+  //   to the 3 axis of fTrackCorr and 2 axis of fEventCorr
+
+  TString canvasNameTmp(GetName());
+  if (canvasName)
+    canvasNameTmp = canvasName;
+
+  TCanvas* canvas = new TCanvas(canvasNameTmp, canvasNameTmp, 1200, 800);
+  canvas->Divide(3, 2);
+
+  if (fTrackCorr) {
+    canvas->cd(1);
+    fTrackCorr->Get1DCorrectionHistogram("x", 0.3, 5, -1, 1)->DrawCopy()->GetYaxis()->SetRangeUser(0, 10);
+
+    canvas->cd(2);
+    fTrackCorr->Get1DCorrectionHistogram("y", 0.3, 5, 0, 0)->DrawCopy()->GetYaxis()->SetRangeUser(0, 10);
+
+    canvas->cd(3);
+    fTrackCorr->Get1DCorrectionHistogram("z", 0, -1, -1, 1)->DrawCopy()->GetYaxis()->SetRangeUser(0, 10);
+  }
+
+  if (fEventCorr)
+  {
+    canvas->cd(4);
+    fEventCorr->Get1DCorrectionHistogram("x")->DrawCopy();
+
+    canvas->cd(5);
+    fEventCorr->Get1DCorrectionHistogram("y")->DrawCopy()->GetXaxis()->SetRangeUser(0, 30);
+  }
+}
+
 //____________________________________________________________________
 void AliCorrection::SetCorrectionToUnity()
 {
@@ -313,3 +346,35 @@ void AliCorrection::Multiply()
   if (fTrackCorr)
     fTrackCorr->Multiply();
 }
+
+//____________________________________________________________________
+void AliCorrection::Scale(Double_t factor)
+{
+  // scales the two contained corrections
+
+  fEventCorr->Scale(factor);
+  fTrackCorr->Scale(factor);
+}
+
+//____________________________________________________________________
+void AliCorrection::PrintInfo(Float_t ptCut)
+{
+  // prints some stats
+
+  TH3* measured = GetTrackCorrection()->GetMeasuredHistogram();
+  TH3* generated = GetTrackCorrection()->GetGeneratedHistogram();
+
+  TH2* measuredEvents = GetEventCorrection()->GetMeasuredHistogram();
+  TH2* generatedEvents = GetEventCorrection()->GetGeneratedHistogram();
+
+  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());
+
+  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());
+
+  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);
+}
index 788c31a..84b8d78 100644 (file)
@@ -40,16 +40,19 @@ public:
   void Divide();
   void Multiply();
   void SetCorrectionToUnity();
+  void Scale(Double_t factor);
 
   void Add(AliCorrection* aCorrectionToAdd, Float_t c=1);
 
   virtual Bool_t LoadHistograms(const Char_t* dir = 0);
   virtual void SaveHistograms();
   virtual void DrawHistograms(const Char_t* name = 0);
+  virtual void DrawOverview(const char* canvasName = 0);
 
   virtual void ReduceInformation();
 
   virtual void Reset(Option_t* option = "");
+  void PrintInfo(Float_t ptCut);
 
 protected:
   AliCorrectionMatrix2D* fEventCorr; // correction on event level
index 28061c3..990a1d4 100644 (file)
@@ -364,3 +364,14 @@ void AliCorrectionMatrix::SetCorrectionToUnity()
       for (Int_t z=1; z<=fhCorr->GetNbinsZ(); ++z)
         fhCorr->SetBinContent(x, y, z, 1);
 }
+
+//____________________________________________________________________
+void AliCorrectionMatrix::Scale(Double_t factor)
+{
+  // scales the generated and measured histogram with the given factor
+
+  Printf("Scaling histograms with %f", factor);
+
+  fhMeas->Scale(factor);
+  fhGene->Scale(factor);
+}
index 359fa82..a7c2c54 100644 (file)
@@ -42,6 +42,7 @@ public:
   void Divide();
   void Multiply();
   void SetCorrectionToUnity();
+  void Scale(Double_t factor);
 
   void Add(AliCorrectionMatrix* aMatrixToAdd, Float_t c=1);
 
index 39659ce..e208f43 100644 (file)
@@ -132,25 +132,25 @@ TH1F* AliCorrectionMatrix2D::Get1DCorrectionHistogram(Char_t* opt, Float_t min,
       gene1D = ((TH2F*)fhGene)->ProjectionX();
     }
     else {
-      AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including y-bins %d to %d \n", binMin, binMax));
+      Printf("Getting 1D map. Including y-bins %d to %d", binMin, binMax);
 
-      meas1D = ((TH2F*)fhMeas)->ProjectionX("pm",binMin,binMax);
-      gene1D = ((TH2F*)fhGene)->ProjectionX("pg",binMin,binMax);
+      meas1D = ((TH2F*)fhMeas)->ProjectionX(Form("%s_pm", GetName()),binMin,binMax);
+      gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_pg", GetName()),binMin,binMax);
     }
   }
   if (strcmp(opt,"y")==0) {
     Int_t binMin = fhMeas->GetXaxis()->FindBin(min);
     Int_t binMax = fhMeas->GetXaxis()->FindBin(max);
-    
+
     if (min>=max) {
       meas1D = ((TH2F*)fhMeas)->ProjectionY();
       gene1D = ((TH2F*)fhGene)->ProjectionY();
     }
     else {
       AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax));
-      
-      meas1D = ((TH2F*)fhMeas)->ProjectionY("pm", binMin, binMax);
-      gene1D = ((TH2F*)fhGene)->ProjectionY("pg", binMin, binMax);
+
+      meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_pm", GetName()), binMin, binMax);
+      gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_pg", GetName()), binMin, binMax);
     }
   }
   gene1D->Sumw2();
index d50d389..b953f3b 100644 (file)
@@ -189,33 +189,46 @@ AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Char_t* opt, Float
 
   TString option = opt;
 
+  // unzoom
+  fhMeas->GetXaxis()->UnZoom();
+  fhMeas->GetYaxis()->UnZoom();
+  fhMeas->GetZaxis()->UnZoom();
+
+  fhGene->GetXaxis()->UnZoom();
+  fhGene->GetYaxis()->UnZoom();
+  fhGene->GetZaxis()->UnZoom();
+
   if (aMin<aMax) {
     if (option.Contains("xy") || option.Contains("yx")) {
       Int_t bMin = fhMeas->GetZaxis()->FindBin(aMin);
       Int_t bMax = fhMeas->GetZaxis()->FindBin(aMax);
-      fhMeas->GetZaxis()->SetRange(bMin, bMax);      
+      fhGene->GetZaxis()->SetRange(bMin, bMax);
+      fhMeas->GetZaxis()->SetRange(bMin, bMax);
     }
     else if (option.Contains("xz") || option.Contains("zx")) {
       Int_t bMin = fhMeas->GetYaxis()->FindBin(aMin);
       Int_t bMax = fhMeas->GetYaxis()->FindBin(aMax);
-      fhMeas->GetYaxis()->SetRange(bMin, bMax);      
+      fhGene->GetYaxis()->SetRange(bMin, bMax);
+      fhMeas->GetYaxis()->SetRange(bMin, bMax);
     }
     else if (option.Contains("yz") || option.Contains("zy")) {
       Int_t bMin = fhMeas->GetXaxis()->FindBin(aMin);
       Int_t bMax = fhMeas->GetXaxis()->FindBin(aMax);
-      fhMeas->GetXaxis()->SetRange(bMin, bMax);      
+      fhGene->GetXaxis()->SetRange(bMin, bMax);
+      fhMeas->GetXaxis()->SetRange(bMin, bMax);
     }
     else {
       AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
       return 0;
     }
   }
+
   AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);
 
-  TH2F* meas = (TH2F*) ((TH3F*)fhMeas)->Project3D(opt);
-  TH2F* gene = (TH2F*) ((TH3F*)fhGene)->Project3D(opt);
+  TH2F* meas = (TH2F*) ((TH3F*)fhMeas)->Project3D(option)->Clone(Form("%s_meas", corr2D->GetName()));
+  TH2F* gene = (TH2F*) ((TH3F*)fhGene)->Project3D(option)->Clone(Form("%s_gene", corr2D->GetName()));
 
-  TH2F* corr = (TH2F*)gene->Clone("corr");
+  TH2F* corr = (TH2F*)gene->Clone(Form("%s_corr", corr2D->GetName()));
   corr->Reset();
 
   corr2D->SetGeneratedHistogram(gene);
@@ -225,9 +238,9 @@ AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Char_t* opt, Float
   corr2D->Divide();
 
   // unzoom
-  fhMeas->GetXaxis()->UnZoom();  
-  fhMeas->GetYaxis()->UnZoom();  
-  fhMeas->GetZaxis()->UnZoom();  
+  fhMeas->GetXaxis()->UnZoom();
+  fhMeas->GetYaxis()->UnZoom();
+  fhMeas->GetZaxis()->UnZoom();
 
   fhGene->GetXaxis()->UnZoom();
   fhGene->GetYaxis()->UnZoom();
@@ -243,16 +256,16 @@ TH1F* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Char_t* opt, Float_t aMin1
   AliDebug(AliLog::kWarning, Form("WARNING: test"));
   
   AliCorrectionMatrix2D* corr2D;
-  if (strcmp(opt,"x")==0) {  
-    corr2D = Get2DCorrection("xy",aMin1,aMax1);
+  if (strcmp(opt,"x")==0) {
+    corr2D = Get2DCorrection("yx",aMin1,aMax1);
     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
   }
   if (strcmp(opt,"y")==0) {
     corr2D = Get2DCorrection("xy",aMin1,aMax1);
-    return corr2D->Get1DCorrectionHistogram("y",aMin2,aMax2);
+    return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
   }  
   if (strcmp(opt,"z")==0) {
-    corr2D = Get2DCorrection("yz",aMin1,aMax1);    
+    corr2D = Get2DCorrection("yz",aMin1,aMax1);
     return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
   }  
   AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s (should be x,y or z)", opt));  
@@ -260,8 +273,6 @@ TH1F* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Char_t* opt, Float_t aMin1
   return 0;
 }
 
-
-
 //____________________________________________________________________
 void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az)
 {
index f1a2826..a9d090c 100644 (file)
@@ -48,7 +48,7 @@ public:
 
   AliCorrectionMatrix2D* Get2DCorrection(Char_t* opt, Float_t aMin, Float_t aMax);
   TH2F* Get2DCorrectionHistogram(Char_t* opt, Float_t aMin, Float_t aMax) {return Get2DCorrection(opt,aMin,aMax)->GetCorrectionHistogram();}
-  TH1F* Get1DCorrectionHistogram(Char_t* opt, Float_t aMins1, Float_t aMax1, Float_t aMins2, Float_t aMax2);
+  TH1F* Get1DCorrectionHistogram(Char_t* opt, Float_t aMins1=0, Float_t aMax1=0, Float_t aMins2=0, Float_t aMax2=0);
 
   void FillMeas(Float_t ax, Float_t ay, Float_t az);
   void FillGene(Float_t ax, Float_t ay, Float_t az);
index 51432ea..7377dbb 100644 (file)
@@ -5,6 +5,7 @@
 #include <TParticle.h>
 #include <TParticlePDG.h>
 #include <TH1.h>
+#include <TH2.h>
 #include <TH3.h>
 #include <TList.h>
 #include <TTree.h>
@@ -41,23 +42,24 @@ Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
   // check if the event was triggered
   //
   // this function needs the branch fTriggerMask
-  //
-  // MB should be
-  // ITS_SPD_GFO_L0  : 32
-  // VZERO_OR_LEFT   : 1
-  // VZERO_OR_RIGHT  : 2
+  
+
+  // definitions from p-p.cfg
+  ULong64_t spdFO = (1 << 14);
+  ULong64_t v0left = (1 << 11);
+  ULong64_t v0right = (1 << 12);
 
   switch (trigger)
   {
     case kMB1:
     {
-      if (triggerMask&32 || ((triggerMask&1) || (triggerMask&2)))
+      if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
         return kTRUE;
       break;
     }
     case kMB2:
     {
-      if (triggerMask&32 && ((triggerMask&1) || (triggerMask&2)))
+      if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
         return kTRUE;
       break;
     }
@@ -89,6 +91,10 @@ Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD)
   if (strcmp(vtxESD->GetName(), "default")==0)
     return kFALSE;
 
+  // check Ncontributors
+  if (vtxESD->GetNContributors() <= 0)
+    return kFALSE;
+
   Double_t vtx_res[3];
   vtx_res[0] = vtxESD->GetXRes();
   vtx_res[1] = vtxESD->GetYRes();
@@ -97,8 +103,6 @@ Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD)
   if (vtx_res[2]==0 || vtx_res[2]>0.1)
     return kFALSE;
 
-  // check Ncontributors, if <0 it means error *gna*
-
   return kTRUE;
 }
 
@@ -111,6 +115,8 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari
   //
   // This function or a equivalent should be available in some common place of AliRoot
   //
+  // WARNING: Call this function only for particles that are among the particles from the event generator!
+  // --> stack->Particle(id) with id < stack->GetNprimary()
 
   // if the particle has a daughter primary, we do not want to count it
   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
@@ -396,3 +402,33 @@ void AliPWG0Helper::SetBranchStatusRecursive(TTree* tree, char *bname, Bool_t st
     }
   }
 }
+
+//____________________________________________________________________
+void AliPWG0Helper::NormalizeToBinWidth(TH1* hist)
+{
+  //
+  // normalizes a 1-d histogram to its bin width
+  //
+
+  for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
+  {
+    hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
+    hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
+  }
+}
+
+//____________________________________________________________________
+void AliPWG0Helper::NormalizeToBinWidth(TH2* hist)
+{
+  //
+  // normalizes a 2-d histogram to its bin width (x width * y width)
+  //
+
+  for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
+    for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
+    {
+      Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
+      hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
+      hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
+    }
+}
index 04fb01d..5e9fdc2 100644 (file)
@@ -10,6 +10,8 @@
 class AliESD;
 class AliESDVertex;
 class TParticle;
+class TH1;
+class TH2;
 class TH3;
 class AliHeader;
 class AliStack;
@@ -36,6 +38,9 @@ class AliPWG0Helper : public TObject
 
     static void SetBranchStatusRecursive(TTree* tree, char *bname, Bool_t status, Bool_t debug = kFALSE);
 
+    static void NormalizeToBinWidth(TH1* hist);
+    static void NormalizeToBinWidth(TH2* hist);
+
   protected:
     ClassDef(AliPWG0Helper, 0)
 
index c96b548..94abd16 100644 (file)
@@ -1224,36 +1224,6 @@ void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSp
   }
 }
 
-//____________________________________________________________________
-void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
-{
-  //
-  // normalizes a 1-d histogram to its bin width
-  //
-
-  for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
-  {
-    hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
-    hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
-  }
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
-{
-  //
-  // normalizes a 2-d histogram to its bin width (x width * y width)
-  //
-
-  for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
-    for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
-    {
-      Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
-      hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
-      hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
-    }
-}
-
 //____________________________________________________________________
 void AliMultiplicityCorrection::DrawHistograms()
 {
index 99ff37b..5a4e4b2 100644 (file)
@@ -85,9 +85,6 @@ class AliMultiplicityCorrection : public TNamed {
     void SetGenMeasFromFunc(TF1* inputMC, Int_t id);
     TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap);
 
-    static void NormalizeToBinWidth(TH1* hist);
-    static void NormalizeToBinWidth(TH2* hist);
-
     void GetComparisonResults(Float_t* mc = 0, Int_t* mcLimit = 0, Float_t* residuals = 0, Float_t* ratioAverage = 0) const;
 
     TH1* GetEfficiency(Int_t inputRange, EventType eventType);
index 48def25..39f8a32 100644 (file)
@@ -239,7 +239,7 @@ void AlidNdEtaCorrection::SaveHistograms()
 void AlidNdEtaCorrection::DrawHistograms()
 {
   //
-  // call the draw histogram method of the correction
+  // call the draw histogram method of the corrections
   //
 
   fTrack2ParticleCorrection     ->DrawHistograms();
@@ -249,6 +249,20 @@ void AlidNdEtaCorrection::DrawHistograms()
   fTriggerBiasCorrectionMBToND  ->DrawHistograms();
 }
 
+//____________________________________________________________________
+void AlidNdEtaCorrection::DrawOverview(const char* canvasName)
+{
+  //
+  // call the DrawOverview histogram method of the corrections
+  //
+
+  fTrack2ParticleCorrection     ->DrawOverview(canvasName);
+  fVertexRecoCorrection         ->DrawOverview(canvasName);
+  fTriggerBiasCorrectionMBToINEL->DrawOverview(canvasName);
+  fTriggerBiasCorrectionMBToNSD ->DrawOverview(canvasName);
+  fTriggerBiasCorrectionMBToND  ->DrawOverview(canvasName);
+}
+
 //____________________________________________________________________
 void AlidNdEtaCorrection::FillMCParticle(Float_t vtx, Float_t eta, Float_t pt, Bool_t trigger, Bool_t vertex, Int_t processType)
 {
index 7d54dc4..59502ff 100644 (file)
@@ -61,6 +61,7 @@ public:
   void    SaveHistograms();
   Bool_t  LoadHistograms(const Char_t* dir = 0);
   void    DrawHistograms();
+  void    DrawOverview(const char* canvasName = 0);
 
   Float_t GetMeasuredFraction(CorrectionType correctionType, Float_t ptCutOff, Float_t eta = -100, Bool_t debug = kFALSE);
 
index 9eb8f98..7a77952 100644 (file)
@@ -270,38 +270,37 @@ void AlidNdEtaTask::Exec(Option_t*)
   // Processing of ESD information (always)
   if (eventTriggered)
   {
-    // control hists
+    // control hist
     fMult->Fill(inputCount);
     if (eventVertex)
+    {
+      // control hist
       fMultVtx->Fill(inputCount);
 
-    // count all events
-    if (!eventVertex)
-      vtx[2] = 0;
-
-    for (Int_t i=0; i<inputCount; ++i)
-    {
-      Float_t eta = etaArr[i];
-      Float_t pt  = ptArr[i];
+      for (Int_t i=0; i<inputCount; ++i)
+      {
+        Float_t eta = etaArr[i];
+        Float_t pt  = ptArr[i];
 
-      fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
+        fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
 
-      if (TMath::Abs(vtx[2]) < 20)
-      {
-        fPartEta[0]->Fill(eta);
+        if (TMath::Abs(vtx[2]) < 20)
+        {
+          fPartEta[0]->Fill(eta);
 
-        if (vtx[2] < 0)
-          fPartEta[1]->Fill(eta);
-        else
-          fPartEta[2]->Fill(eta);
+          if (vtx[2] < 0)
+            fPartEta[1]->Fill(eta);
+          else
+            fPartEta[2]->Fill(eta);
+        }
       }
-    }
 
-    // for event count per vertex
-    fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
+      // for event count per vertex
+      fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
 
-    // control hists
-    fEvents->Fill(vtx[2]);
+      // control hists
+      fEvents->Fill(vtx[2]);
+    }
   }
 
   if (fReadMC)   // Processing of MC information (optional)
index 15baf3f..6a150f2 100644 (file)
@@ -11,8 +11,9 @@ AliESDtrackCuts* CreateTrackCuts(Bool_t hists = kTRUE)
 
   esdTrackCuts->SetMinNClustersTPC(50);
   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
-  esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,1e99);
-  Printf("WARNING: no cut on 1/pt");
+  esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
+  //esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,1e99);
+  //Printf("WARNING: no cut on 1/pt");
   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 
   esdTrackCuts->SetMinNsigmaToVertex(3);
index 08bc84e..f937ab5 100644 (file)
@@ -185,11 +185,11 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
       trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
 
-    /*if (correctionType >= AlidNdEtaCorrection::kVertexReco)
+    if (correctionType >= AlidNdEtaCorrection::kVertexReco)
     {
       trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
       eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
-    }*/
+    }
 
     switch (correctionType)
     {
@@ -361,12 +361,12 @@ 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", bin, dndeta);
+        //Printf("Bin %d has dN/deta = %f", bin, dndeta);
       }
     }
   }
 
-  new TCanvas; fdNdEta[0]->DrawCopy();
+  //new TCanvas; fdNdEta[0]->DrawCopy();
 }
 
 //____________________________________________________________________
index 3723788..933308a 100644 (file)
@@ -20,6 +20,12 @@ Int_t gMax = 5;
 
 extern TSystem* gSystem;
 
+void loadlibs()
+{
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libPWG0base");
+}
+
 void SetRanges(TAxis* axis)
 {
   if (strcmp(axis->GetTitle(), "#eta") == 0)
@@ -94,6 +100,23 @@ void InitPadCOLZ()
   gPad->SetGridy();
 }
 
+// --- end of helpers --- begin functions ---
+
+void DrawOverview(const char* fileName, const char* dirName)
+{
+  loadlibs();
+  if (!TFile::Open(fileName))
+    return;
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
+  if (!dNdEtaCorrection->LoadHistograms())
+    return;
+
+  dNdEtaCorrection->Finish();
+
+  dNdEtaCorrection->DrawOverview();
+}
+
 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");
@@ -202,6 +225,28 @@ void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "corr
   diff->Draw();
 }
 
+Double_t PrintIntegratedDeviation(TH1* histMC, TH1* histESD, const char* desc = "")
+{
+  Double_t avgMC = 0;
+  Double_t avgESD = 0;
+  for (Int_t bin = histMC->FindBin(-0.7999); bin <= histMC->FindBin(0.7999); bin++)
+  {
+    avgMC += histMC->GetBinContent(bin);
+    avgESD += histESD->GetBinContent(bin);
+  }
+  Int_t nBins = histMC->FindBin(0.7999) - histMC->FindBin(-0.7999) + 1;
+
+  avgMC /= nBins;
+  avgESD /= nBins;
+
+  // deviation when integrate in |eta| < 0.8 between mc and esd
+  Double_t diffFullRange = (avgMC - avgESD) / avgMC;
+
+  Printf("%s: Integrated deviation in |eta| < 0.8 is %.2f %%", desc, diffFullRange * 100);
+
+  return diffFullRange;
+}
+
 void dNdEta(Bool_t onlyESD = kFALSE)
 {
   TFile* file = TFile::Open("analysis_esd.root");
@@ -279,7 +324,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   if (onlyESD)
     return;
 
-  gSystem->Load("libPWG0base");
+  loadlibs();
 
   TFile* file2 = TFile::Open("analysis_mc.root");
   TH1* histMC = (TH1*) file2->Get("dndeta/dNdEta_corrected")->Clone("cloned");
@@ -363,6 +408,20 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   ratio->SetLineColor(1);
   ratioNoPt->SetLineColor(2);
 
+  Double_t average = 0;       // average deviation from 1 in ratio (depends on the number of bins if statistical)
+  for (Int_t bin = ratio->FindBin(-0.7999); bin <= ratio->FindBin(0.7999); bin++)
+    average += TMath::Abs(ratio->GetBinContent(bin) - 1);
+  Int_t nBins = ratio->FindBin(0.7999) - ratio->FindBin(-0.7999) + 1;
+  average /= nBins;
+  Printf("Average deviation in |eta| < 0.8 is %.2f %%", average * 100);
+
+  PrintIntegratedDeviation(histMC, histESD, "all events");
+  PrintIntegratedDeviation(histMCTr, histESDMB, "triggered");
+  PrintIntegratedDeviation(histMCTrVtx, histESDMBVtx, "trigger, vertex");
+  PrintIntegratedDeviation(histMCPtCut, histESDNoPt, "all events (no pt corr)");
+  PrintIntegratedDeviation(histMCTrPtCut, histESDMBNoPt, "triggered (no pt corr)");
+  PrintIntegratedDeviation(histMCTrVtxPtCut, histESDMBVtxNoPt, "trigger, vertex (no pt corr)");
+
   TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 600);
   canvas3->Range(0, 0, 1, 1);
   //canvas3->Divide(1, 2, 0, 0);
index 30ecad6..596e040 100644 (file)
@@ -26,6 +26,16 @@ void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", c
 Int_t markers[] = {20,20,21,22,23,28,29};
 Int_t colors[]  = {1,2,3,4,6,8,102};
 
+void loadlibs()
+{
+  gSystem->Load("libTree");
+  gSystem->Load("libVMC");
+
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libPWG0base");
+}
+
 void SetRanges(TAxis* axis)
 {
   if (strcmp(axis->GetTitle(), "#eta") == 0)
@@ -748,21 +758,39 @@ TH1F* Sigma2VertexGaussian()
   return ratio2;
 }
 
-TH1F* Sigma2VertexSimulation()
+TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
 {
-  TFile* file = TFile::Open("systematics.root");
+  TFile* file = TFile::Open(fileName);
 
-  TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex"));
-  if (!sigmavertex)
+  TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
+  TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
+  if (!sigmavertex || !sigmavertexPrim)
   {
-    printf("Could not read histogram\n");
+    printf("Could not read histogram(s)\n");
     return;
   }
 
+  // calculate ratio
   TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 3 #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(3)) / 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 < 4.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);
@@ -775,13 +803,22 @@ TH1F* Sigma2VertexSimulation()
   ratio->SetMarkerStyle(21);
   ratio->DrawCopy("P");
 
-  return ratio;
+  contamination->DrawCopy("SAME");
+
+  TH1F** returnContainer = new TH1F*[2];
+  returnContainer[0] = ratio;
+  returnContainer[1] = contamination;
+
+  return returnContainer;
 }
 
-void Sigma2VertexCompare()
+void Sigma2VertexCompare(const char* fileName = "systematics.root")
 {
   TH1F* ratio1 = Sigma2VertexGaussian();
-  TH1F* ratio2 = Sigma2VertexSimulation();
+
+  TH1F** hists = Sigma2VertexSimulation(fileName);
+  TH1F* ratio2 = hists[0];
+  TH1F* contamination = hists[1];
 
   ratio1->SetStats(kFALSE);
   ratio2->SetStats(kFALSE);
@@ -792,14 +829,15 @@ void Sigma2VertexCompare()
   ratio1->SetLineWidth(2);
   ratio2->SetLineWidth(2);
 
-  TLegend* legend = new TLegend(0.7, 0.8, 0.95, 0.95);
+  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, 4);
+  ratio2->GetXaxis()->SetRangeUser(2, 5);
 
   TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
   InitPad();
@@ -813,6 +851,8 @@ void Sigma2VertexCompare()
   ratio2->Draw("PL");
   ratio1->Draw("SAMEPL");
 
+  contamination->Draw("SAME");
+
   legend->Draw();
 
   canvas->SaveAs("Sigma2VertexCompare.eps");
@@ -923,135 +963,147 @@ void DrawdNdEtaDifferences()
   canvas2->SaveAs("particlecomposition_result.eps");
 }
 
-mergeCorrectionsWithDifferentCrosssections(Char_t* standardCorrectionFileName="correction_map.root",
-                                            Char_t* systematicCorrectionFileName="systematics.root",
-                                            Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
+mergeCorrectionsWithDifferentCrosssections(Char_t* correctionFileName="correction_map.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
   // and SD events.
-  // 
+  //
+  // the dn/deta spectrum is corrected and the ratios
+  // (standard to changed x-section) of the different dN/deta
+  // distributions are saved to a file.
+  //
 
-  gSystem->Load("libPWG0base");
+  loadlibs();
 
   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
-  const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
-  Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5};
-  Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5};
+  const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25"};
+  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};
 
   // cross section from Pythia
   Float_t sigmaND = 55.2;
   Float_t sigmaDD = 9.78;
   Float_t sigmaSD = 14.30;
 
-  AlidNdEtaCorrection* corrections[21];
+  // standard correction
+  TFile::Open(correctionFileName);
+  AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
+  correctionStandard->LoadHistograms();
+
+  // dont take vertexreco from this one
+  correctionStandard->GetVertexRecoCorrection()->Reset();
+  // dont take triggerbias from this one
+  correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
+  correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
+  correctionStandard->GetTriggerBiasCorrectionND()->Reset();
+
+  AlidNdEtaCorrection* corrections[100];
+  TH1F* hRatios[100];
+
+  Int_t counter = 0;
   for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
-    AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
-    correctionStandard->LoadHistograms(standardCorrectionFileName);
 
-    // dont take vertexreco from this one
-    correctionStandard->GetVertexRecoCorrection()->Reset();
-    // dont take triggerbias from this one
-    correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
+    for (Int_t i=0; i<17; i++) {
+      TFile::Open(correctionFileName);
 
-    for (Int_t i=0; i<7; i++) {
       TString name;
       name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
       AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
 
-      name.Form("vertexRecoND");
-      AlidNdEtaCorrection* dNdEtaCorrectionVtxND = new AlidNdEtaCorrection(name,name);
-      dNdEtaCorrectionVtxND->LoadHistograms(systematicCorrectionFileName, name);
-      name.Form("vertexRecoDD");
-      AlidNdEtaCorrection* dNdEtaCorrectionVtxDD = new AlidNdEtaCorrection(name,name);
-      dNdEtaCorrectionVtxDD->LoadHistograms(systematicCorrectionFileName, name);
-      name.Form("vertexRecoSD");
-      AlidNdEtaCorrection* dNdEtaCorrectionVtxSD = new AlidNdEtaCorrection(name,name);
-      dNdEtaCorrectionVtxSD->LoadHistograms(systematicCorrectionFileName, name);
-
-      name.Form("triggerBiasND");
-      AlidNdEtaCorrection* dNdEtaCorrectionTriggerND = new AlidNdEtaCorrection(name,name);
-      dNdEtaCorrectionTriggerND->LoadHistograms(systematicCorrectionFileName, name);
-      name.Form("triggerBiasDD");
-      AlidNdEtaCorrection* dNdEtaCorrectionTriggerDD = new AlidNdEtaCorrection(name,name);
-      dNdEtaCorrectionTriggerDD->LoadHistograms(systematicCorrectionFileName, name);
-      name.Form("triggerBiasSD");
-      AlidNdEtaCorrection* dNdEtaCorrectionTriggerSD = new AlidNdEtaCorrection(name,name);
-      dNdEtaCorrectionTriggerSD->LoadHistograms(systematicCorrectionFileName, name);
+      name.Form("dndeta_correction_ND");
+      AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
+      dNdEtaCorrectionND->LoadHistograms();
+      name.Form("dndeta_correction_DD");
+      AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
+      dNdEtaCorrectionDD->LoadHistograms();
+      name.Form("dndeta_correction_SD");
+      AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
+      dNdEtaCorrectionSD->LoadHistograms();
 
       // calculating relative
       Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
       Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
       Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
 
-      printf(Form("%s : ND=%.1f\%, DD=%.1f\%, SD=%.1f\% \n",changes[i],nd,dd,sd));
+      printf(Form("%s : ND=%.2f\%, DD=%.2f\%, SD=%.2f\% \n",changes[i],nd,dd,sd));
       current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",nd,dd,sd));
       current->SetTitle(name);
 
       // scale
       if (j == 0 || j == 2)
       {
-        dNdEtaCorrectionVtxDD->GetVertexRecoCorrection()->GetMeasuredHistogram()->Scale(scalesDD[i]);
-        dNdEtaCorrectionVtxSD->GetVertexRecoCorrection()->GetMeasuredHistogram()->Scale(scalesSD[i]);
-        dNdEtaCorrectionVtxDD->GetVertexRecoCorrection()->GetGeneratedHistogram()->Scale(scalesDD[i]);
-        dNdEtaCorrectionVtxSD->GetVertexRecoCorrection()->GetGeneratedHistogram()->Scale(scalesSD[i]);
+        dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scalesSD[i]);
       }
       if (j == 1 || j == 2)
       {
-        dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionINEL()->GetMeasuredHistogram()->Scale(scalesDD[i]);
-        dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionINEL()->GetMeasuredHistogram()->Scale(scalesSD[i]);
-        dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionINEL()->GetGeneratedHistogram()->Scale(scalesDD[i]);
-        dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionINEL()->GetGeneratedHistogram()->Scale(scalesSD[i]);
+        dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scalesSD[i]);
+
+        dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scalesSD[i]);
+
+        dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scalesSD[i]);
       }
 
-      //clear track, trigger in Vtx correction
-      dNdEtaCorrectionVtxND->GetTrack2ParticleCorrection()->Reset();
-      dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionNSD()->Reset();
-      dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionND()->Reset();
-      dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionINEL()->Reset();
-      dNdEtaCorrectionVtxSD->GetTrack2ParticleCorrection()->Reset();
-      dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionNSD()->Reset();
-      dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionND()->Reset();
-      dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionINEL()->Reset();
-      dNdEtaCorrectionVtxDD->GetTrack2ParticleCorrection()->Reset();
-      dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionNSD()->Reset();
-      dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionND()->Reset();
-      dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionINEL()->Reset();
-
-      //clear track, vertexreco in trigger correction
-      dNdEtaCorrectionTriggerND->GetTrack2ParticleCorrection()->Reset();
-      dNdEtaCorrectionTriggerND->GetTriggerBiasCorrectionNSD()->Reset();
-      dNdEtaCorrectionTriggerND->GetTriggerBiasCorrectionND()->Reset();
-      dNdEtaCorrectionTriggerND->GetVertexRecoCorrection()->Reset();
-      dNdEtaCorrectionTriggerSD->GetTrack2ParticleCorrection()->Reset();
-      dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionNSD()->Reset();
-      dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionND()->Reset();
-      dNdEtaCorrectionTriggerSD->GetVertexRecoCorrection()->Reset();
-      dNdEtaCorrectionTriggerDD->GetTrack2ParticleCorrection()->Reset();
-      dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionNSD()->Reset();
-      dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionND()->Reset();
-      dNdEtaCorrectionTriggerDD->GetVertexRecoCorrection()->Reset();
+      //clear track in correction
+      dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
+      dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
+      dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
 
       TList collection;
       collection.Add(correctionStandard);
-      collection.Add(dNdEtaCorrectionVtxND);
-      collection.Add(dNdEtaCorrectionVtxDD);
-      collection.Add(dNdEtaCorrectionVtxSD);
-      collection.Add(dNdEtaCorrectionTriggerND);
-      collection.Add(dNdEtaCorrectionTriggerDD);
-      collection.Add(dNdEtaCorrectionTriggerSD);
+      collection.Add(dNdEtaCorrectionND);
+      collection.Add(dNdEtaCorrectionDD);
+      collection.Add(dNdEtaCorrectionSD);
 
       current->Merge(&collection);
       current->Finish();
 
-      corrections[i+j*7] = current;
+      corrections[counter] = current;
+
+      // now correct dNdeta distribution with modified correction map
+      TFile::Open(analysisFileName);
+
+      dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
+      fdNdEtaAnalysis->LoadHistograms();
+
+      fdNdEtaAnalysis->Finish(current, 0.3, AlidNdEtaCorrection::kINEL);
+      //fdNdEtaAnalysis->Finish(current, 0.3, AlidNdEtaCorrection::kNSD);
+
+      name = "ratio";
+      if (j==0) name.Append("_vetexReco_");
+      if (j==1) name.Append("_triggerBias_");
+      if (j==2) name.Append("_vertexReco_triggerBias_");
+      name.Append(changes[i]);
+
+      hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(name);
+
+      name.Append(Form(" (DD #times %0.2f, SD #times %0.2f)",scalesDD[i],scalesSD[i]));
+      hRatios[counter]->SetTitle(name.Data());
+      hRatios[counter]->SetYTitle("ratio (Pythia x-sections)/(changed x-sections)");
+
+      if (counter > 0)
+        hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
+
+      delete fdNdEtaAnalysis;
+
+      counter++;
     }
   }
 
   TFile* fout = new TFile(outputFileName,"RECREATE");
 
-  for (Int_t i=0; i<21; i++)
+  // to make everything consistent
+  hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
+
+  for (Int_t i=0; i<51; i++)
+  {
     corrections[i]->SaveHistograms();
+    hRatios[i]->Write();
+  }
 
   fout->Write();
   fout->Close();
index 11aba80..de420b0 100644 (file)
@@ -2,6 +2,16 @@
 Int_t markers[] = {20,21,22,23,28,29};
 Int_t colors[]  = {1,2,3,4,6,8,102};
 
+void loadlibs()
+{
+  gSystem->Load("libTree");
+  gSystem->Load("libVMC");
+
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libPWG0base");
+}
+
 void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
 {
   //gROOT->ProcessLine(".L drawPlots.C");
@@ -42,35 +52,39 @@ void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
 }
 
 
-void DrawEffectOfChangeInCrossSection() {
+void DrawEffectOfChangeInCrossSection(const char* fileName) {
   
   //get the data
-  TFile* fin = TFile::Open("systematics_xsections.root");
+  TFile* fin = TFile::Open(fileName);
 
-  const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
-  Int_t colors[] = {1,2,103,102,4,6,1};
+  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" };
+  Int_t colors[] = {1,1,kRed,kBlue,kGreen,kPink,1,kRed,kBlue};
 
-  TH1F* hRatios[7];
-  for(Int_t i=0; i<7; i++) {
+  TH1F* hRatios[9];
+  for(Int_t i=0; i<9; i++) {
     hRatios[i] = (TH1F*)fin->Get(Form("ratio_vertexReco_triggerBias_%s",changes[i]));
     hRatios[i]->SetLineWidth(2);
-    hRatios[i]->SetLineColor(colors[i]);                                
+    hRatios[i]->SetLineColor(colors[i]);
     hRatios[i]->SetMarkerStyle(22);
     hRatios[i]->SetMarkerSize(0.8);
+
+    Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
+    Printf("%s: %.2f %%" , hRatios[i]->GetTitle(), (average - 1) * 100);
   }
 
   TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
   p->SetRightMargin(0.2);
   p->SetLeftMargin(0.13);
 
-  TH2F* null = new TH2F("","",100,-1.05,1.05,100,0.93,1.06);
+  TH2F* null = new TH2F("","",100,-1.05,1.05,100,0.9,1.1);
   null->GetXaxis()->SetTitle("#eta");
   null->GetYaxis()->SetTitle("Ratio pythia/modified cross-sections");
   null->Draw();
 
-  TLatex* text[7];
+  TLatex* text[9];
 
-  for(Int_t i=1; i<7; i++) {
+  for(Int_t i=1; i<9; i++) {
     hRatios[i]->Draw("same");
     
     TString str(hRatios[i]->GetTitle());
@@ -169,16 +183,17 @@ TPad* DrawCanvasAndPad(const Char_t* name, Int_t sizeX=600, Int_t sizeY=500) {
   return p1;
 }
 
-void MisalignmentShowRawTrackPlots()
+void MisalignmentShowRawTrackPlots(const char* dirName = "fdNdEtaAnalysisESD")
 {
-  gSystem->Load("libPWG0base");
-  TFile* file = TFile::Open("resC-resC/analysis_esd_raw.root");
-  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
-  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  loadlibs();
+
+  TFile* file = TFile::Open("fullA-simrec/MB2/analysis_esd_raw.root");
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
+  fdNdEtaAnalysis->LoadHistograms();
 
-  TFile* file2 = TFile::Open("resC-fullA/analysis_esd_raw.root");
-  dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta", "dndeta");
-  fdNdEtaAnalysis2->LoadHistograms("dndeta");
+  TFile* file2 = TFile::Open("fullA-sim/analysis_esd_raw.root");
+  dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
+  fdNdEtaAnalysis2->LoadHistograms();
 
   TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
   TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
@@ -224,7 +239,7 @@ void MisalignmentShowRawTrackPlots()
   peripheral->SetMarkerStyle(22);
   peripheral->SetMarkerColor(2);
 
-  TH2* tmp = new TH2F("tmp", ";p_{T} [GeV/c]      ;#frac{tracks residual misalignment}{tracks full misalignment}", 1, 0.1, 10, 1, 0.9, 1.3);
+  TH2* tmp = new TH2F("tmp", ";p_{T} [GeV/c]      ;#frac{tracks full misalignment during rec.}{tracks ideal misalignment during rec.}", 1, 0.1, 10, 1, 0.9, 1.3);
 
   tmp->SetStats(kFALSE);
   //tmp->GetXaxis()->SetNoExponent();
@@ -396,3 +411,68 @@ void vertexShift()
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 }
 
+void CompareRawTrackPlots(Float_t ptCut = 0.0)
+{
+  loadlibs();
+
+  const char* dirName = "fdNdEtaAnalysisESD";
+
+  TFile* file = TFile::Open("fullA-simrec/MB2/analysis_esd_raw.root");
+  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");
+  dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
+  fdNdEtaAnalysis2->LoadHistograms();
+
+  TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
+  TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
+
+  // normalize to number of events;
+  TH2* event1 = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
+  TH2* event2 = fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram();
+  Int_t event1Count = event1->Integral();
+  Int_t event2Count = event2->Integral();
+  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());
+
+  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);
+
+  gROOT->cd();
+
+  AliPWG0Helper::CreateDividedProjections(track1, track2);
+
+  new TCanvas; gROOT->FindObject("track1_yx_div_track2_yx")->Draw("COLZ");
+  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();
+}
+
+void MagnitudeOfCorrection(const char* fileName, const char* dirName = "dndeta", Float_t ptCut = 0.3)
+{
+  loadlibs();
+
+  TFile* file = TFile::Open(fileName);
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
+  fdNdEtaAnalysis->LoadHistograms();
+  fdNdEtaAnalysis->GetData()->PrintInfo(ptCut);
+}
+
index 447c939..fed908e 100644 (file)
@@ -1894,6 +1894,43 @@ Double_t FitPtFunc(Double_t *x, Double_t *par)
   }
 }
 
+void FitPtNew(const char* fileName = "TruePt14TeV.root")
+{
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libPWG0base");
+
+  TFile::Open(fileName);
+
+  TH1* genePt = (TH1*) gFile->Get("fHistPt");
+  genePt->Sumw2();
+
+  // normalize by bin width
+  for (Int_t x=1; x<genePt->GetNbinsX(); x++)
+    genePt->SetBinContent(x, genePt->GetBinContent(x) / genePt->GetBinWidth(x));
+
+  genePt->GetXaxis()->SetRangeUser(0.05, 2.0);
+
+  genePt->Scale(1.0 / genePt->Integral());
+
+  TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
+  //func->SetLineColor(2);
+  func->SetParameters(1, -1);
+
+  genePt->SetMarkerStyle(25);
+  genePt->SetTitle("");
+  genePt->SetStats(kFALSE);
+  genePt->GetYaxis()->SetRangeUser(1e-4, genePt->GetMaximum() * 1.2);
+  //func->Draw("SAME");
+
+  genePt->Fit(func, "0", "", 0.05, 1);
+
+  new TCanvas;
+  genePt->DrawCopy("P");
+  func->SetRange(0.02, 8);
+  func->DrawCopy("SAME");
+  gPad->SetLogy();
+}
+
 void FitPt(const char* fileName = "firstplots100k_truept.root")
 {
   gSystem->Load("libPWG0base");
@@ -1927,6 +1964,9 @@ void FitPt(const char* fileName = "firstplots100k_truept.root")
 
   TH1* genePt = gene->Project3D("z");*/
   TH1* genePt = (TH1*) gFile->Get("fdNdpTTrue");
+  if (!genePt)
+    genePt = (TH1*) gFile->Get("fHistPt");
   genePt->Sumw2();
 
   //genePt->Scale(1.0 / genePt->Integral());
@@ -1940,6 +1980,7 @@ void FitPt(const char* fileName = "firstplots100k_truept.root")
   genePt->GetXaxis()->SetRangeUser(0, 7.9);
   //genePt->GetYaxis()->SetTitle("a.u.");
 
+  //TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x*x)", 0.001, 100);
   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
   //func->SetLineColor(2);
   func->SetParameters(1, -1, 1, 1, 1);
index d16bb0c..fbdd39f 100644 (file)
@@ -90,3 +90,99 @@ void run(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, B
 
   mgr->StartAnalysis((aProof) ? "proof" : "local", chain);
 }
+
+void loadlibs()
+{
+  gSystem->Load("libTree");
+  gSystem->Load("libVMC");
+
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libPWG0base");
+}
+
+void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const char* dataOutput = "analysis_esd.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
+{
+  loadlibs();
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+  TFile::Open(correctionMapFile);
+  dNdEtaCorrection->LoadHistograms();
+
+  TFile* file = TFile::Open(dataInput);
+
+  if (!file)
+  {
+    cout << "Error. File not found" << endl;
+    return;
+  }
+
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
+  fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  TFile* file2 = TFile::Open(dataOutput, "RECREATE");
+  fdNdEtaAnalysis->SaveHistograms();
+
+  file->cd();
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco);
+  fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  file2->cd();
+  fdNdEtaAnalysis->SaveHistograms();
+
+  file->cd();
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle);
+  fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  file2->cd();
+  fdNdEtaAnalysis->SaveHistograms();
+
+  file->cd();
+  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+  fdNdEtaAnalysis->DrawHistograms(kTRUE);
+  file2->cd();
+  fdNdEtaAnalysis->SaveHistograms();
+}
+
+void* FinishAnalysis(const char* analysisFile = "analysis_esd_raw.root", const char* analysisDir = "fdNdEtaAnalysisESD", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", Bool_t useUncorrected = kFALSE, Bool_t simple = kFALSE)
+{
+  loadlibs();
+
+  TFile* file = TFile::Open(analysisFile);
+
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
+  fdNdEtaAnalysis->LoadHistograms();
+
+  if (correctionMapFile)
+  {
+    AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+    TFile::Open(correctionMapFile);
+    dNdEtaCorrection->LoadHistograms();
+
+    fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
+    //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kINEL);
+    //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
+  }
+  else
+    fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone);
+
+  fdNdEtaAnalysis->DrawHistograms(simple);
+
+  TH1* hist = fdNdEtaAnalysis->GetdNdEtaHistogram(1);
+  Int_t binLeft = hist->GetXaxis()->FindBin(-0.5);
+  Int_t binRight = hist->GetXaxis()->FindBin(0.5);
+  Float_t value1 = hist->Integral(binLeft, binRight);
+
+  hist = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
+  Float_t value2 = hist->Integral(binLeft, binRight);
+
+  if (value2 > 0)
+    printf("Ratio is %f, values are %f %f\n", value1 / value2, value1, value2);
+
+  return fdNdEtaAnalysis;
+}
index 2fce071..8d6fe92 100644 (file)
@@ -54,7 +54,7 @@ void runCorrection(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug =
 
   task = new AlidNdEtaCorrectionTask(option);
   task->SetTrackCuts(esdTrackCuts);
-  task->SetAnalysisMode(AlidNdEtaCorrectionTask::kSPD);
+  task->SetAnalysisMode(AlidNdEtaCorrectionTask::kTPC);
 
   mgr->AddTask(task);
 
@@ -77,7 +77,10 @@ void runCorrection(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug =
 
   // Enable debug printouts
   if (aDebug)
+  {
     mgr->SetDebugLevel(2);
+    AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kDebug+2);
+  }
 
   // Run analysis
   mgr->InitAnalysis();
index 7acd8e4..7c97279 100644 (file)
@@ -70,98 +70,3 @@ void rundNdEtaAnalysis(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC
   }
 }
 
-void loadlibs()
-{
-  gSystem->Load("libTree");
-  gSystem->Load("libVMC");
-
-  gSystem->Load("libSTEERBase");
-  gSystem->Load("libANALYSIS");
-  gSystem->Load("libPWG0base");
-}
-
-void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const char* dataOutput = "analysis_esd.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
-{
-  loadlibs();
-
-  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
-  TFile::Open(correctionMapFile);
-  dNdEtaCorrection->LoadHistograms();
-
-  TFile* file = TFile::Open(dataInput);
-
-  if (!file)
-  {
-    cout << "Error. File not found" << endl;
-    return;
-  }
-
-  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
-  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
-  fdNdEtaAnalysis->DrawHistograms(kTRUE);
-  TFile* file2 = TFile::Open(dataOutput, "RECREATE");
-  fdNdEtaAnalysis->SaveHistograms();
-
-  file->cd();
-  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
-  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco);
-  fdNdEtaAnalysis->DrawHistograms(kTRUE);
-  file2->cd();
-  fdNdEtaAnalysis->SaveHistograms();
-
-  file->cd();
-  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
-  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle);
-  fdNdEtaAnalysis->DrawHistograms(kTRUE);
-  file2->cd();
-  fdNdEtaAnalysis->SaveHistograms();
-
-  file->cd();
-  fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
-  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
-  fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
-  fdNdEtaAnalysis->DrawHistograms(kTRUE);
-  file2->cd();
-  fdNdEtaAnalysis->SaveHistograms();
-}
-
-void* FinishAnalysis(const char* analysisFile = "analysis_esd.root", const char* analysisDir = "dndeta", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", Bool_t useUncorrected = kFALSE, Bool_t simple = kFALSE)
-{
-  loadlibs();
-
-  TFile* file = TFile::Open(analysisFile);
-
-  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
-  fdNdEtaAnalysis->LoadHistograms();
-
-  if (correctionMapFile)
-  {
-    AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
-    TFile::Open(correctionMapFile);
-    dNdEtaCorrection->LoadHistograms();
-
-    //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.2, AlidNdEtaCorrection::kINEL);
-    fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kINEL);
-    //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
-  }
-  else
-    fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone);
-
-  fdNdEtaAnalysis->DrawHistograms(simple);
-
-  TH1* hist = fdNdEtaAnalysis->GetdNdEtaHistogram(1);
-  Int_t binLeft = hist->GetXaxis()->FindBin(-0.5);
-  Int_t binRight = hist->GetXaxis()->FindBin(0.5);
-  Float_t value1 = hist->Integral(binLeft, binRight);
-
-  hist = fdNdEtaAnalysis->GetdNdEtaHistogram(2);
-  Float_t value2 = hist->Integral(binLeft, binRight);
-
-  if (value2 > 0)
-    printf("Ratio is %f, values are %f %f\n", value1 / value2, value1, value2);
-
-  return fdNdEtaAnalysis;
-}
index 8e5eba0..559b225 100644 (file)
@@ -16,6 +16,7 @@
 #include <TGraph.h>
 #include <TLegend.h>
 #include <TLine.h>
+#include <TMath.h>
 
 #include <AliLog.h>
 #include <AliESD.h>
@@ -863,6 +864,7 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
 
   /*
 
+  gSystem->Load("libANALYSIS");
   gSystem->Load("libPWG0base");
   .L AliHighMultiplicitySelector.cxx+g
   x = new AliHighMultiplicitySelector();
@@ -974,6 +976,7 @@ void AliHighMultiplicitySelector::Ntrigger(Bool_t relative)
       if (downScale < 1)
         downScale = 1;
       Long64_t nTrigger = (Long64_t) (normalRate / downScale * lengthRun);
+      nTrigger = TMath::Nint(((Float_t) nTrigger) / 1000) * 1000;
 
       Printf("Normal rate is %f, downscale: %f, Simulating %lld triggers", normalRate, downScale, nTrigger);
       proj2->FillRandom(proj, nTrigger);