updated multiplicity analysis for the INEL>0 case
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Apr 2010 14:36:51 +0000 (14:36 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Apr 2010 14:36:51 +0000 (14:36 +0000)
PWG0/multiplicity/AliMultiplicityCorrection.cxx
PWG0/multiplicity/AliMultiplicityCorrection.h
PWG0/multiplicity/AliMultiplicityTask.cxx
PWG0/multiplicity/AliMultiplicityTask.h
PWG0/multiplicity/correct.C
PWG0/multiplicity/plots.C
PWG0/multiplicity/run.C

index c87e4f653ee68d4d6226599af015fb86310ff471..f58c6110de5d14f50216842bd8df3f9368fc758a 100644 (file)
 #include <TRandom.h>
 #include <TProfile.h>
 #include <TProfile2D.h>
+#include <AliLog.h>
 
 ClassImp(AliMultiplicityCorrection)
 
+// Defined where the efficiency drops below 1/3
+// |eta| < 1.4 --> -0.3 ... 0.8
+// |eta| < 1.3 --> -1.9 ... 2.4
+// |eta| < 1.0 --> -5.6 ... 6.1
+//Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.6, -1.9 };
+//Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] =   {  10.0,  6.1,  2.4 };
+Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.5, -1.9 };
+Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] =   {  10.0,  5.5,  2.4 };
+
 // These are the areas where the quality of the unfolding results are evaluated
 // Default defined here, call SetQualityRegions to change them
 // unit is in multiplicity (not in bin!)
-
 // SPD:   peak area - flat area - low stat area
 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {1,  20, 70};
 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {10, 65, 80};
@@ -96,7 +105,11 @@ AliMultiplicityCorrection::AliMultiplicityCorrection() :
   //
 
   for (Int_t i = 0; i < kESDHists; ++i)
+  {
     fMultiplicityESD[i] = 0;
+    fTriggeredEvents[i] = 0;
+    fNoVertexEvents[i] = 0;
+  }
 
   for (Int_t i = 0; i < kMCHists; ++i)
   {
@@ -128,6 +141,8 @@ AliMultiplicityCorrection* AliMultiplicityCorrection::Open(const char* fileName,
     return 0;
   }
   
+  Printf("AliMultiplicityCorrection::Open: Reading file %s", fileName);
+  
   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folderName, folderName);
   mult->LoadHistograms();
   
@@ -146,12 +161,14 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C
   fLastChi2MC(0),
   fLastChi2MCLimit(0),
   fLastChi2Residuals(0),
-  fRatioAverage(0)
+  fRatioAverage(0),
+  fVtxBegin(0),
+  fVtxEnd(0)
 {
   //
   // named constructor
   //
-
+  
   // do not add this hists to the directory
   Bool_t oldStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
@@ -176,10 +193,12 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C
 
   #define NBINNING 201, -0.5, 200.5
   
-  Double_t vtxRange[] = { 15, 6, 2 };
-
   for (Int_t i = 0; i < kESDHists; ++i)
-    fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 1, -vtxRange[i], vtxRange[i], NBINNING);
+  {
+    fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;measured multiplicity;Count", 1, fgVtxRangeBegin[i], fgVtxRangeEnd[i], NBINNING);
+    fTriggeredEvents[i] = new TH1F(Form("fTriggeredEvents%d", i), "fTriggeredEvents;measured multiplicity;Count", NBINNING);
+    fNoVertexEvents[i] = new TH1F(Form("fNoVertexEvents%d", i), "fNoVertexEvents;generated multiplicity;Count", NBINNING);
+  }
 
   for (Int_t i = 0; i < kMCHists; ++i)
   {
@@ -198,15 +217,134 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C
 
   for (Int_t i = 0; i < kCorrHists; ++i)
   {
-    fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 1, -vtxRange[i%3], vtxRange[i%3], NBINNING, NBINNING);
-    fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
+    fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;true multiplicity;measured multiplicity", 1, fgVtxRangeBegin[i%3], fgVtxRangeEnd[i%3], NBINNING, NBINNING);
+    fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;true multiplicity;Count", NBINNING);
   }
 
   TH1::AddDirectory(oldStatus);
 
   AliUnfolding::SetNbins(120, 120);
   AliUnfolding::SetSkipBinsBegin(1);
-  AliUnfolding::SetNormalizeInput(kTRUE);
+  //AliUnfolding::SetNormalizeInput(kTRUE);
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::Rebin2DY(TH2F*& hist, Int_t nBins, Double_t* newBins) const
+{
+  // 
+  // rebins the y axis of a two-dimensional histogram giving variable size binning (missing in ROOT v5/25/02)
+  //
+  
+  TH2F* temp = new TH2F(hist->GetName(), Form("%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle()), hist->GetNbinsX(), hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax(), nBins, newBins);
+  
+  for (Int_t x=0; x<=hist->GetNbinsX()+1; x++)
+    for (Int_t y=0; y<=hist->GetNbinsY()+1; y++)
+      temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetBinContent(x, y));
+  
+  for (Int_t x=0; x<=temp->GetNbinsX()+1; x++)
+    for (Int_t y=0; y<=temp->GetNbinsY()+1; y++)
+      temp->SetBinError(x, y, TMath::Sqrt(temp->GetBinContent(x, y)));
+      
+  delete hist;
+  hist = temp;
+} 
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::Rebin3DY(TH3F*& hist, Int_t nBins, Double_t* newBins) const
+{
+  // 
+  // rebins the y axis of a three-dimensional histogram giving variable size binning (missing in ROOT v5/25/02)
+  // this function is a mess - and it should have been Fons who should have gone through the pain of writing it! (JF)
+  //
+  
+  // construct variable size arrays for fixed size binning axes because TH3 lacks some constructors
+  Double_t* xBins = new Double_t[hist->GetNbinsX()+1];
+  Double_t* zBins = new Double_t[hist->GetNbinsZ()+1];
+  
+  for (Int_t x=1; x<=hist->GetNbinsX()+1; x++)
+  {
+    xBins[x-1] = hist->GetXaxis()->GetBinLowEdge(x);
+    //Printf("%d %f", x, xBins[x-1]);
+  }
+  
+  for (Int_t z=1; z<=hist->GetNbinsZ()+1; z++)
+  {
+    zBins[z-1] = hist->GetZaxis()->GetBinLowEdge(z);
+    //Printf("%d %f", y, yBins[y-1]);
+  }
+  
+  TH3F* temp = new TH3F(hist->GetName(), Form("%s;%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle(), hist->GetZaxis()->GetTitle()), hist->GetNbinsX(), xBins, nBins, newBins, hist->GetNbinsZ(), zBins);
+  
+  for (Int_t x=0; x<=hist->GetNbinsX()+1; x++)
+    for (Int_t y=0; y<=hist->GetNbinsY()+1; y++)
+      for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++)
+        temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z), hist->GetBinContent(x, y, z));
+  
+  for (Int_t x=0; x<=temp->GetNbinsX()+1; x++)
+    for (Int_t y=0; y<=temp->GetNbinsY()+1; y++)
+      for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++)
+        temp->SetBinError(x, y, z, TMath::Sqrt(temp->GetBinContent(x, y, z)));
+      
+  delete[] xBins;
+  delete[] zBins;
+      
+  delete hist;
+  hist = temp;
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::RebinGenerated(Int_t nBins05, Double_t* newBins05, Int_t nBins10, Double_t* newBins10, Int_t nBins13, Double_t* newBins13)
+{
+  //
+  // Rebins the (and only the) generated multiplicity axis 
+  //
+  
+  Printf("Rebinning generated-multiplicity axis...");
+
+  // do not add this hists to the directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+  
+  if (kESDHists != 3)
+    AliFatal("This function only works for three ESD hists!");
+  
+  for (Int_t i = 0; i < kESDHists; ++i)
+  {
+    Int_t nBins = -1;
+    Double_t* newBins = 0;
+    
+    switch (i)
+    {
+      case 0:
+        nBins = nBins05;
+        newBins = newBins05;
+        break;
+      case 1:
+        nBins = nBins10;
+        newBins = newBins10;
+        break;
+      case 2:
+        nBins = nBins13;
+        newBins = newBins13;
+        break;
+    }
+  
+    // 1D
+    // TODO mem leak
+    fNoVertexEvents[i] = (TH1F*) fNoVertexEvents[i]->Rebin(nBins, fNoVertexEvents[i]->GetName(), newBins);
+    fMultiplicityESDCorrected[i] = (TH1F*) fMultiplicityESDCorrected[i]->Rebin(nBins, fMultiplicityESDCorrected[i]->GetName(), newBins);
+  
+    // 2D
+    Rebin2DY(fMultiplicityVtx[i], nBins, newBins);
+    Rebin2DY(fMultiplicityMB[i], nBins, newBins);
+    Rebin2DY(fMultiplicityINEL[i], nBins, newBins);
+    Rebin2DY(fMultiplicityNSD[i], nBins, newBins);
+
+    // 3D
+    Rebin3DY(fCorrelation[i], nBins, newBins);
+  }
+
+  TH1::AddDirectory(oldStatus);
 }
 
 //____________________________________________________________________
@@ -223,6 +361,14 @@ AliMultiplicityCorrection::~AliMultiplicityCorrection()
     if (fMultiplicityESD[i])
       delete fMultiplicityESD[i];
     fMultiplicityESD[i] = 0;
+    
+    if (fTriggeredEvents[i])
+      delete fTriggeredEvents[i];
+    fTriggeredEvents[i]= 0;
+  
+    if (fNoVertexEvents[i])
+      delete fNoVertexEvents[i];
+    fNoVertexEvents[i]= 0;
   }
 
   for (Int_t i = 0; i < kMCHists; ++i)
@@ -273,7 +419,7 @@ Long64_t AliMultiplicityCorrection::Merge(const TCollection* list)
   TObject* obj;
 
   // collections of all histograms
-  TList collections[kESDHists+kMCHists*4+kCorrHists*2];
+  TList collections[3*kESDHists+kMCHists*4+kCorrHists*2];
 
   Int_t count = 0;
   while ((obj = iter->Next())) {
@@ -283,41 +429,49 @@ Long64_t AliMultiplicityCorrection::Merge(const TCollection* list)
       continue;
 
     for (Int_t i = 0; i < kESDHists; ++i)
+    {
       collections[i].Add(entry->fMultiplicityESD[i]);
+      collections[kESDHists+i].Add(entry->fTriggeredEvents[i]);
+      collections[kESDHists*2+i].Add(entry->fNoVertexEvents[i]);
+    }
 
     for (Int_t i = 0; i < kMCHists; ++i)
     {
-      collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
-      collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
-      collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
-      collections[kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
+      collections[3*kESDHists+i].Add(entry->fMultiplicityVtx[i]);
+      collections[3*kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
+      collections[3*kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
+      collections[3*kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
     }
 
     for (Int_t i = 0; i < kCorrHists; ++i)
-      collections[kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
+      collections[3*kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
 
     for (Int_t i = 0; i < kCorrHists; ++i)
-      collections[kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
+      collections[3*kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
 
     count++;
   }
 
   for (Int_t i = 0; i < kESDHists; ++i)
+  {
     fMultiplicityESD[i]->Merge(&collections[i]);
-
+    fTriggeredEvents[i]->Merge(&collections[kESDHists+i]);
+    fNoVertexEvents[i]->Merge(&collections[2*kESDHists+i]);
+  }
+  
   for (Int_t i = 0; i < kMCHists; ++i)
   {
-    fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
-    fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
-    fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
-    fMultiplicityNSD[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
+    fMultiplicityVtx[i]->Merge(&collections[3*kESDHists+i]);
+    fMultiplicityMB[i]->Merge(&collections[3*kESDHists+kMCHists+i]);
+    fMultiplicityINEL[i]->Merge(&collections[3*kESDHists+kMCHists*2+i]);
+    fMultiplicityNSD[i]->Merge(&collections[3*kESDHists+kMCHists*3+i]);
   }
 
   for (Int_t i = 0; i < kCorrHists; ++i)
-    fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*4+i]);
+    fCorrelation[i]->Merge(&collections[3*kESDHists+kMCHists*4+i]);
 
   for (Int_t i = 0; i < kCorrHists; ++i)
-    fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*4+kCorrHists+i]);
+    fMultiplicityESDCorrected[i]->Merge(&collections[3*kESDHists+kMCHists*4+kCorrHists+i]);
 
   delete iter;
 
@@ -342,9 +496,15 @@ Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
   TList oldObjects;
   oldObjects.SetOwner(1);
   for (Int_t i = 0; i < kESDHists; ++i)
+  {
     if (fMultiplicityESD[i])
       oldObjects.Add(fMultiplicityESD[i]);
-
+    if (fTriggeredEvents[i])
+      oldObjects.Add(fTriggeredEvents[i]);
+    if (fNoVertexEvents[i])
+      oldObjects.Add(fNoVertexEvents[i]);
+  }
+  
   for (Int_t i = 0; i < kMCHists; ++i)
   {
     if (fMultiplicityVtx[i])
@@ -368,7 +528,9 @@ Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
   for (Int_t i = 0; i < kESDHists; ++i)
   {
     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
-    if (!fMultiplicityESD[i])
+    fTriggeredEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fTriggeredEvents[i]->GetName()));
+    fNoVertexEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fNoVertexEvents[i]->GetName()));
+    if (!fMultiplicityESD[i] || !fTriggeredEvents[i] || !fNoVertexEvents[i])
       success = kFALSE;
   }
 
@@ -414,11 +576,17 @@ void AliMultiplicityCorrection::SaveHistograms(const char* dir)
   gDirectory->cd(dir);
 
   for (Int_t i = 0; i < kESDHists; ++i)
+  {
     if (fMultiplicityESD[i])
     {
       fMultiplicityESD[i]->Write();
       fMultiplicityESD[i]->ProjectionY(Form("%s_px", fMultiplicityESD[i]->GetName()), 1, fMultiplicityESD[i]->GetNbinsX())->Write();
     }
+    if (fTriggeredEvents[i])
+      fTriggeredEvents[i]->Write();
+    if (fNoVertexEvents[i])
+      fNoVertexEvents[i]->Write();
+  }
 
   for (Int_t i = 0; i < kMCHists; ++i)
   {
@@ -504,6 +672,35 @@ void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_
   fMultiplicityESD[2]->Fill(vtx, measured14);
 }
 
+//____________________________________________________________________
+void AliMultiplicityCorrection::FillTriggeredEvent(Int_t measured05, Int_t measured10, Int_t measured14)
+{
+  //
+  // fills raw distribution of triggered events
+  //
+  
+  fTriggeredEvents[0]->Fill(measured05);
+  fTriggeredEvents[1]->Fill(measured10);
+  fTriggeredEvents[2]->Fill(measured14);
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::FillNoVertexEvent(Float_t vtx, Bool_t vertexReconstructed, Int_t generated05, Int_t generated10, Int_t generated14, Int_t measured05, Int_t measured10, Int_t measured14)
+{
+  //
+  // fills raw distribution of triggered events
+  //
+  
+  if (vtx > fgVtxRangeBegin[0] && vtx < fgVtxRangeEnd[0] && (!vertexReconstructed || measured05 == 0))
+    fNoVertexEvents[0]->Fill(generated05);
+    
+  if (vtx > fgVtxRangeBegin[1] && vtx < fgVtxRangeEnd[1] && (!vertexReconstructed || measured10 == 0))
+    fNoVertexEvents[1]->Fill(generated10);
+    
+  if (vtx > fgVtxRangeBegin[2] && vtx < fgVtxRangeEnd[2] && (!vertexReconstructed || measured14 == 0))
+    fNoVertexEvents[2]->Fill(generated14);
+}
+
 //____________________________________________________________________
 void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured14)
 {
@@ -534,7 +731,14 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
   fMultiplicityESDCorrected[correlationID]->Sumw2();
 
   // project without under/overflow bins
-  fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
+  Int_t begin = 1;
+  Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins();
+  if (fVtxEnd > fVtxBegin)
+  {
+    begin = fVtxBegin;
+    end = fVtxEnd;
+  }
+  fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", begin, end);
   fCurrentESD->Sumw2();
 
   // empty under/overflow bins in x, otherwise Project3D takes them into account
@@ -548,6 +752,9 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
     }
   }
 
+  if (fVtxEnd > fVtxBegin)
+    hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd);
+  
   fCurrentCorrelation = (TH2*) hist->Project3D("zy");
   fCurrentCorrelation->Sumw2();
   
@@ -582,15 +789,21 @@ TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventT
   // calculates efficiency for given event type
   //
   
+  TString name1;
+  name1.Form("divisor%d", inputRange);
+  
+  TString name2;
+  name2.Form("CurrentEfficiency%d", inputRange);
+  
   TH1* divisor = 0;
   switch (eventType)
   {
     case kTrVtx : break;
-    case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
-    case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
-    case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY("divisor", 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
+    case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY(name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
+    case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
+    case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY(name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
   }
-  TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
+  TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY(name2, 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
   
   if (eventType == kTrVtx)
   {
@@ -604,20 +817,34 @@ TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventT
 }
 
 //____________________________________________________________________
-TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange)
+TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange, EventType eventType)
 {
   //
   // calculates efficiency for given event type
   //
 
-  TH1* divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e");
-  TH1* eff = fMultiplicityMB[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
+  TString name1;
+  name1.Form("divisor%d", inputRange);
+  
+  TString name2;
+  name2.Form("CurrentEfficiency%d", inputRange);
+
+  TH1* divisor = 0;
+  switch (eventType)
+  {
+    case kTrVtx : AliFatal("Not supported!"); break;
+    case kMB: divisor =   fMultiplicityMB[inputRange]->ProjectionY  (name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
+    case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
+    case kNSD: divisor =  fMultiplicityNSD[inputRange]->ProjectionY (name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
+  }
+  TH1* eff = fMultiplicityMB[inputRange]->ProjectionY(name2, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
+  
   eff->Divide(eff, divisor, 1, 1, "B");
   return eff;
 }
 
 //____________________________________________________________________
-Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
+Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t check, TH1* initialConditions, Bool_t errorAsBias)
 {
   //
   // correct spectrum using minuit chi2 method
@@ -627,23 +854,106 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha
 
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   
-  AliUnfolding::SetCreateOverflowBin(5);
+  //AliUnfolding::SetCreateOverflowBin(5);
   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+  AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1);
+  
+  // use here only vtx efficiency (to MB sample) which is always needed if we use the 0 bin
+  SetupCurrentHists(inputRange, fullPhaseSpace, (eventType == kTrVtx) ? kTrVtx : kMB);
+  
+  // TODO set errors on measured with 0.5 * TMath::ChisquareQuantile(0.1, 20000)
+  // see PDG: Statistics / Poission or binomial data / Eq. 32.49a/b in 2004 edition
+  
+  Calculate0Bin(inputRange, eventType, zeroBinEvents);
+
+  Int_t resultCode = -1;
+  if (errorAsBias == kFALSE)
+  {
+    resultCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
+  }
+  else
+  {
+    resultCode = AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]);
+  }
+  
+  // HACK store new vertex reco efficiency for bin 0, changing number of events with trigger and vertex in MC map
+  if (zeroBinEvents > 0)
+  {
+    Printf("WARNING: Stored vertex reco efficiency from unfolding for bin 0.");
+    fMultiplicityVtx[inputRange]->SetBinContent(1, 1, fMultiplicityMB[inputRange]->GetBinContent(1, 1) * fCurrentEfficiency->GetBinContent(1));
+  }
+  
+  // correct for the trigger bias if requested
+  if (eventType > kMB)
+  {
+    Printf("Applying trigger efficiency");
+    TH1* eff = GetTriggerEfficiency(inputRange, eventType);
+    for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); i++)
+    {
+      fMultiplicityESDCorrected[correlationID]->SetBinContent(i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i) / eff->GetBinContent(i));
+      fMultiplicityESDCorrected[correlationID]->SetBinError(i, fMultiplicityESDCorrected[correlationID]->GetBinError(i) / eff->GetBinContent(i));
+    }
+  }
+  
+  return resultCode;
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::Calculate0Bin(Int_t inputRange, EventType eventType, Int_t zeroBinEvents)
+{
+  // fills the 0 bin
   
-  if (!initialConditions)
+  if (eventType == kTrVtx)
+    return;
+  
+  Double_t fractionEventsInVertexRange = fMultiplicityESD[inputRange]->Integral(1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityESD[inputRange]->Integral(0, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins() + 1);
+  
+  // difference of fraction that is inside the considered range between triggered events and events with vertex
+  // Extension to NSD not needed, INEL and NSD vertex distributions are nature-given and unbiased!
+  Double_t differenceVtxDist = (fMultiplicityINEL[inputRange]->Integral(1, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityINEL[inputRange]->Integral(0, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins() + 1)) / (fMultiplicityVtx[inputRange]->Integral(1, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityVtx[inputRange]->Integral(0, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins() + 1));
+
+  Printf("Enabling 0 bin estimate for triggered events without vertex.");
+  Printf("  Events in 0 bin: %d", zeroBinEvents);
+  Printf("  Fraction in range: %.1f%%", fractionEventsInVertexRange * 100);
+  Printf("  Difference Vtx Dist: %f", differenceVtxDist);
+  
+  AliUnfolding::SetNotFoundEvents(zeroBinEvents * fractionEventsInVertexRange * differenceVtxDist);
+}
+
+//____________________________________________________________________
+void AliMultiplicityCorrection::FixTriggerEfficiencies(Int_t start)
+{
+  //
+  // sets trigger and vertex efficiencies to 1 for large multiplicities where no event was simulated
+  //
+  
+  for (Int_t etaRange = 0; etaRange < kMCHists; etaRange++)
   {
-    initialConditions = (TH1*) fCurrentESD->Clone("initialConditions");
-    initialConditions->Scale(1.0 / initialConditions->Integral());
-    if (!check)
+    for (Int_t i = 1; i <= fMultiplicityINEL[etaRange]->GetNbinsY(); i++)
     {
-      // set minimum value to prevent MINUIT just staying in the small value
-      for (Int_t i=1; i<=initialConditions->GetNbinsX(); i++)
-        initialConditions->SetBinContent(i, TMath::Max(initialConditions->GetBinContent(i), 1e-3));
+      if (fMultiplicityINEL[etaRange]->GetYaxis()->GetBinCenter(i) < start)
+        continue;
+        
+      if (fMultiplicityINEL[etaRange]->GetBinContent(1, i) > 0)
+        continue;
+        
+      fMultiplicityINEL[etaRange]->SetBinContent(1, i, 1);
+      fMultiplicityNSD[etaRange]->SetBinContent(1, i, 1);
+      fMultiplicityMB[etaRange]->SetBinContent(1, i, 1);
+      fMultiplicityVtx[etaRange]->SetBinContent(1, i, 1);
     }
   }
+}
 
-  return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
+//____________________________________________________________________
+Float_t AliMultiplicityCorrection::GetFraction0Generated(Int_t inputRange)
+{
+  //
+  // returns the fraction of events that have 0 generated particles in the given range among all events without vertex OR 0 tracklets
+  //
+  
+  TH1* multMB = GetNoVertexEvents(inputRange);
+  return multMB->GetBinContent(1) / multMB->Integral();
 }
 
 //____________________________________________________________________
@@ -723,12 +1033,9 @@ void AliMultiplicityCorrection::DrawHistograms()
 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple, EventType eventType)
 {
   // draw comparison plots
-
-
-  //mcHist->Rebin(2);
-
+  
   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
-
+  
   TString tmpStr;
   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
 
@@ -737,111 +1044,54 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     printf("ERROR. Unfolded histogram is empty\n");
     return;
   }
-
-  //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
-  fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
-  fCurrentESD->Sumw2();
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
-
-  // normalize unfolded result to 1
-  fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
-
-  //fCurrentESD->Scale(mcHist->Integral(2, 200));
-
-  //new TCanvas;
-  /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
-  ratio->Divide(mcHist);
-  ratio->Draw("HIST");
-  ratio->Fit("pol0", "W0", "", 20, 120);
-  Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
-  delete ratio;
-  mcHist->Scale(scalingFactor);*/
-
-  // find bin with <= 50 entries. this is used as limit for the chi2 calculation
-  Int_t mcBinLimit = 0;
-  for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
+  
+  Int_t begin = 1;
+  Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins();
+  if (fVtxEnd > fVtxBegin)
   {
-    if (mcHist->GetBinContent(i) > 50)
-    {
-      mcBinLimit = i;
-    }
-    else
-      break;
+    begin = fVtxBegin;
+    end = fVtxEnd;
   }
-  Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
-  
-  // scale to 1
-  mcHist->Sumw2();
-  if (mcHist->Integral() > 0)
-    mcHist->Scale(1.0 / mcHist->Integral());
-
-  // calculate residual
+  fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", begin, end);
+  fCurrentESD->Sumw2();
 
-  // for that we convolute the response matrix with the gathered result
-  // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
-  TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
-  
-  // undo trigger, vertex efficiency correction
-  fCurrentEfficiency = GetEfficiency(inputRange, eventType);
-  tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
-  
-  TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
-  TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
-  if (convolutedProj->Integral() > 0)
+  mcHist->Sumw2();
+  Int_t mcMax = 0;
+  for (Int_t i=5; i<=mcHist->GetNbinsX(); ++i)
   {
-    convolutedProj->Scale(1.0 / convolutedProj->Integral());
+    if (mcHist->GetBinContent(i) > 0)
+      mcMax = mcHist->GetXaxis()->GetBinCenter(i) + 2;
   }
-  else
-    printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
-
-  TH1* residual = (TH1*) convolutedProj->Clone("residual");
-  residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
-
-  residual->Add(fCurrentESD, -1);
-  //residual->Divide(residual, fCurrentESD, 1, 1, "B");
+  if (mcMax == 0)
+  {
+    for (Int_t i=5; i<=fMultiplicityESDCorrected[esdCorrId]->GetNbinsX(); ++i)
+      if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) > 1)
+        mcMax = fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->GetBinCenter(i) + 2;
+  }  
+  Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcMax);
+  // calculate residual
+  Float_t tmp;
+  TH1* convolutedProj = (TH1*) GetConvoluted(esdCorrId, eventType)->Clone("convolutedProj");
+  TH1* residual = GetResiduals(esdCorrId, eventType, tmp);
 
   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
 
-  // find bin limit
-  Int_t lastBin = 0;
-  for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
-  {
-    if (fCurrentESD->GetBinContent(i) <= 5)
-    {
-      lastBin = i;
-      break;
-    }
-  }
-  
-  // TODO fix errors
   Float_t chi2 = 0;
-  for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
+  for (Int_t i=1; i<=TMath::Min(residual->GetNbinsX(), 75); ++i)
   {
-    if (fCurrentESD->GetBinError(i) > 0)
-    {
-      Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
-      if (i > 1 && i <= lastBin)
-        chi2 += value * value;
-      Printf("%d --> %f (%f)", i, value * value, chi2);
-      residual->SetBinContent(i, value);
-      residualHist->Fill(value);
-    }
-    else
-    {
-      //printf("Residual bin %d set to 0\n", i);
-      residual->SetBinContent(i, 0);
-    }
+    Float_t value = residual->GetBinContent(i);
+    // TODO has to get a parameter (used in Chi2Function, GetResiduals, and here)
+    if (i > 1)
+      chi2 += value * value;
+    Printf("%d --> %f (%f)", i, value * value, chi2);
+    residualHist->Fill(value);
     convolutedProj->SetBinError(i, 0);
-    residual->SetBinError(i, 0);
   }
   fLastChi2Residuals = chi2;
 
-  new TCanvas; residualHist->DrawCopy();
-
-  //residualHist->Fit("gaus", "N");
-  //delete residualHist;
+  //new TCanvas; residualHist->DrawCopy();
 
-  printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, lastBin);
+  printf("Difference (Residuals) is %f\n", fLastChi2Residuals);
 
   TCanvas* canvas1 = 0;
   if (simple)
@@ -863,31 +1113,33 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   canvas1->cd(1)->SetLeftMargin(0.12);
   canvas1->cd(1)->SetBottomMargin(0.12);
   TH1* proj = (TH1*) mcHist->Clone("proj");
+  if (proj->GetEntries() > 0)
+    AliPWG0Helper::NormalizeToBinWidth(proj);
 
-  // normalize without 0 bin
-  proj->Scale(1.0 / proj->Integral(2, proj->GetNbinsX()));
-  Printf("Normalized without 0 bin!");
-  proj->GetXaxis()->SetRangeUser(0, 200);
+  proj->GetXaxis()->SetRangeUser(0, mcMax);
   proj->GetYaxis()->SetTitleOffset(1.4);
-  //proj->SetLabelSize(0.05, "xy");
-  //proj->SetTitleSize(0.05, "xy");
   proj->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
   proj->SetStats(kFALSE);
 
+  fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, mcMax);
+  fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetRangeUser(0.1, fMultiplicityESDCorrected[esdCorrId]->GetMaximum() * 1.5);
+  fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetTitleOffset(1.4);
+  fMultiplicityESDCorrected[esdCorrId]->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
+  fMultiplicityESDCorrected[esdCorrId]->SetStats(kFALSE);
+  
   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
   fMultiplicityESDCorrected[esdCorrId]->SetMarkerColor(2);
-  //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(5);
-  // normalize without 0 bin
-  fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(2, fMultiplicityESDCorrected[esdCorrId]->GetNbinsX()));
-  Printf("Normalized without 0 bin!");
 
+  TH1* esdCorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("esdCorrected");
+  AliPWG0Helper::NormalizeToBinWidth(esdCorrected);
+  
   if (proj->GetEntries() > 0) {
     proj->DrawCopy("HIST");
-    fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
+    esdCorrected->DrawCopy("SAME HIST E");
   }
   else
-    fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
-
+    esdCorrected->DrawCopy("HIST E");
+    
   gPad->SetLogy();
 
   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
@@ -896,7 +1148,6 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   legend->SetFillColor(0);
   legend->SetTextSize(0.04);
   legend->Draw();
-  // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
 
   canvas1->cd(2);
   canvas1->cd(2)->SetGridx();
@@ -907,20 +1158,15 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   canvas1->cd(2)->SetBottomMargin(0.12);
 
   gPad->SetLogy();
-  fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
-  //fCurrentESD->SetLineColor(2);
+  fCurrentESD->GetXaxis()->SetRangeUser(0, mcMax);
   fCurrentESD->SetTitle(Form(";Measured multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
   fCurrentESD->SetStats(kFALSE);
   fCurrentESD->GetYaxis()->SetTitleOffset(1.4);
-  //fCurrentESD->SetLabelSize(0.05, "xy");
-  //fCurrentESD->SetTitleSize(0.05, "xy");
   fCurrentESD->DrawCopy("HIST E");
 
   convolutedProj->SetLineColor(2);
   convolutedProj->SetMarkerColor(2);
   convolutedProj->SetMarkerStyle(5);
-  //proj2->SetMarkerColor(2);
-  //proj2->SetMarkerStyle(5);
   convolutedProj->DrawCopy("HIST SAME P");
 
   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
@@ -930,224 +1176,28 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   legend->SetTextSize(0.04);
   legend->Draw();
 
-  //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
-  //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
-
-  /*Float_t chi2 = 0;
-  Float_t chi = 0;
-  fLastChi2MCLimit = 0;
-  Int_t limit = 0;
-  for (Int_t i=2; i<=200; ++i)
-  {
-    if (proj->GetBinContent(i) != 0)
-    {
-      Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
-      chi2 += value * value;
-      chi += TMath::Abs(value);
-
-      //printf("%d %f\n", i, chi);
-
-      if (chi2 < 0.2)
-        fLastChi2MCLimit = i;
-
-      if (chi < 3)
-        limit = i;
-
-    }
-  }*/
-
-  /*chi2 = 0;
-  Float_t chi = 0;
-  Int_t limit = 0;
-  for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
-  {
-    if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
-    {
-      Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
-      if (value > 1e8)
-        value = 1e8; //prevent arithmetic exception
-      else if (value < -1e8)
-        value = -1e8;
-      if (i > 1)
-      {
-        chi2 += value * value;
-        chi += TMath::Abs(value);
-      }
-      diffMCUnfolded->SetBinContent(i, value);
-    }
-    else
-    {
-      //printf("diffMCUnfolded bin %d set to 0\n", i);
-      diffMCUnfolded->SetBinContent(i, 0);
-    }
-    if (chi2 < 1000)
-      fLastChi2MCLimit = i;
-    if (chi < 1000)
-      limit = i;
-    if (i == 150)
-      fLastChi2MC = chi2;
-  }
-
-  printf("limits %d %d\n", fLastChi2MCLimit, limit);
-  fLastChi2MCLimit = limit;
-
-  printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
-
   if (!simple)
   {
-    /*canvas1->cd(3);
-
-    diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
-    //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
-    diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
-    diffMCUnfolded->DrawCopy("HIST");
-
-    TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
-    for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
-      fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
-
-    //new TCanvas; fluctuation->DrawCopy();
-    delete fluctuation;*/
-
-    /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
-    legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
-    legend->AddEntry(fMultiplicityMC, "MC");
-    legend->AddEntry(fMultiplicityESD, "ESD");
-    legend->Draw();*/
-
     canvas1->cd(4);
     residual->GetYaxis()->SetRangeUser(-5, 5);
-    residual->GetXaxis()->SetRangeUser(0, 200);
+    residual->GetXaxis()->SetRangeUser(0, mcMax);
+    residual->SetStats(kFALSE);
     residual->DrawCopy();
 
     canvas1->cd(5);
-
     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
     ratio->Divide(mcHist);
     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
-    ratio->GetXaxis()->SetRangeUser(0, 200);
+    ratio->GetXaxis()->SetRangeUser(0, mcMax);
     ratio->SetStats(kFALSE);
     ratio->Draw("HIST");
 
-    Double_t ratioChi2 = 0;
-    fRatioAverage = 0;
-    fLastChi2MCLimit = 0;
-    for (Int_t i=2; i<=150; ++i)
-    {
-      Float_t value = ratio->GetBinContent(i) - 1;
-      if (value > 1e8)
-        value = 1e8; //prevent arithmetic exception
-      else if (value < -1e8)
-        value = -1e8;
-
-      ratioChi2 += value * value;
-      fRatioAverage += TMath::Abs(value);
-
-      if (ratioChi2 < 0.1)
-        fLastChi2MCLimit = i;
-    }
-    fRatioAverage /= 149;
-
-    printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
-    // TODO FAKE
-    fLastChi2MC = ratioChi2;
-
-    // FFT of ratio
-    canvas1->cd(6);
-    const Int_t kFFT = 128;
-    Double_t fftReal[kFFT];
-    Double_t fftImag[kFFT];
-
-    for (Int_t i=0; i<kFFT; ++i)
-    {
-      fftReal[i] = ratio->GetBinContent(i+1+10);
-      // test: ;-)
-      //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
-      fftImag[i] = 0;
-    }
-
-    FFT(-1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
-
-    TH1* fftResult = (TH1*) ratio->Clone("fftResult");
-    fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
-    TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
-    TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
-    fftResult->Reset();
-    fftResult2->Reset();
-    fftResult3->Reset();
-    fftResult->GetYaxis()->UnZoom();
-    fftResult2->GetYaxis()->UnZoom();
-    fftResult3->GetYaxis()->UnZoom();
-
-    //Printf("AFTER FFT");
-    for (Int_t i=0; i<kFFT; ++i)
-    {
-      //Printf("%d: %f", i, fftReal[i]);
-      fftResult->SetBinContent(i+1, fftReal[i]);
-      /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
-      {
-        Printf("Nulled %d", i);
-        fftReal[i] = 0;
-      }*/
-    }
-
-    fftResult->SetLineColor(1);
-    fftResult->DrawCopy();
-
-
-    //Printf("IMAG");
-    for (Int_t i=0; i<kFFT; ++i)
-    {
-      //Printf("%d: %f", i, fftImag[i]);
-      fftResult2->SetBinContent(i+1, fftImag[i]);
-
-      fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
-    }
-
-    fftResult2->SetLineColor(2);
-    fftResult2->DrawCopy("SAME");
-
-    fftResult3->SetLineColor(4);
-    fftResult3->DrawCopy("SAME");
-
-    for (Int_t i=1; i<kFFT - 1; ++i)
-    {
-      if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
-      {
-        fftReal[i-1] = 0;
-        fftReal[i] = 0;
-        fftReal[i+1] = 0;
-        fftImag[i-1] = 0;
-        fftImag[i] = 0;
-        fftImag[i+1] = 0;
-        //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
-        //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
-        //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
-      }
-    }
-
-    FFT(1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
-
-    TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
-    fftResult4->Reset();
-
-    //Printf("AFTER BACK-TRAFO");
-    for (Int_t i=0; i<kFFT; ++i)
-    {
-      //Printf("%d: %f", i, fftReal[i]);
-      fftResult4->SetBinContent(i+1+10, fftReal[i]);
-    }
-
-    canvas1->cd(5);
-    fftResult4->SetLineColor(4);
-    fftResult4->DrawCopy("SAME");
-
     // plot (MC - Unfolded) / error (MC)
     canvas1->cd(3);
 
     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
-    diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
+    diffMCUnfolded2->Add(esdCorrected, -1);
 
     Int_t ndfQual[kQualityRegions];
     for (Int_t region=0; region<kQualityRegions; ++region)
@@ -1156,7 +1206,6 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
       ndfQual[region] = 0;
     }
 
-
     Double_t newChi2 = 0;
     Double_t newChi2Limit150 = 0;
     Int_t ndf = 0;
@@ -1167,7 +1216,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
       {
         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
         newChi2 += value * value;
-        if (i > 1 && i <= mcBinLimit)
+        if (i > 1 && i <= mcMax)
           newChi2Limit150 += value * value;
         ++ndf;
 
@@ -1190,24 +1239,29 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
     }
 
-    if (mcBinLimit > 1)
+    if (mcMax > 1)
     {
-      // TODO another FAKE
-      fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
-      Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
+      fLastChi2MC = newChi2Limit150 / (mcMax - 1);
+      Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcMax, newChi2Limit150, mcMax - 1, fLastChi2MC);
     }
     else
       fLastChi2MC = -1;
 
     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
 
-    diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
+    diffMCUnfolded2->SetTitle("#chi^{2};true multiplicity;(MC - Unfolded) / e(MC)");
     diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5);
-    diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
+    diffMCUnfolded2->GetXaxis()->SetRangeUser(0, mcMax);
     diffMCUnfolded2->DrawCopy("HIST");
+    
+    canvas1->cd(6);
+    // draw penalty factor
+    
+    TH1* penalty = AliUnfolding::GetPenaltyPlot(fMultiplicityESDCorrected[esdCorrId]);
+    penalty->SetStats(0);
+    penalty->GetXaxis()->SetRangeUser(0, mcMax);
+    penalty->DrawCopy("HIST");
   }
-
-  canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
 }
 
 //____________________________________________________________________
@@ -1391,7 +1445,7 @@ TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
 }
 
 //____________________________________________________________________
-TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo)
+TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo)
 {
   //
   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
@@ -1409,10 +1463,14 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType
   gRandom->SetSeed(0);
   
   if (methodType == AliUnfolding::kChi2Minimization)
-    AliUnfolding::SetCreateOverflowBin(5);
+  {
+    Calculate0Bin(inputRange, eventType, zeroBinEvents);
+    AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1);
+  }
+  
   AliUnfolding::SetUnfoldingMethod(methodType);
 
-  const Int_t kErrorIterations = 150;
+  const Int_t kErrorIterations = 20;
 
   TH1* maxError = 0;
   TH1* firstResult = 0;
@@ -1493,7 +1551,10 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType
       Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result);
 
       if (returnCode != 0)
-        return 0;
+      {
+       n--;
+       continue;
+      }
     }
 
     // normalize
@@ -1628,17 +1689,21 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType
     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
 
   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
-    fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
+    fMultiplicityESDCorrected[correlationID]->SetBinError(i, standardDeviation->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
 
   return standardDeviation;
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
+void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Int_t determineError)
 {
   //
   // correct spectrum using bayesian method
   //
+  // determineError: 
+  //   0 = no errors
+  //   1 = from randomizing
+  //   2 = with UnfoldGetBias
 
   // initialize seed with current time
   gRandom->SetSeed(0);
@@ -1676,10 +1741,19 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   AliUnfolding::SetBayesianParameters(regPar, nIterations);
   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
-  if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
+  
+  if (determineError <= 1)
+  {
+    if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
+      return;
+  }
+  else if (determineError == 2)
+  {
+    AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]);
     return;
+  }  
 
-  if (!determineError)
+  if (determineError == 0)
   {
     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
     return;
@@ -1692,7 +1766,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   const Int_t kErrorIterations = 20;
 
-  printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
+  Printf("Spectrum unfolded. Determining error (%d iterations)...", kErrorIterations);
 
   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
   TH1* resultArray[kErrorIterations+1];
@@ -1710,9 +1784,10 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
     TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
     result2->Reset();
     if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0)
-      return;
-
-    result2->Scale(1.0 / result2->Integral());
+    {
+      n--;
+      continue;
+    }
 
     resultArray[n+1] = result2;
   }
@@ -1724,8 +1799,12 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
   for (Int_t n=0; n<kErrorIterations; ++n)
     delete resultArray[n+1];
 
+  Printf("Comparing bias and error:");
   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
+  {
+    Printf("Bin %d: Content: %f Error: %f Bias: %f", i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i), error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i), fMultiplicityESDCorrected[correlationID]->GetBinError(i));
     fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
+  }
 
   delete error;
 }
@@ -1952,13 +2031,74 @@ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t ful
   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
 }
 
+//____________________________________________________________________
+TH1* AliMultiplicityCorrection::GetConvoluted(Int_t i, EventType eventType)
+{
+  // convolutes the corrected histogram i with the response matrix
+  
+  TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected");
+  
+  // undo efficiency correction (hist must not be deleted, is reused)
+  TH1* efficiency = GetEfficiency(i, eventType);
+  //new TCanvas; efficiency->DrawCopy();
+  corrected->Multiply(efficiency);
+  
+  TH2* convoluted = CalculateMultiplicityESD(corrected, i);
+  TH1* convolutedProj = convoluted->ProjectionY("GetConvoluted_convolutedProj", 1, convoluted->GetNbinsX());
+  
+  delete convoluted;
+  delete corrected;
+  
+  return convolutedProj;
+}
+
+//____________________________________________________________________
+TH1* AliMultiplicityCorrection::GetResiduals(Int_t i, EventType eventType, Float_t& residualSum)
+{
+  // creates the residual histogram from the corrected histogram i corresponding to an eventType event sample using the corresponding correlation matrix
+  // residual is : M - UT / eM
+  // residualSum contains the squared sum of the residuals
+  
+  TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected");
+  TH1* convolutedProj = GetConvoluted(i, eventType);
+  
+  Int_t begin = 1;
+  Int_t end = fMultiplicityESD[i]->GetNbinsX();
+  if (fVtxEnd > fVtxBegin)
+  {
+    begin = fVtxBegin;
+    end = fVtxEnd;
+  }
+  TH1* measuredProj = fMultiplicityESD[i]->ProjectionY("measuredProj", begin, end);
+  
+  TH1* residuals = (TH1*) measuredProj->Clone("GetResiduals_residuals");
+  residuals->SetTitle(";measured multiplicity;residuals (M-Ut)/e");
+  residuals->Add(convolutedProj, -1);
+  
+  residualSum = 0;
+  for (Int_t i=1; i<=residuals->GetNbinsX(); i++)
+  {
+    if (measuredProj->GetBinContent(i) > 0)
+      residuals->SetBinContent(i, residuals->GetBinContent(i) / TMath::Sqrt(measuredProj->GetBinContent(i)));
+    residuals->SetBinError(i, 0);
+    
+    if (i > 1)
+      residualSum += residuals->GetBinContent(i) * residuals->GetBinContent(i);
+  }
+  
+  delete corrected;
+  delete convolutedProj;
+  delete measuredProj;
+  
+  return residuals;
+}
+
 //____________________________________________________________________
 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
 {
   // runs the distribution given in inputMC through the response matrix identified by
   // correlationMap and produces a measured distribution
   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
-  // if normalized is set, inputMC is expected to be normalized to the bin width
 
   if (!inputMC)
     return 0;
@@ -1977,6 +2117,9 @@ TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t co
     }
   }
 
+  if (fVtxEnd > fVtxBegin)
+    hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd);
+  
   TH2* corr = (TH2*) hist->Project3D("zy");
   //corr->Rebin2D(2, 1);
   corr->Sumw2();
@@ -2010,6 +2153,7 @@ TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t co
       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
 
       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
+      // TODO fix error
       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
     }
 
index 7ff343e1918ab706a98ff2d539efa5d67718d906..9a376d0cdbcab9c88b5f2e405dacab639844264d 100644 (file)
@@ -43,57 +43,79 @@ class AliMultiplicityCorrection : public TNamed {
     virtual Long64_t Merge(const TCollection* list);
 
     void FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured14);
+    void FillTriggeredEvent(Int_t measured05, Int_t measured10, Int_t measured14);
     void FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, AliPWG0Helper::MCProcessType processType, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll);
 
     void FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured14);
+    void FillNoVertexEvent(Float_t vtx, Bool_t vertexReconstructed, Int_t generated05, Int_t generated10, Int_t generated14, Int_t measured05, Int_t measured10, Int_t measured14);
 
     Bool_t LoadHistograms(const Char_t* dir = 0);
     void SaveHistograms(const char* dir = 0);
     void DrawHistograms();
+    void Rebin2DY(TH2F*& hist, Int_t nBins, Double_t* newBins) const;
+    void Rebin3DY(TH3F*& hist, Int_t nBins, Double_t* newBins) const;
+    void RebinGenerated(Int_t nBins05, Double_t* newBins05, Int_t nBins10, Double_t* newBins10, Int_t nBins13, Double_t* newBins13);
     void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE, EventType eventType = kTrVtx);
 
-    Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* initialConditions = 0);
+    Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t check = kFALSE, TH1* initialConditions = 0, Bool_t errorAsBias = kFALSE);
 
-    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Bool_t determineError = kTRUE);
+    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Int_t determineError = 1);
 
     static TH1* CalculateStdDev(TH1** results, Int_t max);
-    TH1* StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo = 0);
+    TH1* StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo = 0);
 
     Int_t ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
 
     void ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
 
+    void Calculate0Bin(Int_t inputRange, EventType eventType, Int_t zeroBinEvents);
+    Float_t GetFraction0Generated(Int_t inputRange);
+    
     TH2F* GetMultiplicityESD(Int_t i) const { return fMultiplicityESD[i]; }
+    TH1F* GetTriggeredEvents(Int_t i) const { return fTriggeredEvents[i]; }
     TH2F* GetMultiplicityVtx(Int_t i) const { return fMultiplicityVtx[i]; }
     TH2F* GetMultiplicityMB(Int_t i) const { return fMultiplicityMB[i]; }
     TH2F* GetMultiplicityINEL(Int_t i) const { return fMultiplicityINEL[i]; }
     TH2F* GetMultiplicityNSD(Int_t i) const { return fMultiplicityNSD[i]; }
     TH2F* GetMultiplicityMC(Int_t i, EventType eventType) const;
     TH3F* GetCorrelation(Int_t i) const { return fCorrelation[i]; }
+    TH1F* GetNoVertexEvents(Int_t i) const { return fNoVertexEvents[i]; }
     TH1F* GetMultiplicityESDCorrected(Int_t i) const { return fMultiplicityESDCorrected[i]; }
 
     void SetMultiplicityESD(Int_t i, TH2F* const hist)  { fMultiplicityESD[i]  = hist; }
+    void SetTriggeredEvents(Int_t i, TH1F* hist)  { fTriggeredEvents[i]  = hist; }
     void SetMultiplicityVtx(Int_t i, TH2F* const hist)  { fMultiplicityVtx[i]  = hist; }
     void SetMultiplicityMB(Int_t i, TH2F* const hist)   { fMultiplicityMB[i]   = hist; }
     void SetMultiplicityINEL(Int_t i, TH2F* const hist) { fMultiplicityINEL[i] = hist; }
     void SetMultiplicityNSD(Int_t i, TH2F* const hist) { fMultiplicityNSD[i] = hist; }
     void SetMultiplicityMC(Int_t i, EventType eventType, TH2F* const hist);
     void SetCorrelation(Int_t i, TH3F* const hist) { fCorrelation[i] = hist; }
+    void SetNoVertexEvents(Int_t i, TH1F* hist) { fNoVertexEvents[i] = hist; }
     void SetMultiplicityESDCorrected(Int_t i, TH1F* const hist) { fMultiplicityESDCorrected[i] = hist; }
 
     void SetGenMeasFromFunc(const TF1* inputMC, Int_t id);
     TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap);
+    
+    void FixTriggerEfficiencies(Int_t start);
+    
+    TH1* GetConvoluted(Int_t i, EventType eventType);
+    TH1* GetResiduals(Int_t i, EventType eventType, Float_t& residualSum);
 
     void GetComparisonResults(Float_t* const mc = 0, Int_t* const mcLimit = 0, Float_t* const residuals = 0, Float_t* const ratioAverage = 0) const;
 
     TH1* GetEfficiency(Int_t inputRange, EventType eventType);
-    TH1* GetTriggerEfficiency(Int_t inputRange);
+    TH1* GetTriggerEfficiency(Int_t inputRange, EventType eventType);
 
     static void SetQualityRegions(Bool_t SPDStudy);
     Float_t GetQuality(Int_t region) const { return fQuality[region]; }
 
     void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y) const;
+    
+    Float_t GetVertexBegin(Int_t i) { return fgVtxRangeBegin[i]; }
+    Float_t GetVertexEnd(Int_t i) { return fgVtxRangeEnd[i]; }
+    
+    void SetVertexRange(Int_t begin, Int_t end) { fVtxBegin = begin; fVtxEnd = end; }
 
   protected:
     void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
@@ -105,6 +127,8 @@ class AliMultiplicityCorrection : public TNamed {
     TH1* fCurrentEfficiency;  //! current efficiency
 
     TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4 (0..2)
+    TH1F* fTriggeredEvents[kESDHists]; // (raw) multiplicity distribution of triggered events; array: |eta| < 0.5, 1.0, 1.4 (0..2)
+    TH1F* fNoVertexEvents[kESDHists];  // distribution of true multiplicity just of triggered events without vertex or with 0 tracklets; array: |eta| < 0.5, 1.0, 1.4 (0..2)
 
     TH2F* fMultiplicityVtx[kMCHists];  // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4, inf (0..3)
     TH2F* fMultiplicityMB[kMCHists];   // multiplicity histogram of triggered events                        : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4, inf (0..3)
@@ -114,12 +138,18 @@ class AliMultiplicityCorrection : public TNamed {
     TH3F* fCorrelation[kCorrHists];              // vtx vs. (gene multiplicity (trig+vtx)) vs. (meas multiplicity); array: |eta| < 0.5, 1, 1.4, (0..2 and 3..5), the first corrects to the eta range itself, the second to full phase space
 
     TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms
-
+    
     Int_t fLastBinLimit;        //! last bin limit, determined in SetupCurrentHists()
     Float_t fLastChi2MC;        //! last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
     Int_t   fLastChi2MCLimit;   //! bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
     Float_t fLastChi2Residuals; //! last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
     Float_t fRatioAverage;      //! last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
+    
+    Int_t fVtxBegin;            //! vertex range for analysis
+    Int_t fVtxEnd;              //! vertex range for analysis
+    
+    static Double_t fgVtxRangeBegin[kESDHists]; //! begin of allowed vertex range for this eta bin
+    static Double_t fgVtxRangeEnd[kESDHists];   //! end of allowed vertex range for this eta bin
 
     static Int_t   fgQualityRegionsB[kQualityRegions]; //! begin, given in multiplicity units
     static Int_t   fgQualityRegionsE[kQualityRegions]; //! end
@@ -129,7 +159,7 @@ class AliMultiplicityCorrection : public TNamed {
     AliMultiplicityCorrection(const AliMultiplicityCorrection&);
     AliMultiplicityCorrection& operator=(const AliMultiplicityCorrection&);
 
-  ClassDef(AliMultiplicityCorrection, 5);
+  ClassDef(AliMultiplicityCorrection, 7);
 };
 
 #endif
index 6e24595aab31805b3dec48b48139112ff926e4bb..e30aa5b4be493eb9ae1acc5e30414cafa52b6e4c 100644 (file)
@@ -34,6 +34,7 @@
 #include "multiplicity/AliMultiplicityCorrection.h"
 #include "AliCorrection.h"
 #include "AliCorrectionMatrix3D.h"
+#include "AliPhysicsSelection.h"
 #include "AliTriggerAnalysis.h"
 
 ClassImp(AliMultiplicityTask)
@@ -45,6 +46,7 @@ AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn)),
   fTrigger(AliTriggerAnalysis::kMB1),
   fDeltaPhiCut(-1),
+  fDiffTreatment(AliPWG0Helper::kMCFlags),
   fReadMC(kFALSE),
   fUseMCVertex(kFALSE),
   fMultiplicity(0),
@@ -54,6 +56,8 @@ AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
   fParticleSpecies(0),
   fdNdpT(0),
   fPtSpectrum(0),
+  fTemp1(0),
+  fTemp2(0),
   fOutput(0)
 {
   //
@@ -62,10 +66,25 @@ AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
 
   for (Int_t i = 0; i<8; i++)
     fParticleCorrection[i] = 0;
+    
+  for (Int_t i=0; i<3; i++)
+    fEta[i] = 0;
 
   // Define input and output slots here
   DefineInput(0, TChain::Class());
   DefineOutput(0, TList::Class());
+
+  if (fOption.Contains("only-process-type-nd"))
+    fSelectProcessType = 1;
+
+  if (fOption.Contains("only-process-type-sd"))
+    fSelectProcessType = 2;
+
+  if (fOption.Contains("only-process-type-dd"))
+    fSelectProcessType = 3;
+
+  if (fSelectProcessType != 0)
+    AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
 }
 
 AliMultiplicityTask::~AliMultiplicityTask()
@@ -138,12 +157,6 @@ void AliMultiplicityTask::CreateOutputObjects()
   fdNdpT->Sumw2();
   fOutput->Add(fdNdpT);
 
-  if (fOption.Contains("skip-particles"))
-  {
-    fSystSkipParticles = kTRUE;
-    AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
-  }
-
   if (fOption.Contains("particle-efficiency"))
     for (Int_t i = 0; i<8; i++)
     {
@@ -151,18 +164,6 @@ void AliMultiplicityTask::CreateOutputObjects()
       fOutput->Add(fParticleCorrection[i]);
     }
 
-  if (fOption.Contains("only-process-type-nd"))
-    fSelectProcessType = 1;
-
-  if (fOption.Contains("only-process-type-sd"))
-    fSelectProcessType = 2;
-
-  if (fOption.Contains("only-process-type-dd"))
-    fSelectProcessType = 3;
-
-  if (fSelectProcessType != 0)
-    AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
-
   if (fOption.Contains("pt-spectrum-hist"))
   {
     TFile* file = TFile::Open("ptspectrum_fit.root");
@@ -201,6 +202,15 @@ void AliMultiplicityTask::CreateOutputObjects()
     fOutput->Add(fParticleSpecies);
   }
 
+  fTemp1 = new TH2F("fTemp1", "fTemp1", 100, -0.5, 99.5, 100, -0.5, 99.5);
+  fOutput->Add(fTemp1);
+  
+  for (Int_t i=0; i<3; i++)
+  {
+    fEta[i] = new TH1F(Form("fEta_%d", i), ";#eta", 100, -2, 2);
+    fOutput->Add(fEta[i]);
+  }
+  
   // TODO set seed for random generator
 }
 
@@ -215,11 +225,32 @@ void AliMultiplicityTask::Exec(Option_t*)
     return;
   }
 
-  static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
-  Bool_t eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
-  //Printf("%lld", fESD->GetTriggerMask());
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+  if (!inputHandler)
+  {
+    Printf("ERROR: Could not receive input handler");
+    return;
+  }
+    
+  Bool_t eventTriggered = inputHandler->IsEventSelected();
 
+  static AliTriggerAnalysis* triggerAnalysis = 0;
+  if (!triggerAnalysis)
+  {
+    AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+    triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+  }
+  if (eventTriggered)
+    eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
+    
   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+  if (vtxESD && !AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
+    vtxESD = 0;
+    
+  // remove vertices outside +- 15 cm
+  if (vtxESD && TMath::Abs(vtxESD->GetZv()) > 15)
+    vtxESD = 0;
+  
   Bool_t eventVertex = (vtxESD != 0);
 
   Double_t vtx[3];
@@ -237,36 +268,60 @@ void AliMultiplicityTask::Exec(Option_t*)
   Float_t* etaArr = 0;
   if (fAnalysisMode & AliPWG0Helper::kSPD)
   {
-    // get tracklets
-    const AliMultiplicity* mult = fESD->GetMultiplicity();
-    if (!mult)
-    {
-      AliDebug(AliLog::kError, "AliMultiplicity not available");
-      return;
-    }
-
-    labelArr = new Int_t[mult->GetNumberOfTracklets()];
-    etaArr = new Float_t[mult->GetNumberOfTracklets()];
-
-    // get multiplicity from ITS tracklets
-    for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+    if (vtxESD)
     {
-      //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
-
-      Float_t deltaPhi = mult->GetDeltaPhi(i);
-      
-      if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
-        continue;
-
-      etaArr[inputCount] = mult->GetEta(i);
-      if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
+      // get tracklets
+      const AliMultiplicity* mult = fESD->GetMultiplicity();
+      if (!mult)
       {
-        labelArr[inputCount] = mult->GetLabel(i, 0);
+        AliDebug(AliLog::kError, "AliMultiplicity not available");
+        return;
       }
-      else
-        labelArr[inputCount] = -1;
+  
+      labelArr = new Int_t[mult->GetNumberOfTracklets()];
+      etaArr = new Float_t[mult->GetNumberOfTracklets()];
+      
+      Bool_t foundInEta10 = kFALSE;
+      
+      // get multiplicity from ITS tracklets
+      for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+      {
+        //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+  
+        Float_t deltaPhi = mult->GetDeltaPhi(i);
+        
+        if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
+          continue;
+  
+        if (fSystSkipParticles && gRandom->Uniform() < (0.0153))
+        {
+          Printf("Skipped tracklet!");
+          continue;
+        }
+          
+        etaArr[inputCount] = mult->GetEta(i);
+        if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
+        {
+          labelArr[inputCount] = mult->GetLabel(i, 0);
+        }
+        else
+          labelArr[inputCount] = -1;
+          
+        for (Int_t i=0; i<3; i++)
+        {
+          if (vtx[2] > fMultiplicity->GetVertexBegin(i) && vtx[2] < fMultiplicity->GetVertexEnd(i))
+            fEta[i]->Fill(etaArr[inputCount]);
+        }
         
-      ++inputCount;
+        // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
+        if (TMath::Abs(etaArr[inputCount]) < 1)
+          foundInEta10 = kTRUE;
+          
+        ++inputCount;
+      }
+      
+      if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
+        eventTriggered = kFALSE;
     }
   }
   else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
@@ -295,6 +350,16 @@ void AliMultiplicityTask::Exec(Option_t*)
           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
           continue;
         }
+        
+        if (esdTrack->Pt() < 0.15)
+          continue;
+        
+        Float_t d0z0[2],covd0z0[3];
+        esdTrack->GetImpactParameters(d0z0,covd0z0);
+        Float_t sigma= 0.0050+0.0060/TMath::Power(esdTrack->Pt(),0.9);
+        Float_t d0max = 7.*sigma;
+        if (TMath::Abs(d0z0[0]) > d0max) 
+          continue;
   
         etaArr[inputCount] = esdTrack->Eta();
         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
@@ -318,28 +383,44 @@ void AliMultiplicityTask::Exec(Option_t*)
 
   if (!fReadMC) // Processing of ESD information
   {
-    if (eventTriggered && eventVertex)
+    Int_t nESDTracks05 = 0;
+    Int_t nESDTracks10 = 0;
+    Int_t nESDTracks14 = 0;
+    
+    for (Int_t i=0; i<inputCount; ++i)
     {
-      Int_t nESDTracks05 = 0;
-      Int_t nESDTracks10 = 0;
-      Int_t nESDTracks14 = 0;
+      Float_t eta = etaArr[i];
 
-      for (Int_t i=0; i<inputCount; ++i)
-      {
-        Float_t eta = etaArr[i];
+      if (TMath::Abs(eta) < 0.5)
+        nESDTracks05++;
 
-        if (TMath::Abs(eta) < 0.5)
-          nESDTracks05++;
+      if (TMath::Abs(eta) < 1.0)
+        nESDTracks10++;
 
-        if (TMath::Abs(eta) < 1.0)
-          nESDTracks10++;
+      if (TMath::Abs(eta) < 1.3)
+        nESDTracks14++;
+    }
+    
+    // kick out randomly for combinatorics
+    /*
+    if (gRandom->Uniform() < 1.3e-4 * nESDTracks05 * nESDTracks05)
+      nESDTracks05--;
 
-        if (TMath::Abs(eta) < 1.4)
-          nESDTracks14++;
-      }
+    if (gRandom->Uniform() < 8.7e-5 * nESDTracks10 * nESDTracks10)
+      nESDTracks10--;
 
+    if (gRandom->Uniform() < 9.6e-5 * nESDTracks14 * nESDTracks14)
+      nESDTracks14--;
+    */
+
+    //if (nESDTracks05 >= 20 || nESDTracks10 >= 30 || nESDTracks14 >= 32)
+    //  Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d; Tracks: %d %d %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber(), nESDTracks05, nESDTracks10, nESDTracks14);
+
+    if (eventTriggered)
+      fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
+    
+    if (eventTriggered && eventVertex)
       fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
-    }
   }
   else if (fReadMC)   // Processing of MC information
   {
@@ -381,7 +462,7 @@ void AliMultiplicityTask::Exec(Option_t*)
     }
     
     // get process information
-    AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(header);
+    AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
 
     Bool_t processEvent = kTRUE;
     if (fSelectProcessType > 0)
@@ -488,7 +569,7 @@ void AliMultiplicityTask::Exec(Option_t*)
         if (TMath::Abs(particle->Eta()) < 1.0)
           nMCTracks10 += particleWeight;
 
-        if (TMath::Abs(particle->Eta()) < 1.4)
+        if (TMath::Abs(particle->Eta()) < 1.3)
           nMCTracks14 += particleWeight;
 
         nMCTracksAll += particleWeight;
@@ -517,276 +598,284 @@ void AliMultiplicityTask::Exec(Option_t*)
 
       fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
 
-      if (eventTriggered && eventVertex)
-      {
-        Int_t nESDTracks05 = 0;
-        Int_t nESDTracks10 = 0;
-        Int_t nESDTracks14 = 0;
+      // ESD processing
+      Int_t nESDTracks05 = 0;
+      Int_t nESDTracks10 = 0;
+      Int_t nESDTracks14 = 0;
 
-        // tracks per particle species, in |eta| < 2 (systematic study)
-        Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
-        for (Int_t i = 0; i<7; ++i)
-          nESDTracksSpecies[i] = 0;
+      // tracks per particle species, in |eta| < 2 (systematic study)
+      Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
+      for (Int_t i = 0; i<7; ++i)
+        nESDTracksSpecies[i] = 0;
 
-        Bool_t* foundPrimaries = new Bool_t[nPrim];   // to prevent double counting
-        for (Int_t i=0; i<nPrim; i++)
-          foundPrimaries[i] = kFALSE;
+      Bool_t* foundPrimaries = new Bool_t[nPrim];   // to prevent double counting
+      for (Int_t i=0; i<nPrim; i++)
+        foundPrimaries[i] = kFALSE;
 
-        Bool_t* foundPrimaries2 = new Bool_t[nPrim];   // to prevent double counting
-        for (Int_t i=0; i<nPrim; i++)
-          foundPrimaries2[i] = kFALSE;
+      Bool_t* foundPrimaries2 = new Bool_t[nPrim];   // to prevent double counting
+      for (Int_t i=0; i<nPrim; i++)
+        foundPrimaries2[i] = kFALSE;
 
-        Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
-        for (Int_t i=0; i<nMCPart; i++)
-          foundTracks[i] = kFALSE;
+      Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
+      for (Int_t i=0; i<nMCPart; i++)
+        foundTracks[i] = kFALSE;
 
-        for (Int_t i=0; i<inputCount; ++i)
-        {
-          Float_t eta = etaArr[i];
-          Int_t label = labelArr[i];
+      for (Int_t i=0; i<inputCount; ++i)
+      {
+        Float_t eta = etaArr[i];
+        Int_t label = labelArr[i];
 
-          Int_t particleWeight = 1;
+        Int_t particleWeight = 1;
 
-          // systematic study: 5% lower efficiency
-          if (fSystSkipParticles && (gRandom->Uniform() < 0.05))
-            continue;
-      
-          // in case of systematic study, weight according to the change of the pt spectrum
-          if (fPtSpectrum)
-          {
-            TParticle* mother = 0;
+        // in case of systematic study, weight according to the change of the pt spectrum
+        if (fPtSpectrum)
+        {
+          TParticle* mother = 0;
 
-            // preserve label for later
-            Int_t labelCopy = label;
-            if (labelCopy >= 0)
-              labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
-            if (labelCopy >= 0)
-              mother = stack->Particle(labelCopy);
+          // preserve label for later
+          Int_t labelCopy = label;
+          if (labelCopy >= 0)
+            labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
+          if (labelCopy >= 0)
+            mother = stack->Particle(labelCopy);
 
-            // in case of pt study we do not count particles w/o label, because they cannot be scaled
-            if (!mother)
-              continue;
+          // in case of pt study we do not count particles w/o label, because they cannot be scaled
+          if (!mother)
+            continue;
 
-            // it cannot be just multiplied because we cannot count "half of a particle"
-            // instead a random generator decides if the particle is counted twice (if value > 1) 
-            // or not (if value < 0)
-            Int_t bin = fPtSpectrum->FindBin(mother->Pt());
-            if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
+          // it cannot be just multiplied because we cannot count "half of a particle"
+          // instead a random generator decides if the particle is counted twice (if value > 1) 
+          // or not (if value < 0)
+          Int_t bin = fPtSpectrum->FindBin(mother->Pt());
+          if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
+          {
+            Float_t factor = fPtSpectrum->GetBinContent(bin);
+            if (factor > 0)
             {
-              Float_t factor = fPtSpectrum->GetBinContent(bin);
-              if (factor > 0)
+              Float_t random = gRandom->Uniform();
+              if (factor > 1 && random < factor - 1)
               {
-                Float_t random = gRandom->Uniform();
-                if (factor > 1 && random < factor - 1)
-                {
-                  particleWeight = 2;
-                }
-                else if (factor < 1 && random < 1 - factor)
-                  particleWeight = 0;
+                particleWeight = 2;
               }
+              else if (factor < 1 && random < 1 - factor)
+                particleWeight = 0;
             }
           }
+        }
 
-          //Printf("ESD weight is: %d", particleWeight);
+        //Printf("ESD weight is: %d", particleWeight);
 
-          if (TMath::Abs(eta) < 0.5)
-            nESDTracks05 += particleWeight;
+        if (TMath::Abs(eta) < 0.5)
+          nESDTracks05 += particleWeight;
 
-          if (TMath::Abs(eta) < 1.0)
-            nESDTracks10 += particleWeight;
+        if (TMath::Abs(eta) < 1.0)
+          nESDTracks10 += particleWeight;
 
-          if (TMath::Abs(eta) < 1.4)
-            nESDTracks14 += particleWeight;
+        if (TMath::Abs(eta) < 1.3)
+          nESDTracks14 += particleWeight;
 
-          if (fParticleSpecies)
+        if (fParticleSpecies)
+        {
+          Int_t motherLabel = -1;
+          TParticle* mother = 0;
+
+          // find mother
+          if (label >= 0)
+            motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
+          if (motherLabel >= 0)
+            mother = stack->Particle(motherLabel);
+
+          if (!mother)
+          {
+            // count tracks that did not have a label
+            if (TMath::Abs(eta) < etaRange)
+              nESDTracksSpecies[4]++;
+          }
+          else
           {
-            Int_t motherLabel = -1;
-            TParticle* mother = 0;
+            // get particle type (pion, proton, kaon, other) of mother
+            Int_t idMother = -1;
+            switch (TMath::Abs(mother->GetPdgCode()))
+            {
+              case 211: idMother = 0; break;
+              case 321: idMother = 1; break;
+              case 2212: idMother = 2; break;
+              default: idMother = 3; break;
+            }
+
+            // double counting is ok for particle ratio study
+            if (TMath::Abs(eta) < etaRange)
+              nESDTracksSpecies[idMother]++;
 
-            // find mother
-            if (label >= 0)
-              motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
-            if (motherLabel >= 0)
-              mother = stack->Particle(motherLabel);
+            // double counting is not ok for efficiency study
 
-            if (!mother)
+            // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
+            if (foundTracks[label])
             {
-              // count tracks that did not have a label
               if (TMath::Abs(eta) < etaRange)
-                nESDTracksSpecies[4]++;
+                nESDTracksSpecies[6]++;
             }
             else
             {
-              // get particle type (pion, proton, kaon, other) of mother
-              Int_t idMother = -1;
-              switch (TMath::Abs(mother->GetPdgCode()))
-              {
-                case 211: idMother = 0; break;
-                case 321: idMother = 1; break;
-                case 2212: idMother = 2; break;
-                default: idMother = 3; break;
-              }
-
-              // double counting is ok for particle ratio study
-              if (TMath::Abs(eta) < etaRange)
-                nESDTracksSpecies[idMother]++;
+              foundTracks[label] = kTRUE;
 
-              // double counting is not ok for efficiency study
-
-              // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
-              if (foundTracks[label])
+              // particle (primary) already counted?
+              if (foundPrimaries[motherLabel])
               {
                 if (TMath::Abs(eta) < etaRange)
-                  nESDTracksSpecies[6]++;
+                  nESDTracksSpecies[5]++;
               }
               else
-              {
-                foundTracks[label] = kTRUE;
-
-                // particle (primary) already counted?
-                if (foundPrimaries[motherLabel])
-                {
-                  if (TMath::Abs(eta) < etaRange)
-                    nESDTracksSpecies[5]++;
-                }
-                else
-                  foundPrimaries[motherLabel] = kTRUE;
-              }
+                foundPrimaries[motherLabel] = kTRUE;
             }
           }
+        }
 
-          if (fParticleCorrection[0])
+        if (fParticleCorrection[0])
+        {
+          if (label >= 0 && stack->IsPhysicalPrimary(label))
           {
-            if (label >= 0 && stack->IsPhysicalPrimary(label))
-            {
-              TParticle* particle = stack->Particle(label);
+            TParticle* particle = stack->Particle(label);
 
-              // get particle type (pion, proton, kaon, other)
-              Int_t id = -1;
-              switch (TMath::Abs(particle->GetPdgCode()))
-              {
-                case 211: id = 0; break;
-                case 321: id = 1; break;
-                case 2212: id = 2; break;
-                default: id = 3; break;
-              }
+            // get particle type (pion, proton, kaon, other)
+            Int_t id = -1;
+            switch (TMath::Abs(particle->GetPdgCode()))
+            {
+              case 211: id = 0; break;
+              case 321: id = 1; break;
+              case 2212: id = 2; break;
+              default: id = 3; break;
+            }
 
-              // todo check if values are not completely off??
+            // todo check if values are not completely off??
 
-              // particle (primary) already counted?
-              if (!foundPrimaries2[label])
-              {
-                foundPrimaries2[label] = kTRUE;
-                fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
-              }
+            // particle (primary) already counted?
+            if (!foundPrimaries2[label])
+            {
+              foundPrimaries2[label] = kTRUE;
+              fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
             }
           }
         }
-          
-        if (fParticleCorrection[0])
+      }
+        
+      if (fParticleCorrection[0])
+      {
+        // if the particle decays/stops before this radius we do not see it
+        // 8cm larger than SPD layer 2
+        // 123cm TPC radius where a track has about 50 clusters (cut limit)          
+        const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
+                
+        // loop over all primaries that have not been found
+        for (Int_t i=0; i<nPrim; i++)
         {
-          // if the particle decays/stops before this radius we do not see it
-          // 8cm larger than SPD layer 2
-          // 123cm TPC radius where a track has about 50 clusters (cut limit)          
-          const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
-                  
-          // loop over all primaries that have not been found
-          for (Int_t i=0; i<nPrim; i++)
-          {
-            // already found
-            if (foundPrimaries2[i])
-              continue;
-              
-            TParticle* particle = 0;
-            TClonesArray* trackrefs = 0;
-            mcEvent->GetParticleAndTR(i, particle, trackrefs);
-            
-            // true primary and charged
-            if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
-              continue;              
-            
-            //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
-            if (TMath::Abs(particle->Eta()) > 3)
-              continue;
-            
-            // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
+          // already found
+          if (foundPrimaries2[i])
+            continue;
             
-            // get particle type (pion, proton, kaon, other)
-            Int_t id = -1;
-            switch (TMath::Abs(particle->GetPdgCode()))
-            {
-              case 211: id = 4; break;
-              case 321: id = 5; break;
-              case 2212: id = 6; break;
-              default: id = 7; break;
-            }            
+          TParticle* particle = 0;
+          TClonesArray* trackrefs = 0;
+          mcEvent->GetParticleAndTR(i, particle, trackrefs);
+          
+          // true primary and charged
+          if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
+            continue;              
+          
+          //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
+          if (TMath::Abs(particle->Eta()) > 3)
+            continue;
+          
+          // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
+          
+          // get particle type (pion, proton, kaon, other)
+          Int_t id = -1;
+          switch (TMath::Abs(particle->GetPdgCode()))
+          {
+            case 211: id = 4; break;
+            case 321: id = 5; break;
+            case 2212: id = 6; break;
+            default: id = 7; break;
+          }            
+          
+          if (!fParticleCorrection[id])
+            continue;
             
-            if (!fParticleCorrection[id])
-              continue;
-              
-            // get last track reference
-            AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
+          // get last track reference
+          AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
+          
+          if (!trackref)
+          {
+            Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
+            particle->Print();
+            continue;
+          }
             
-            if (!trackref)
+          // particle in tracking volume long enough...
+          if (trackref->R() > endRadius)
+            continue;  
+          
+          if (particle->GetLastDaughter() >= 0)
+          {
+            Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
+            //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
+            if (uID == kPDecay)
             {
-              Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
-              particle->Print();
-              continue;
-            }
+              // decayed
               
-            // particle in tracking volume long enough...
-            if (trackref->R() > endRadius)
-              continue;  
-            
-            if (particle->GetLastDaughter() >= 0)
-            {
-              Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
-              //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
-              if (uID == kPDecay)
-              {
-                // decayed
-                
-                Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
-                particle->Print();
-                Printf("Daughers:");
-                for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
-                  stack->Particle(d)->Print();
-                Printf("");
-                
-                fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
-                continue;
-              }
-            }
-            
-            if (trackref->DetectorId() == -1)
-            {
-              // stopped
-              Printf("Particle %d stopped at %f:", i, trackref->R());
+              Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
               particle->Print();
+              Printf("Daughers:");
+              for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
+                stack->Particle(d)->Print();
               Printf("");
               
-              fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
+              fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
               continue;
             }
-            
-            Printf("Particle %d simply not tracked", i);
+          }
+          
+          if (trackref->DetectorId() == -1)
+          {
+            // stopped
+            Printf("Particle %d stopped at %f:", i, trackref->R());
             particle->Print();
             Printf("");
+            
+            fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
+            continue;
           }
+          
+          Printf("Particle %d simply not tracked", i);
+          particle->Print();
+          Printf("");
         }
-        
-        delete[] foundTracks;
-        delete[] foundPrimaries;
-        delete[] foundPrimaries2;
+      }
+      
+      delete[] foundTracks;
+      delete[] foundPrimaries;
+      delete[] foundPrimaries2;
 
-        if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
-        {
-            TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
-            printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
-        }
+//         if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
+//         {
+//             TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+//             printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
+//         }
 
+      if (eventTriggered)
+      {
+        fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
+        fMultiplicity->FillNoVertexEvent(vtxMC[2], eventVertex, nMCTracks05, nMCTracks10, nMCTracks14, nESDTracks05, nESDTracks10, nESDTracks14);
+//         if (!eventVertex)
+//         {
+//           if (nESDTracks05 == 0)
+//             fTemp1->Fill(nMCTracks05, nESDTracks05);
+//         }
+      }
+      
+      if (eventTriggered && eventVertex)
+      {
         // fill response matrix using vtxMC (best guess)
-        fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks14,  nMCTracksAll,  nESDTracks05,  nESDTracks10, nESDTracks14);
+        fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks14,  nMCTracksAll,  nESDTracks05, nESDTracks10, nESDTracks14);
 
         fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 
@@ -827,7 +916,15 @@ void AliMultiplicityTask::Terminate(Option_t *)
     return;
   }
 
-  TFile* file = TFile::Open("multiplicity.root", "RECREATE");
+  TString fileName("multiplicity");
+  if (fSelectProcessType == 1)
+    fileName += "ND";
+  if (fSelectProcessType == 2)
+    fileName += "SD";
+  if (fSelectProcessType == 3)
+    fileName += "DD";
+  fileName += ".root";
+  TFile* file = TFile::Open(fileName, "RECREATE");
 
   fMultiplicity->SaveHistograms();
   for (Int_t i = 0; i < 8; ++i)
@@ -838,10 +935,28 @@ void AliMultiplicityTask::Terminate(Option_t *)
   if (fdNdpT)
     fdNdpT->Write();
 
+  fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
+  if (fTemp1)
+    fTemp1->Write();
+    
+  for (Int_t i=0; i<3; i++)
+  {
+    fEta[i] = dynamic_cast<TH1*> (fOutput->FindObject(Form("fEta_%d", i)));
+    if (fEta[i])
+    {
+      fEta[i]->Sumw2();
+      Float_t events = fMultiplicity->GetMultiplicityESD(i)->Integral(1, fMultiplicity->GetMultiplicityESD(i)->GetNbinsX());
+      if (events > 0)
+        fEta[i]->Scale(1.0 / events);
+      fEta[i]->Scale(1.0 / fEta[i]->GetXaxis()->GetBinWidth(1));
+      fEta[i]->Write();
+    }
+  }
+  
   TObjString option(fOption);
   option.Write();
 
   file->Close();
 
-  Printf("Written result to multiplicity.root");
+  Printf("Written result to %s", fileName.Data());
 }
index b450ec446aa90532026915a9ab68b160dbb52f0b..e4ec4e1df7224ffe9fc881452bc988d1c6e471cd 100644 (file)
@@ -35,7 +35,9 @@ class AliMultiplicityTask : public AliAnalysisTask {
 
     void SetReadMC(Bool_t flag = kTRUE) { fReadMC = flag; }
     void SetUseMCVertex(Bool_t flag = kTRUE) { fUseMCVertex = flag; }
-
+    void SetSkipParticles(Bool_t flag = kTRUE) { fSystSkipParticles = flag; }
+    void SetDiffTreatment(AliPWG0Helper::DiffTreatment diffTreatment) { fDiffTreatment = diffTreatment; }
+  
  protected:
     AliESDEvent *fESD;    //! ESD object
 
@@ -43,24 +45,30 @@ class AliMultiplicityTask : public AliAnalysisTask {
     AliPWG0Helper::AnalysisMode fAnalysisMode; // detector that is used for analysis
     AliTriggerAnalysis::Trigger fTrigger;      // trigger that is used
     Float_t fDeltaPhiCut;                      // cut in delta phi (only SPD)
-
+    AliPWG0Helper::DiffTreatment  fDiffTreatment;  // how to identify SD events (see AliPWG0Helper::GetEventProcessType)
+    
     Bool_t  fReadMC;       // if true reads MC data (to build correlation maps)
     Bool_t  fUseMCVertex;  // the MC vtx is used instead of the ESD vertex (for syst. check)
 
     AliMultiplicityCorrection* fMultiplicity; //! object containing the extracted data
     AliESDtrackCuts* fEsdTrackCuts;           // Object containing the parameters of the esd track cuts
 
-    Bool_t fSystSkipParticles;     //! if true skips particles (systematic study)
+    Bool_t fSystSkipParticles;          // if true skips particles (systematic study)
     AliCorrection* fParticleCorrection[8]; //! correction from measured to generated particles for different particles for trigger, vertex sample in |eta| < 2; switch on with "particle-efficiency"
                                            // for each of the species (0..3): pi, k, p, other; for systematic study of pt cut off
                                            // 4..7 counts for the same species the decayed particles (in generated) and stopped (in measured)
-    Int_t fSelectProcessType;        //! 0 = all (default), 1 = ND, 2 = SD, 3 = DD (for systematic study)
+    Int_t fSelectProcessType;        // 0 = all (default), 1 = ND, 2 = SD, 3 = DD (for systematic study)
     TNtuple *fParticleSpecies;       //! per event: vtx_mc, (pi, k, p, rest (in |eta| < 1)) X (true, recon) + (nolabel,
                                      // doubleTracks, doublePrimaries) [doubleTracks + doublePrimaries are already part of
                                      // rec. particles!); enable with: particle-species
     TH1* fdNdpT;                     //! true pT spectrum (MC)
 
     TH1D* fPtSpectrum;               // function that modifies the pt spectrum (syst. study)
+    
+    TH1* fTemp1;                                 //! temp histogram for quick study of variables
+    TH1* fTemp2;                                 //! temp histogram for quick study of variables
+    
+    TH1* fEta[3];                    //! eta histogram of events in the acceptance region for each of the eta-bins (control histogram)
 
     TList* fOutput;                  //! list send on output slot 0
     
index 7bff843a0faff6c97a693857f8923721d2136042..97809d010a0a5f45f539b4afcd01719f069e3910 100644 (file)
@@ -4,6 +4,21 @@
 // script to correct the multiplicity spectrum + helpers
 //
 
+Bool_t is900GeV = 0;
+Bool_t is2360TeV = 0;
+
+// 900 GeV, MC
+//const Int_t kBinLimits[]   = { 42, 57, 60 };
+//const Int_t kTrustLimits[] = { 26, 42, 54 };
+
+// 2.36 TeV
+//const Int_t kBinLimits[]   = { 43, 67, 83 };
+//const Int_t kTrustLimits[] = { 33, 57, 73 };
+
+// 7 TeV
+const Int_t kBinLimits[]   = { 60, 120, 60 };
+const Int_t kTrustLimits[] = { 26, 70, 54 };
+
 void SetTPC()
 {
   gSystem->Load("libPWG0base");
@@ -60,32 +75,134 @@ void loadlibs()
   gSystem->Load("libPWG0base");
 }
 
-void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = kFALSE, Int_t histID = 1, Bool_t fullPhaseSpace = kFALSE, Float_t beta  = 1e5, Int_t eventType = 0 /* AliMultiplicityCorrection::kTrVtx */, const char* targetFile = "unfolded.root", const char* folderESD = "Multiplicity")
+void LoadAndInitialize(void* multVoid, void* esdVoid, void* multTriggerVoid, Int_t histID, Bool_t fullPhaseSpace, Int_t* geneLimits)
 {
-  loadlibs();
-
-  AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
-  AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
+  AliMultiplicityCorrection* mult = (AliMultiplicityCorrection*) multVoid;
+  AliMultiplicityCorrection* esd = (AliMultiplicityCorrection*) esdVoid;
+  AliMultiplicityCorrection* multTrigger = (AliMultiplicityCorrection*) multTriggerVoid;
   
-  //AliUnfolding::SetNbins(150, 150);
-  //AliUnfolding::SetNbins(65, 65); 
-  //AliUnfolding::SetNbins(35, 35);
-  //AliUnfolding::SetDebug(1);
-
-  for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
+  for (Int_t i=0; i<3; i++)
+    geneLimits[i] = kBinLimits[i];
+  
+  // REBINNING
+  if (1)
   {
-    switch (hID)
+    // 900 GeV
+    if (is900GeV)
+    {
+      if (1)
+      {
+        Int_t bins05 = 31;
+        Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+                  10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+                  20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5, 
+                  100.5 };
+      }
+
+      if (0)
+      {
+        Int_t bins05 = 29;
+        Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+                10.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+                20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5, 
+                100.5 };
+      }
+      
+      if (0)
+      {
+        Int_t bins05 = 25;
+        Double_t* newBinsEta05 = new Double_t[bins05+1];
+      
+        //newBinsEta05[0] = -0.5;
+        //newBinsEta05[1] = 0.5;
+        //newBinsEta05[2] = 1.5;
+        //newBinsEta05[3] = 2.5;
+        
+        for (Int_t i=0; i<bins05+1; i++)
+          newBinsEta05[i] = -0.5 + i*2;
+          //newBinsEta05[i] = 4.5 + (i-4)*2;
+      }
+
+      Int_t bins10 = 54;
+      Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+                                  10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+                                  20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
+                                  30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
+                                  40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5, 
+                                  60.5, 65.5, 70.5, 100.5 };
+                        
+      geneLimits[0] = bins05;
+      geneLimits[1] = bins10;
+      geneLimits[2] = bins10;
+    }
+    
+    // 2.36 TeV
+    if (is2360TeV)
     {
-      case 0: AliUnfolding::SetNbins(35, 35); break;
-      case 1: AliUnfolding::SetNbins(60, 60); break;
-      case 2: AliUnfolding::SetNbins(70, 70); beta *= 3; break;
+      Int_t bins05 = 36;
+      Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+                                 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+                                 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5, 
+                                 40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
+      Int_t bins10 = 64;
+      Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+                                 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+                                 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
+                                 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
+                                 40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5, 
+                                 60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5, 
+                                 80.5, 85.5, 90.5, 100.5 };
+
+      geneLimits[0] = bins05;
+      geneLimits[1] = bins10;
+      geneLimits[2] = bins10;
     }
+    
+    // 7 TeV
+    if (!is900GeV && !is2360TeV)
+    {
+      Int_t bins05 = 36;
+      Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+                                 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+                                 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5, 
+                                 40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
+      
+      Int_t bins10 = 85;
+      Double_t newBinsEta10[] = { 
+          -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
+                                 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
+                                 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
+                                 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
+                                 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 
+                                 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 
+                                 60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5, 
+                                 80.5, 82.5, 84.5, 86.5, 88.5, 90.5, 92.5, 94.5, 96.5, 98.5, 
+                                 100.5, 105.5, 110.5, 115.5, 120.5 };
+      
+      geneLimits[0] = bins05;
+      geneLimits[1] = bins10;
+      geneLimits[2] = bins10;
+    }
+
+    esd->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
+    mult->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
+    multTrigger->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
+  }
   
+  for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
+  {
     TH2F* hist = esd->GetMultiplicityESD(hID);
-    TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
   
     mult->SetMultiplicityESD(hID, hist);
-  
+    mult->SetTriggeredEvents(hID, esd->GetTriggeredEvents(hID));
+    
+    // insert trigger efficiency in flat response matrix
+    for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
+      mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
+    
+    mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
+    mult->FixTriggerEfficiencies(10);
+      
     // small hack to get around charge conservation for full phase space ;-)
     if (fullPhaseSpace)
     {
@@ -98,74 +215,371 @@ void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder
           corr->SetBinError(i, j, corr->GetBinError(i-1, j));
         }
     }
+  }
+}
+
+void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = 1, Int_t histID = 1, Bool_t fullPhaseSpace = kFALSE, Float_t beta = -25, Int_t eventType = 2 /* AliMultiplicityCorrection::kTrVtx */, const char* targetFile = "unfolded.root", const char* folderESD = "Multiplicity", Bool_t calcBias = kFALSE)
+{
+  loadlibs();
+  
+  Int_t geneLimits[] = { 0, 0, 0 };
   
+  AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
+  AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
+  AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
+
+  LoadAndInitialize(mult, esd, multTrigger, histID, fullPhaseSpace, geneLimits);
+  
+  //AliUnfolding::SetDebug(1);
+
+  for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
+  {
+    AliUnfolding::SetNbins(kBinLimits[hID], geneLimits[hID]);
+  
+    TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
+
     if (chi2)
     {
-      AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, beta);
+      if (is900GeV)
+      {
+        //Float_t regParamPol0[] = { 2, 2, 10 };   // TPCITS
+        Float_t regParamPol0[] = { 5, 15, 25 }; // SPD
+        Float_t regParamPol1[] =  { 0.15, 0.25, 0.5 };
+      }
+      else if (is2360TeV)
+      {
+        Float_t regParamPol0[] = { 10, 12, 40 };
+        Float_t regParamPol1[] =  { 0.25, 0.25, 2 };
+      }
+      else
+      {
+        Float_t regParamPol0[] = { 1, 25, 10 };
+        Float_t regParamPol1[] =  { 0.15, 0.5, 0.5 };
+        AliUnfolding::SetSkipBinsBegin(3);
+      }
+      
+      Int_t reg = AliUnfolding::kPol0;
+      if (beta > 0)
+        reg = AliUnfolding::kPol1;
+        
+      Float_t regParam = TMath::Abs(beta);
+      if (histID == -1)
+      {
+        if (beta > 0)
+          regParam = regParamPol1[hID];
+        else
+          regParam = regParamPol0[hID];
+      }
+      AliUnfolding::SetChi2Regularization(reg, regParam);
+      
+      //AliUnfolding::SetChi2Regularization(AliUnfolding::kLog, 1000000);
+      //AliUnfolding::SetChi2Regularization(AliUnfolding::kRatio, 10);
+      //TVirtualFitter::SetDefaultFitter("Minuit2");
+      
       //AliUnfolding::SetCreateOverflowBin(kFALSE);
       //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
       //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
      // AliUnfolding::SetChi2Regularization(AliUnfolding::kEntropy, beta);
       //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
       
-      //mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kTRUE, mcCompare);
-      //mult->SetMultiplicityESDCorrected(histID, (TH1F*) mcCompare);
+      if (0)
+      {
+        // part for checking
+        TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
+        mcCompare->Sumw2();
+        mult->ApplyMinuitFit(hID, fullPhaseSpace, AliMultiplicityCorrection::kMB, 0, kTRUE, mcCompare);
+        mult->SetMultiplicityESDCorrected(hID, (TH1F*) mcCompare);
+      }
+      else
+      {
+        // Zero Bin
+        Int_t zeroBin = 0;
+        if (is900GeV) // from data
+        {
+          // background subtraction
+          Int_t background = 0;
+          
+          //background = 1091 + 4398; // MB1 for run 104824 - 52 (SPD)
+          //background = 1087 + 4398; // MB1 for run 104824 - 52 (TPCITS)
+          
+          //background = 417 + 1758; // MB1 for run 104867 - 92 (SPD)
+          //background = 1162+422;     // MB1 for run 104892 (TPCITS)
+          //background = 5830 + 1384;    // MB1 for run 104824,25,45,52,67,92 (TPC runs) (TPCITS)
+          
+          background = 20 + 0; // V0AND for run 104824 - 52
+          //background = 10 + 0; // V0AND for run 104867 - 92
+          
+          Printf("NOTE: Subtracting %d background events", background);
+          gSystem->Sleep(1000);
+      
+          zeroBin = mult->GetTriggeredEvents(hID)->GetBinContent(1) - background - mult->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1);
+        }
+        else if (is2360TeV)
+        {
+          // from MC
+          Float_t fractionZeroBin = (multTrigger->GetTriggeredEvents(hID)->GetBinContent(1) - multTrigger->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1)) / multTrigger->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
+          zeroBin = fractionZeroBin * mult->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
+          
+          Printf("Zero bin from MC: Estimating %d events with trigger but without vertex", zeroBin);
+          gSystem->Sleep(1000);
+        }
+        else
+        {
+          AliUnfolding::SetSkip0BinInChi2(kTRUE);
+        }
       
-      mult->ApplyMinuitFit(hID, fullPhaseSpace, eventType, kFALSE, 0, kFALSE); //hist2->ProjectionY("mymchist"));
-      //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType);
+        //mult->SetVertexRange(3, 4);
+        mult->ApplyMinuitFit(hID, fullPhaseSpace, eventType, zeroBin, kFALSE, 0, calcBias);
+      }
     }
     else
     {
-      mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, 10, 0, kTRUE);
+      // HACK
+      //mult->GetMultiplicityESD(hID)->SetBinContent(1, 1, 0);
+      //for (Int_t bin=1; bin<=mult->GetCorrelation(hID)->GetNbinsY(); bin++)
+      //  mult->GetCorrelation(hID)->SetBinContent(1, bin, 1, 0);
+      AliUnfolding::SetChi2Regularization(AliUnfolding::kNone, 0);
+      mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, beta, 0, kTRUE);
       //mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, 5, 0, kFALSE);
     }
-  
-    mult->SetMultiplicityMC(hID, eventType, hist2);
   }
   
   Printf("Writing result to %s", targetFile);
   TFile* file = TFile::Open(targetFile, "RECREATE");
-  mult->SaveHistograms();
+  mult->SaveHistograms("Multiplicity");
   file->Write();
   file->Close();
 
+  if (histID == -1)
+    return;
+
   for (hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
   {
     TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
     TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
-    mult->DrawComparison(Form("%s_%d", (chi2) ? "MinuitChi2" : "Bayesian", hID), hID, fullPhaseSpace, kTRUE, mcCompare);
+    mult->DrawComparison(Form("%s_%d_%f", (chi2) ? "MinuitChi2" : "Bayesian", hID, beta), hID, fullPhaseSpace, kTRUE, mcCompare, kFALSE, eventType);
+
+    Printf("<n> = %f", mult->GetMultiplicityESDCorrected(hID)->GetMean());
   }
 }
 
-TH1* GetChi2Bias(Float_t alpha)
+void correctAll(Int_t eventType)
+{
+  const char* suffix = "";
+  switch (eventType)
+  {
+    case 0: suffix = "trvtx"; break;
+    case 1: suffix = "mb"; break;
+    case 2: suffix = "inel"; break;
+    case 3: suffix = "nsd"; break;
+  }
+
+  correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0, -1, eventType, TString(Form("chi2a_%s.root", suffix)));
+  correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0,  1, eventType, TString(Form("chi2b_%s.root", suffix)));
+  correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 0, -1, 0, 40, eventType, TString(Form("bayes_%s.root", suffix)));
+  
+  if (eventType == 3)
+    drawAll(1);
+  else if (eventType == 2)
+    drawAllINEL();
+}
+
+void drawAll(Bool_t showUA5 = kFALSE)
+{
+  const char* files[] = { "chi2a_nsd.root", "chi2b_nsd.root", "bayes_nsd.root" };
+  drawAll(files, showUA5);
+}
+
+void drawAllINEL()
+{
+  const char* files[] = { "chi2a_inel.root", "chi2b_inel.root", "bayes_inel.root" };
+  drawAll(files);
+}
+
+void drawAllMB()
+{
+  const char* files[] = { "chi2a_mb.root", "chi2b_mb.root", "bayes_mb.root" };
+  drawAll(files);
+}
+
+void drawAll(const char** files, Bool_t showUA5 = kFALSE, Bool_t normalize = kFALSE)
 {
   loadlibs();
   
-  const char* fileNameMC = "multiplicityMC.root";
-  const char* folder = "Multiplicity";
-  const char* fileNameESD = "multiplicityESD.root";
-  const char* folderESD = "Multiplicity";
+  Int_t colors[] = { 1, 2, 4 };
+  
+  c = new TCanvas(Form("c%d", gRandom->Uniform(100)), "c", 1800, 600);
+  c->Divide(3, 1);
+  
+  l = new TLegend(0.6, 0.6, 0.99, 0.99);
+  l->SetFillColor(0);
+  
+  TH1* hist0 = 0;
+  
+  for (Int_t hID=0; hID<3; hID++)
+  {
+    c->cd(hID+1)->SetLogy();
+    gPad->SetGridx();
+    gPad->SetGridy();
+    for (Int_t i=0; i<3; i++)
+    {
+      TFile::Open(files[i]);
+      
+      hist = (TH1*) gFile->Get(Form("Multiplicity/fMultiplicityESDCorrected%d", hID));
 
-  AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
-  AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
+      Float_t average0 = hist->GetMean();
+      hist2 = (TH1*) hist->Clone("temp");
+      hist2->SetBinContent(1, 0);
+      Float_t average1 = hist2->GetMean();
+      Printf("%d: <N> = %.2f <N>/(eta) = %.2f | without 0: <N> = %.2f <N>/(eta) = %.2f", hID, average0, average0 / ((hID+1) - 0.4 * (hID / 2)), average1, average1 / ((hID+1) - 0.4 * (hID / 2)));
+      
+      hist->SetLineColor(colors[i]);
+      
+      hist->SetStats(0);
+      hist->GetXaxis()->SetRangeUser(0, kBinLimits[hID]);
+      
+      Float_t min = 0.1;
+      Float_t max = hist->GetMaximum() * 1.5;
+      
+      if (normalize)
+        min = 1e-6;
+      
+      hist->GetYaxis()->SetRangeUser(min, max);
+      hist->SetTitle(";unfolded multiplicity;events");
+      
+      if (hID == 0)
+      {
+        l->AddEntry(hist, files[i]);
+      }
+      
+      if (hist->Integral() <= 0)
+        continue;
+      
+      if (normalize)
+        hist->Scale(1.0 / hist->Integral());
+      
+      AliPWG0Helper::NormalizeToBinWidth(hist);
+      
+      hist->DrawCopy((i == 0) ? "" : "SAME");
+      
+      if (!hist0)
+        hist0 = (TH1*) hist->Clone("hist0");
+    }
+      
+    if (hist0)
+    {
+      line = new TLine(kTrustLimits[hID], hist0->GetMinimum(), kTrustLimits[hID], hist0->GetMaximum());
+      line->SetLineWidth(3);
+      line->Draw();
+    }
+  }
   
-  //Int_t hID = 0; const Int_t kMaxBins = 35;
-  //Int_t hID = 1;  const Int_t kMaxBins = 60;
-  Int_t hID = 2;  const Int_t kMaxBins = 70;
-  AliUnfolding::SetNbins(kMaxBins, kMaxBins);
-  //AliUnfolding::SetDebug(1);
-  //AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 50);
-  AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, alpha);
+  c->cd(1);
+  l->Draw();
   
-  TH2F* hist = esd->GetMultiplicityESD(hID);
-  mult->SetMultiplicityESD(hID, hist);  
+  if (showUA5)
+  {
+    TGraphErrors *gre = new TGraphErrors(23);
+    gre->SetName("Graph");
+    gre->SetTitle("UA5");
+    gre->SetFillColor(1);
+    gre->SetMarkerStyle(24);
+    gre->SetPoint(0,0,0.15);
+    gre->SetPointError(0,0.5,0.01);
+    gre->SetPoint(1,1,0.171);
+    gre->SetPointError(1,0.5,0.007);
+    gre->SetPoint(2,2,0.153);
+    gre->SetPointError(2,0.5,0.007);
+    gre->SetPoint(3,3,0.124);
+    gre->SetPointError(3,0.5,0.006);
+    gre->SetPoint(4,4,0.099);
+    gre->SetPointError(4,0.5,0.005);
+    gre->SetPoint(5,5,0.076);
+    gre->SetPointError(5,0.5,0.005);
+    gre->SetPoint(6,6,0.057);
+    gre->SetPointError(6,0.5,0.004);
+    gre->SetPoint(7,7,0.043);
+    gre->SetPointError(7,0.5,0.004);
+    gre->SetPoint(8,8,0.032);
+    gre->SetPointError(8,0.5,0.003);
+    gre->SetPoint(9,9,0.024);
+    gre->SetPointError(9,0.5,0.003);
+    gre->SetPoint(10,10,0.018);
+    gre->SetPointError(10,0.5,0.002);
+    gre->SetPoint(11,11,0.013);
+    gre->SetPointError(11,0.5,0.002);
+    gre->SetPoint(12,12,0.01);
+    gre->SetPointError(12,0.5,0.002);
+    gre->SetPoint(13,13,0.007);
+    gre->SetPointError(13,0.5,0.001);
+    gre->SetPoint(14,14,0.005);
+    gre->SetPointError(14,0.5,0.001);
+    gre->SetPoint(15,15,0.004);
+    gre->SetPointError(15,0.5,0.001);
+    gre->SetPoint(16,16,0.0027);
+    gre->SetPointError(16,0.5,0.0009);
+    gre->SetPoint(17,17,0.0021);
+    gre->SetPointError(17,0.5,0.0008);
+    gre->SetPoint(18,18,0.0015);
+    gre->SetPointError(18,0.5,0.0007);
+    gre->SetPoint(19,19,0.0011);
+    gre->SetPointError(19,0.5,0.0006);
+    gre->SetPoint(20,20,0.0008);
+    gre->SetPointError(20,0.5,0.0005);
+    gre->SetPoint(21,22,0.00065);
+    gre->SetPointError(21,1,0.0003);
+    gre->SetPoint(22,27.5,0.00013);
+    gre->SetPointError(22,3.5,7.1e-05);
+    
+    for (Int_t i=0; i<gre->GetN(); i++)
+    {
+      gre->GetY()[i] *= hist0->Integral();
+      gre->GetEY()[i] *= hist0->Integral();
+    }
+    
+    gre->Draw("P");
+    
+    ratio = (TGraphErrors*) gre->Clone("ratio");
+    
+    for (Int_t i=0; i<gre->GetN(); i++)
+    {
+      //Printf("%f %d", hist0->GetBinContent(hist0->FindBin(ratio->GetX()[i])), hist0->FindBin(ratio->GetX()[i]));
+      Int_t bin = hist0->FindBin(gre->GetX()[i]);
+      if (hist0->GetBinContent(bin) > 0)
+      {
+        ratio->GetEY()[i] = TMath::Sqrt((ratio->GetEY()[i] / ratio->GetY()[i])**2 + (hist0->GetBinError(bin) / hist0->GetBinContent(bin))**2);
+        ratio->GetY()[i] /= hist0->GetBinContent(bin);
+        ratio->GetEY()[i] *= ratio->GetY()[i];
+      }
+    }
+    new TCanvas;
+    gPad->SetGridx();
+    gPad->SetGridy();
+    //ratio->Print();
+    ratio->Draw("AP");
+    ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
+  }
+  
+  c->SaveAs("draw_all.png");
+}
+
+TH1* GetChi2Bias(Float_t alpha = -5)
+{
+  loadlibs();
+
+  AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kMB;
+  Int_t hID = 1;
   
+  correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "nobias.root", "Multiplicity", kFALSE);
+
+  correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "withbias.root", "Multiplicity", kTRUE);
+
   // without bias calculation
-  mult->ApplyMinuitFit(hID, kFALSE, 0, kFALSE);
+  mult = AliMultiplicityCorrection::Open("nobias.root");
   baseold = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("baseold");
   
   // with bias calculation
-  mult->ApplyMinuitFit(hID, kFALSE, 0, kFALSE, 0, kTRUE);
+  mult = AliMultiplicityCorrection::Open("withbias.root");
   base = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("base");
  
   // relative error plots
@@ -176,7 +590,7 @@ TH1* GetChi2Bias(Float_t alpha)
   ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("hist1");
   ratio2->Reset();
 
-  for (Int_t t = 0; t<kMaxBins; t++)
+  for (Int_t t = 0; t<baseold->GetNbinsX(); t++)
   {
     Printf("Bin %d; Content: %f; Chi2 Error: %f; Bias: %f; In percent: %.2f %.2f", t+1, base->GetBinContent(t+1), baseold->GetBinError(t+1), base->GetBinError(t+1), (base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1, (base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
     if (base->GetBinContent(t+1) <= 0)
@@ -187,9 +601,9 @@ TH1* GetChi2Bias(Float_t alpha)
   
   //new TCanvas; base->Draw(); gPad->SetLogy();
 
-  new TCanvas;
+  c = new TCanvas;
   ratio1->GetYaxis()->SetRangeUser(0, 1);
-  ratio1->GetXaxis()->SetRangeUser(0, kMaxBins);
+  ratio1->GetXaxis()->SetRangeUser(0, 50);
   ratio1->Draw();
   ratio2->SetLineColor(2);
   ratio2->Draw("SAME");
@@ -243,7 +657,14 @@ void CheckBayesianBias()
   hist2->Draw("SAME");
 }
 
-void DataScan(Bool_t redo = kTRUE)
+void DrawUnfoldingLimit(Int_t hID, Float_t min, Float_t max)
+{
+  line = new TLine(kTrustLimits[hID], min, kTrustLimits[hID], max);
+  line->SetLineWidth(2);
+  line->Draw();
+}
+
+void DataScan(Int_t hID, Bool_t redo = kTRUE)
 {
   // makes a set of unfolded spectra and compares
   // don't forget FindUnfoldedLimit in plots.C
@@ -251,14 +672,13 @@ void DataScan(Bool_t redo = kTRUE)
   loadlibs();
 
   // files...
-  Bool_t mc = kTRUE;
+  Bool_t mc = 1;
   const char* fileNameMC = "multiplicityMC.root";
   const char* folder = "Multiplicity";
   const char* fileNameESD = "multiplicityESD.root";
   const char* folderESD = "Multiplicity";
-  Int_t hID = 0;  const Int_t kMaxBins = 35;
-  //Int_t hID = 1;  const Int_t kMaxBins = 60;
-  //Int_t hID = 2;  const Int_t kMaxBins = 70;
+  const Int_t kMaxBins = kBinLimits[hID];
+  
   AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx;
   Bool_t evaluteBias = kFALSE;
   
@@ -266,15 +686,15 @@ void DataScan(Bool_t redo = kTRUE)
   
   // chi2 range
   AliUnfolding::RegularizationType regBegin = AliUnfolding::kPol0;
-  AliUnfolding::RegularizationType regEnd = AliUnfolding::kPol0;
-  Float_t regParamBegin[] = { 0, 1, 10 }; 
-  Float_t regParamEnd[] = { 0, 40, 101 }; 
-  Float_t regParamStep[] = { 0, TMath::Sqrt(10), TMath::Sqrt(10) }; 
+  AliUnfolding::RegularizationType regEnd = AliUnfolding::kPol1;
+  Float_t regParamBegin[] = { 0, 1, 0.2, 3 }; 
+  Float_t regParamEnd[] = { 0, 11, 0.5, 31 }; 
+  Float_t regParamStep[] = { 0, 2, 2, TMath::Sqrt(10) }; 
 
   // bayesian range
-  Int_t iterBegin = 5;
-  Int_t iterEnd = 21;
-  Int_t iterStep = 5;
+  Int_t iterBegin = 40;
+  Int_t iterEnd = 41;
+  Int_t iterStep = 20;
   
   TList labels;
   
@@ -283,6 +703,16 @@ void DataScan(Bool_t redo = kTRUE)
   
   mult->SetMultiplicityESD(hID, esd->GetMultiplicityESD(hID));
   
+  // insert trigger efficiency in flat response matrix
+  AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
+  for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
+    mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
+  
+  mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
+  mult->FixTriggerEfficiencies(10);
+  
+  Float_t fraction0Generated = multTrigger->GetFraction0Generated(hID);
+  
   if (mc)
     mult->SetMultiplicityMC(hID, eventType, esd->GetMultiplicityMC(hID, eventType));
   
@@ -303,7 +733,7 @@ void DataScan(Bool_t redo = kTRUE)
     
       AliUnfolding::SetChi2Regularization(reg, regParam);
       
-      mult->ApplyMinuitFit(hID, kFALSE, eventType, kFALSE, 0, evaluteBias);
+      mult->ApplyMinuitFit(hID, kFALSE, eventType, fraction0Generated, kFALSE, 0, evaluteBias);
       
       file = TFile::Open(Form("datascan_%d.root", count), "RECREATE");
       mult->SaveHistograms();
@@ -347,17 +777,21 @@ void DataScan(Bool_t redo = kTRUE)
   
   TH1* residualSum = new TH1F("residualSum", ";;residual squared sum", count + 1, 0.5, count + 1.5);
   
+  Int_t allowedCount = count;
   count = 0;
   while (1)
   {
     mult = AliMultiplicityCorrection::Open(Form("datascan_%d.root", count++));
     if (!mult)
       break;
+    if (count > allowedCount+1)
+      break;
       
     hist = (TH1*) mult->GetMultiplicityESDCorrected(hID);
     hist->SetLineColor((count-1) % 4 + 1);
     hist->SetLineStyle((count-1) / 4 + 1);
     hist->GetXaxis()->SetRangeUser(0, kMaxBins);
+    hist->GetYaxis()->SetRangeUser(0.1, hist->GetMaximum() * 1.5);
     hist->SetStats(kFALSE);
     hist->SetTitle("");
     
@@ -366,6 +800,8 @@ void DataScan(Bool_t redo = kTRUE)
     c->cd(1);
     hist->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
     
+    DrawUnfoldingLimit(hID, 0.1, hist->GetMaximum());
+    
     if (mc)
     {
       TH2* mcHist2d = mult->GetMultiplicityMC(hID, eventType);
@@ -374,7 +810,7 @@ void DataScan(Bool_t redo = kTRUE)
       c->cd(1);
       mcHist->SetMarkerStyle(24);
       mcHist->Draw("P SAME");
-    
+      
       c->cd(2);
       // calculate ratio using only the error on the mc bin
       ratio = (TH1*) hist->Clone("ratio");
@@ -389,6 +825,8 @@ void DataScan(Bool_t redo = kTRUE)
       ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
       ratio->GetYaxis()->SetTitle("Ratio unfolded / MC");
       ratio->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
+    
+      DrawUnfoldingLimit(hID, 0.5, 1.5);
     }
     
     c->cd(3);
@@ -406,6 +844,8 @@ void DataScan(Bool_t redo = kTRUE)
     ratio->GetYaxis()->SetTitle("Ratio unfolded / unfolded reference case");
     ratio->DrawCopy((count == 1) ? "" : "SAME");
     
+    DrawUnfoldingLimit(hID, 0.5, 1.5);
+    
     c->cd(4);
     ratio = (TH1*) hist->Clone("ratio");
     for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
@@ -419,6 +859,8 @@ void DataScan(Bool_t redo = kTRUE)
     ratio->GetYaxis()->SetTitle("Relative error");
     ratio->DrawCopy((count == 1) ? "" : "SAME");
     
+    DrawUnfoldingLimit(hID, 0, 1);
+    
     c->cd(5);
     Float_t sum;
     residuals = mult->GetResiduals(hID, AliMultiplicityCorrection::kTrVtx, sum);
@@ -1855,6 +2297,7 @@ void Merge(Int_t n, const char** files, const char* output)
       name.Form("Multiplicity%d", i);
 
     TFile::Open(files[i]);
+    Printf("Loading from file %s", files[i]);
     data[i] = new AliMultiplicityCorrection(name, name);
     data[i]->LoadHistograms("Multiplicity");
     if (i > 0)
@@ -2266,6 +2709,7 @@ void BuildResponseFromTree(const char* fileName, const char* target)
 {
   //
   // builds several response matrices with different particle ratios (systematic study)
+  // particle-species study
   // 
   // WARNING doesn't work uncompiled, see test.C
 
@@ -2282,8 +2726,10 @@ void BuildResponseFromTree(const char* fileName, const char* target)
   Int_t secondaries = 0;
   Int_t doubleCount = 0;
 
-  for (Int_t num = 0; num < 9; num++)
+  for (Int_t num = 0; num < 7; num++)
   {
+    Printf("%d", num);
+  
     TString name;
     name.Form("Multiplicity_%d", num);
     AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(name, name);
@@ -2292,21 +2738,25 @@ void BuildResponseFromTree(const char* fileName, const char* target)
     for (Int_t i = 0; i < 4; i++)
       ratio[i] = 1;
 
-    switch (num)
-    {
-      case 1 : ratio[1] = 0.5; break;
-      case 2 : ratio[1] = 1.5; break;
-      case 3 : ratio[2] = 0.5; break;
-      case 4 : ratio[2] = 1.5; break;
-      case 5 : ratio[1] = 0.5; ratio[2] = 0.5; break;
-      case 6 : ratio[1] = 1.5; ratio[2] = 1.5; break;
-      case 7 : ratio[1] = 0.5; ratio[2] = 1.5; break;
-      case 8 : ratio[1] = 1.5; ratio[2] = 0.5; break;
-    }
-
+    if (num == 1)
+      ratio[1] = 0.7;
+    if (num == 2)
+      ratio[1] = 1.3;
+    
+    if (num == 3)
+      ratio[2] = 0.7;
+    if (num == 4)
+      ratio[2] = 1.3;
+    
+    if (num == 5)
+      ratio[3] = 0.7;
+    if (num == 6)
+      ratio[3] = 1.3;
+      
     for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
     {
       fParticleSpecies->GetEvent(i);
+      
 
       Float_t* f = fParticleSpecies->GetArgs();
 
@@ -2315,8 +2765,48 @@ void BuildResponseFromTree(const char* fileName, const char* target)
 
       for (Int_t j = 0; j < 4; j++)
       {
-        gene += ratio[j] * f[j+1];
-        meas += ratio[j] * f[j+1+4];
+        if (ratio[j] == 1)
+        {
+          gene += f[j+1];
+        }
+        else
+        {
+          for (Int_t k=0; k<f[j+1]; k++)
+          {
+            if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
+            {
+              gene += 1;
+            }
+            else if (ratio[j] > 1)
+            {
+              gene += 1;
+              if (gRandom->Uniform() < ratio[j] - 1)
+                gene += 1;
+            }
+          }
+        }
+          
+        if (ratio[j] == 1)
+        {
+          meas += f[j+1+4];
+        }
+        else
+        {
+          for (Int_t k=0; k<f[j+1+4]; k++)
+          {
+            if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
+            {
+              meas += 1;
+            }
+            else if (ratio[j] > 1)
+            {
+              meas += 1;
+              if (gRandom->Uniform() < ratio[j] - 1)
+                meas += 1;
+            }
+          }
+        }
+        
         tracks += f[j+1+4];
       }
 
@@ -2335,10 +2825,10 @@ void BuildResponseFromTree(const char* fileName, const char* target)
 
       //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
 
-      fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, 0, meas, meas, meas, meas);
+      fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, meas, meas, meas);
       // HACK all as kND = 1
-      fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 1, gene, gene, gene, gene, 0);
-      fMultiplicity->FillMeasured(f[0], meas, meas, meas, meas);
+      fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 1, gene, gene, gene, gene);
+      fMultiplicity->FillMeasured(f[0], meas, meas, meas);
     }
 
     //fMultiplicity->DrawHistograms();
@@ -2360,86 +2850,110 @@ void ParticleRatioStudy()
 {
   loadlibs();
 
-  for (Int_t num = 0; num < 9; num++)
+  for (Int_t num = 0; num < 7; num++)
   {
     TString target;
-    target.Form("chi2_species_%d.root", num);
-    correct("species.root", Form("Multiplicity_%d", num), "species.root", 1, 1, 0, 1e5, 0, target, "Multiplicity_0");
+    target.Form("chi2a_inel_species_%d.root", num);
+    
+    correct("compositions.root", Form("Multiplicity_%d", num), "multiplicityESD.root", 1, -1, 0, -1, 2, target);
   }
 }
 
-void MergeModifyCrossSection(const char* output = "multiplicityMC_xsection.root")
+void MergeCrossSection(Int_t xsectionID, const char* output = "multiplicityMC_xsection.root")
 {
-  const char* files[] = { "multiplicityMC_nd.root", "multiplicityMC_sd.root", "multiplicityMC_dd.root" };
+  const char* files[] = { "multiplicitySD.root", "multiplicityDD.root", "multiplicityND.root" };
 
   loadlibs();
 
   TFile::Open(output, "RECREATE");
   gFile->Close();
 
-  for (Int_t num = 0; num < 7; num++)
-  {
-    AliMultiplicityCorrection* data[3];
-    TList list;
-
-    Float_t ratio[3];
-    switch (num)
-    {
-      case 0: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.0; break;
-      case 1: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.0; break;
-      case 2: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 1.0; break;
-      case 3: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.5; break;
-      case 4: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 0.5; break;
-      case 5: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.5; break;
-      case 6: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 0.5; break;
-      default: return;
-    }
+  AliMultiplicityCorrection* data[3];
+  TList list;
 
-    for (Int_t i=0; i<3; ++i)
-    {
-      TString name;
-      name.Form("Multiplicity_%d", num);
-      if (i > 0)
-        name.Form("Multiplicity_%d_%d", num, i);
+  Float_t ref_SD = -1;
+  Float_t ref_DD = -1;
+  Float_t ref_ND = -1;
 
-      TFile::Open(files[i]);
-      data[i] = new AliMultiplicityCorrection(name, name);
-      data[i]->LoadHistograms("Multiplicity");
+  Float_t error_SD = -1;
+  Float_t error_DD = -1;
+  Float_t error_ND = -1;
 
-      // modify x-section
-      for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
-      {
-        data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
-        data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
-        data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
-        data[i]->GetMultiplicityNSD(j)->Scale(ratio[i]);
-      }
+  gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawSystematics.C"));
+  GetRelativeFractions(xsectionID, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
+  
+  for (Int_t i=0; i<3; ++i)
+  {
+    TString name;
+    name.Form("Multiplicity");
+    if (i > 0)
+      name.Form("Multiplicity_%d", i);
 
-      for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
-        data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
+    TFile::Open(files[i]);
+    data[i] = new AliMultiplicityCorrection(name, name);
+    data[i]->LoadHistograms("Multiplicity");
+  }
+  
+  // TODO is the under/overflow bin scaled as well? --> seems like, verify anyway!
 
-      for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
-        data[i]->GetCorrelation(j)->Scale(ratio[i]);
+  // calculating relative
+  Float_t sd = data[0]->GetMultiplicityINEL(0)->Integral(0, data[0]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+  Float_t dd = data[1]->GetMultiplicityINEL(0)->Integral(0, data[1]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+  Float_t nd = data[2]->GetMultiplicityINEL(0)->Integral(0, data[2]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+  Float_t total = nd + dd + sd;
+  
+  nd /= total;
+  sd /= total;
+  dd /= total;
+  
+  Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
+  
+  Float_t ratio[3];
+  ratio[0] = ref_SD / sd;
+  ratio[1] = ref_DD / dd;
+  ratio[2] = ref_ND / nd;
+  
+  Printf("SD=%.2f, DD=%.2f, ND=%.2f",ratio[0], ratio[1], ratio[2]);
+  
+  for (Int_t i=0; i<3; ++i)
+  {
+    // modify x-section
+    for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
+    {
+      data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
+      data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
+      data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
+      data[i]->GetMultiplicityNSD(j)->Scale(ratio[i]);
+    }
 
-      if (i > 0)
-        list.Add(data[i]);
+    for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
+    {
+      data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
+      data[i]->GetTriggeredEvents(j)->Scale(ratio[i]);
+      data[i]->GetNoVertexEvents(j)->Scale(ratio[i]);
     }
+    
+    for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
+      data[i]->GetCorrelation(j)->Scale(ratio[i]);
+
+    if (i > 0)
+      list.Add(data[i]);
+  }
 
-    printf("Case %d, %s: Entries in response matrix 3: ND: %.2f SD: %.2f DD: %.2f", num, data[0]->GetName(), data[0]->GetCorrelation(3)->Integral(), data[1]->GetCorrelation(3)->Integral(), data[2]->GetCorrelation(3)->Integral());
+  printf("Case %d, %s: Entries in response matrix: SD: %.2f DD: %.2f ND: %.2f", xsectionID, data[0]->GetName(), data[0]->GetCorrelation(0)->Integral(), data[1]->GetCorrelation(0)->Integral(), data[2]->GetCorrelation(0)->Integral());
 
-    data[0]->Merge(&list);
+  data[0]->Merge(&list);
 
-    Printf(" Total: %.2f", data[0]->GetCorrelation(3)->Integral());
+  Printf(" Total: %.2f", data[0]->GetCorrelation(0)->Integral());
 
-    TFile::Open(output, "UPDATE");
-    data[0]->SaveHistograms();
-    gFile->Close();
+  TFile::Open(output, "RECREATE");
+  data[0]->SaveHistograms();
+  gFile->Close();
 
-    list.Clear();
+  list.Clear();
 
-    for (Int_t i=0; i<3; ++i)
-      delete data[i];
-  }
+  for (Int_t i=0; i<3; ++i)
+    delete data[i];
 }
 
 void Rebin(const char* fileName = "multiplicityMC_3M.root", Int_t corrMatrix = 3)
@@ -2597,3 +3111,219 @@ void EvaluateRegularizationEffect(Int_t step, const char* fileNameRebinned = "mu
     canvas->SaveAs(Form("%s.eps", canvas->GetName()));
   }
 }
+
+void MergeDistributions()
+{
+  loadlibs();
+
+  const char* dir1 = "run000104824-52_pass4";
+  const char* dir2 = "run000104867_90_92_pass4";
+  
+  for (Int_t evType = 0; evType < 2; evType++)
+  {
+    Printf("%d", evType);
+    
+    const char* evTypeStr = ((evType == 0) ? "inel" : "nsd");
+
+    const char* id = "chi2a";
+    //const char* id = "bayes";
+    
+    TString suffix;
+    suffix.Form("/all/spd/%s_", id);
+    if (evType == 1)
+      suffix.Form("/v0and/spd/%s_", id);
+  
+    TString file1, file2;
+    file1.Form("%s%s%%s.root", dir1, suffix.Data());
+    file2.Form("%s%s%%s.root", dir2, suffix.Data());
+    
+    if (1)
+    {
+      const char* files[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr) };
+      Merge(2, files, Form("merged/%s_%s.root", id, evTypeStr));
+    }
+    else
+    {
+      AliMultiplicityCorrection* mult1 = AliMultiplicityCorrection::Open(Form(file1.Data(), evTypeStr));
+      AliMultiplicityCorrection* mult2 = AliMultiplicityCorrection::Open(Form(file2.Data(), evTypeStr));
+      
+      AliMultiplicityCorrection* target = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+      
+      for (Int_t etaRange = 0; etaRange < 3; etaRange++)
+      {
+        hist1 = mult1->GetMultiplicityESDCorrected(etaRange);
+        hist2 = mult2->GetMultiplicityESDCorrected(etaRange);
+        targetHist = target->GetMultiplicityESDCorrected(etaRange);
+        
+        //hist1->Scale(1.0 / hist1->Integral());
+        //hist2->Scale(1.0 / hist2->Integral());
+        
+        //targetHist->SetBinContent(1, hist1->GetBinContent(1) * 0.5 + hist2->GetBinContent(1) * 0.5);
+        targetHist->SetBinContent(1, hist1->GetBinContent(1) + hist2->GetBinContent(1));
+        for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
+        {
+          if (hist1->GetBinError(i) > 0 && hist2->GetBinError(i) > 0)
+          {
+            Float_t w1 = 1.0 / hist1->GetBinError(i) / hist1->GetBinError(i);
+            Float_t w2 = 1.0 / hist2->GetBinError(i) / hist2->GetBinError(i);
+            
+            Float_t average = (hist1->GetBinContent(i) * w1 + hist2->GetBinContent(i) * w2) / (w1 + w2);
+            
+            //targetHist->SetBinContent(i, average);
+            //targetHist->SetBinError(i, TMath::Max(mult1->GetBinError(i), mult2->GetBinError(i)));
+            
+            targetHist->SetBinContent(i, hist1->GetBinContent(i) + hist2->GetBinContent(i));
+            targetHist->SetBinError(i, hist1->GetBinError(i) + hist2->GetBinError(i));
+          }
+        }
+      }
+      
+      file = TFile::Open(Form("merged/%s_%s.root", id, evTypeStr), "RECREATE");
+      target->SaveHistograms();
+      file->Close();
+    }
+
+    const char* files2[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr), Form("merged/%s_%s.root", id, evTypeStr) };
+    drawAll(files2, (evType == 1), kTRUE);
+  }
+}
+
+void CompareDistributions(Int_t evType)
+{
+  loadlibs();
+  gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C"));
+
+  const char* dir1 = "run000104824-52_pass4";
+  const char* dir2 = "run000104867_90_92_pass4";
+  
+  const char* evTypeStr = (evType == 0) ? "inel" : "nsd";
+
+  const char* suffix = "/all/spd/chi2a_";
+  if (evType == 1)
+    suffix = "/v0and/spd/chi2a_";
+  
+  TString file1, file2;
+  file1.Form("%s%s%s.root", dir1, suffix, evTypeStr);
+  file2.Form("%s%s%s.root", dir2, suffix, evTypeStr);
+    
+  for (Int_t hID = 0; hID < 3; hID++)
+    CompareQualityHists(file1, file2, Form("Multiplicity/fMultiplicityESDCorrected%d", hID), 1, 1);
+}
+
+void StatisticalUncertainty(Int_t methodType, Int_t etaRange, Bool_t mc = kFALSE)
+{
+  loadlibs();
+  
+  AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
+  AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open("multiplicityESD.root");
+  AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
+
+  TH1* mcHist = esd->GetMultiplicityMB(etaRange)->ProjectionY("mymc", 1, 1);
+
+  Int_t geneLimits[] = { 0, 0, 0 };
+
+  LoadAndInitialize(mult, esd, multTrigger, etaRange, kFALSE, geneLimits);
+
+  AliUnfolding::SetNbins(kBinLimits[etaRange], geneLimits[etaRange]);
+  AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 5);
+  AliUnfolding::SetBayesianParameters(1, 10);
+
+  // background subtraction
+  Int_t background = 0;
+  
+  //background = 1091 + 4398; // MB1 for run 104824 - 52
+  background = 417 + 1758; // MB1 for run 104867 - 92
+  
+  //background = 20 + 0; // V0AND for run 104824 - 52
+  //background = 10 + 0; // V0AND for run 104867 - 92
+  
+  Printf("NOTE: Subtracting %d background events", background);
+
+  Int_t zeroBin = mult->GetTriggeredEvents(etaRange)->GetBinContent(1) - background - mult->GetMultiplicityESD(etaRange)->Integral(0, 2, 1, 1);
+  
+  TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
+  
+  new TCanvas; errorMeasured->Draw();
+  new TCanvas; 
+
+  AliPWG0Helper::NormalizeToBinWidth(mult->GetMultiplicityESDCorrected(etaRange));
+  mult->GetMultiplicityESDCorrected(etaRange)->Draw();
+  
+  if (mc)
+  {
+    AliPWG0Helper::NormalizeToBinWidth(mcHist);
+    mcHist->Scale(mult->GetMultiplicityESDCorrected(etaRange)->Integral() / mcHist->Integral()); 
+    mcHist->SetLineColor(2);
+    mcHist->Draw("SAME");
+  }
+  gPad->SetLogy();
+  
+  TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
+
+  TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
+
+  if (mc)
+  {
+    TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+    DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
+  }
+
+  TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
+  errorResponse->Write();
+  errorMeasured->Write();
+  errorBoth->Write();
+  file->Close();
+}
+
+void ChangeResponseMatrixEfficiency(Bool_t reduceEff = kTRUE, Float_t factor = 0.01625)
+{
+  loadlibs();
+  
+  AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicity.root");
+  
+  for (Int_t etaR = 0; etaR < 3; etaR++)
+  {
+    TH3* corr = mult->GetCorrelation(etaR);
+    TH3* corrOld = (TH3*) corr->Clone("corrOld");
+    
+    corr->Reset();
+    
+    Int_t total = 0;
+    Int_t change = 0;
+    for (Int_t x=0; x<=corr->GetNbinsX()+1; x++)
+    {
+      for (Int_t y=0; y<=corr->GetNbinsY()+1; y++)
+      {
+        Printf("%d", y);
+        for (Int_t z=0; z<=corr->GetNbinsZ()+1; z++)
+        {
+          Float_t tracklets = corr->GetZaxis()->GetBinCenter(z);
+          
+          for (Int_t n=0; n<corrOld->GetBinContent(x, y, z); n++)
+          {
+            total += tracklets;
+            Float_t trackletsNew = tracklets;
+            
+            for (Int_t i=0; i<tracklets; i++)
+            {
+              if (gRandom->Uniform() < factor)
+              {
+                trackletsNew += (reduceEff) ? -1 : 1;
+                change += (reduceEff) ? -1 : 1;
+              }
+            }
+                
+            corr->Fill(corr->GetXaxis()->GetBinCenter(x), corr->GetYaxis()->GetBinCenter(y), trackletsNew);
+          }
+        }
+      }
+    }
+    
+    Printf("%d: Changed %d out of %d total tracklets", etaR, change, total);
+  }
+
+  file = TFile::Open("out.root", "RECREATE");
+  mult->SaveHistograms();
+  file->Close();
+}
+
index 7b933074e709615801916274620eb80d86f8048f..c68d706acbe9d5fcdb77bdd178789eb88d2ed366 100644 (file)
@@ -71,6 +71,10 @@ const char* GetMultLabel(Int_t etaR = -1, Bool_t trueM = kTRUE)
        {
                tmpStr += " (full phase space)";
        }
+       else if (etaR == 2)
+       {
+               tmpStr += " in |#eta| < 1.3";
+       }
        else
                tmpStr += Form(" in |#eta| < %.1f", (etaR+1)* 0.5);
        return Form("%s", tmpStr.Data());
@@ -92,20 +96,20 @@ void Smooth(TH1* hist, Int_t windowWidth = 20)
   delete clone;
 }
 
-void responseMatrixPlot(const char* fileName = 0)
+TH1* responseMatrixPlot(const char* fileName = 0, const char* folder = "Multiplicity")
 {
   loadlibs();
 
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
 
   if (fileName == 0)
     fileName = correctionFile;
   
   TFile::Open(fileName);
-  mult->LoadHistograms("Multiplicity");
+  mult->LoadHistograms();
 
   // empty under/overflow bins in x, otherwise Project3D takes them into account
-  TH1* hist = mult->GetCorrelation(etaRange);
+  TH2* hist = (TH2*) mult->GetCorrelation(etaRange);
   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
   {
     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
@@ -114,9 +118,29 @@ void responseMatrixPlot(const char* fileName = 0)
       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
     }
   }
-  hist = ((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy");
+  
+  hist = (TH2*) (((TH3*) mult->GetCorrelation(etaRange))->Project3D("zy"));
   hist->SetStats(kFALSE);
+  
+  if (0)
+  {
+    // normalize
+    // normalize correction for given nPart
+    for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
+    {
+      Double_t sum = hist->Integral(i, i, 1, hist->GetNbinsY());
+      if (sum <= 0)
+        continue;
 
+      for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
+      {
+        // npart sum to 1
+        hist->SetBinContent(i, j, hist->GetBinContent(i, j) / sum);// * correlation->GetXaxis()->GetBinWidth(i));
+        hist->SetBinError(i, j, hist->GetBinError(i, j) / sum);
+      }
+    }
+  }
+  
   hist->SetTitle(Form(";%s;%s;Entries", GetMultLabel(), GetMultLabel(etaRange, kFALSE)));
   hist->GetXaxis()->SetRangeUser(0, longDisplayRange);
   hist->GetYaxis()->SetRangeUser(0, longDisplayRange);
@@ -132,6 +156,8 @@ void responseMatrixPlot(const char* fileName = 0)
   hist->Draw("COLZ");
 
   canvas->SaveAs("responsematrix.eps");
+  
+  return hist;
 }
 
 void multPythiaPhojet()
@@ -543,13 +569,55 @@ TCanvas* Draw2ResultRatios(TH1* mcHist, TH1* result1, TH1* result2, TH1* ratio2,
   return canvas;
 }
 
-TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE)
+void DrawEfficiencyChange()
+{
+  loadlibs();
+
+  const char* fileName = "chi2a_inel.root";
+
+  base = AliMultiplicityCorrection::Open(Form("spd/%s", fileName));
+  low = AliMultiplicityCorrection::Open(Form("spd-loweff/%s", fileName));
+  high = AliMultiplicityCorrection::Open(Form("spd-higheff/%s", fileName));
+  
+  const char* legendStrings[] = { "low efficiency", "high efficiency" };
+  
+  file = TFile::Open("systunc_detectorefficiency.root", "RECREATE");
+  
+  for (Int_t etaR = 1; etaR < 2; etaR++)
+  {
+    base->GetMultiplicityESDCorrected(etaR)->Scale(1.0 / base->GetMultiplicityESDCorrected(etaR)->Integral(2, base->GetMultiplicityESDCorrected(etaR)->GetNbinsX()));
+  
+    TH1* hists[2];
+    hists[0] = low->GetMultiplicityESDCorrected(etaR);
+    hists[0]->Scale(1.0 / hists[0]->Integral(2, hists[0]->GetNbinsX()));
+
+    if (high)
+    {
+      hists[1] = high->GetMultiplicityESDCorrected(etaR);
+      hists[1]->Scale(1.0 / hists[1]->Integral(2, hists[1]->GetNbinsX()));
+    }
+    
+    largestError = (TH1*) hists[0]->Clone(Form("detectorefficiency_%d", etaR));
+    largestError->Reset();
+  
+    DrawRatio(base->GetMultiplicityESDCorrected(etaR), (high) ? 2 : 1, hists, Form("eff_%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError);
+    
+    largestError->Write();
+    
+    new TCanvas;
+    largestError->DrawCopy();
+  }
+  
+  file->Close();
+}
+
+TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0, Bool_t errors = kFALSE, TH1* largestErrorLow = 0, TH1* largestErrorHigh = 0)
 {
   // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
   // a systematic effect
 
   // normalize results
-  result->Scale(1.0 / result->Integral(2, 200));
+  //result->Scale(1.0 / result->Integral(1, 200));
 
   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 500);
   canvas->SetTopMargin(0.05);
@@ -572,13 +640,13 @@ TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString eps
 
   TLegend* legend = new TLegend(0.2, 0.7, 0.7, 0.93);
   legend->SetFillColor(0);
-  legend->SetTextSize(0.04);
+  //legend->SetTextSize(0.04);
   if (nResultSyst > 6)
     legend->SetNColumns(2);
 
   for (Int_t n=0; n<nResultSyst; ++n)
   {
-    resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
+    //resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(1, 200));
 
     // calculate ratio
     TH1* ratio = (TH1*) result->Clone("ratio");
@@ -602,6 +670,26 @@ TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString eps
 
     if (legendStrings && legendStrings[n])
       legend->AddEntry(ratio, legendStrings[n], "L");
+      
+    /*
+    s = new TSpline3(ratio);
+    s->SetNpx(5);
+    s->SetLineColor(kRed);
+    s->Draw("same");
+    */
+    
+    if (largestErrorLow && largestErrorHigh)
+      for (Int_t i=1; i<=ratio->GetNbinsX(); ++i)
+      {
+        if (ratio->GetBinContent(i) - 1 > 0)
+          largestErrorHigh->SetBinContent(i, TMath::Max(ratio->GetBinContent(i) - 1, largestErrorHigh->GetBinContent(i)));
+        if (ratio->GetBinContent(i) - 1 < 0)
+          largestErrorLow->SetBinContent(i, TMath::Min(ratio->GetBinContent(i) - 1, largestErrorLow->GetBinContent(i)));
+      }
+
+    if (largestErrorLow && !largestErrorHigh)
+      for (Int_t i=1; i<=ratio->GetNbinsX(); ++i)
+        largestErrorLow->SetBinContent(i, TMath::Max(TMath::Abs(ratio->GetBinContent(i) - 1), largestErrorLow->GetBinContent(i)));
 
     // get average of ratio
     Float_t sum = 0;
@@ -629,7 +717,6 @@ TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString eps
   line->Draw();
 
   canvas->SaveAs(epsName);
-  canvas->SaveAs(Form("%s.gif", epsName.Data()));
 
   return canvas;
 }
@@ -851,16 +938,34 @@ void DrawResiduals(const char* fileName, const char* epsName)
 {
   loadlibs();
 
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  TFile::Open(fileName);
-  mult->LoadHistograms("Multiplicity");
+  AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileName);
 
   TH1* measured = mult->GetMultiplicityESD(etaRange)->ProjectionY("myesd", 1, 1);
-  TH1* unfoldedFolded = mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1);
+  
+  if (0)
+  {
+    // test how far we are from a normal exponential in the unfolded
+    corrected = mult->GetMultiplicityESDCorrected(etaRange);
+    new TCanvas; corrected->DrawCopy(); gPad->SetLogy();
+    
+    expo = new TF1("exp", "expo(0)", 0, 50);
+    //expo->SetParameters(10, -0.18);
+    expo->SetParameters(9.96, -0.176);
+    expo->Draw("SAME");  
+    
+    for (Int_t i=21; i<=50; i++)
+      corrected->SetBinContent(i, expo->Eval(corrected->GetXaxis()->GetBinCenter(i)));
+    
+    corrected->SetLineColor(2);
+    corrected->Draw("SAME");
+  }
+  
+  TH1* unfoldedFolded = mult->GetConvoluted(etaRange, AliMultiplicityCorrection::kINEL);
+  //mult->CalculateMultiplicityESD(mult->GetMultiplicityESDCorrected(etaRange), etaRange)->ProjectionY("myfolded", 1, 1);
   
   // normalize
-  unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1));
-  unfoldedFolded->Scale(measured->Integral(2, displayRange+1));
+  //unfoldedFolded->Scale(1.0 / unfoldedFolded->Integral(2, displayRange+1));
+  //unfoldedFolded->Scale(measured->Integral(2, displayRange+1));
 
   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
   canvas->Range(0, 0, 1, 1);
@@ -928,10 +1033,12 @@ void DrawResiduals(const char* fileName, const char* epsName)
   residual->Add(unfoldedFolded, -1);
 
   // projection
-  TH1* residualHist = new TH1F("residualHist", ";", 11, -3, 3);
+  TH1* residualHist = new TH1F("residualHist", ";", 16, -3, 3);
   residualHist->Sumw2();
 
   Float_t chi2 = 0;
+  Float_t chi2_hump = 0;
+  
   for (Int_t i=1; i<=displayRange+1; ++i)
   {
     if (measured->GetBinError(i) > 0)
@@ -940,6 +1047,9 @@ void DrawResiduals(const char* fileName, const char* epsName)
       residual->SetBinError(i, 1);
 
       residualHist->Fill(residual->GetBinContent(i));
+      if (i >= 15 && i <= 23)
+        chi2_hump += residual->GetBinContent(i) * residual->GetBinContent(i);
+
       chi2 += residual->GetBinContent(i) * residual->GetBinContent(i);
     }
     else
@@ -950,6 +1060,7 @@ void DrawResiduals(const char* fileName, const char* epsName)
   }
   
   Printf("chi2 / ndf = %f / %d = %f", chi2, displayRange+1, chi2 / (displayRange+1));
+  Printf("chi2 (hump) / ndf = %f / %d = %f", chi2_hump, 23-15+1, chi2_hump / (23-15+1));
 
   residual->GetYaxis()->SetTitle("Residuals:   (1/e) (M - R  #otimes U)");
   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
@@ -1476,9 +1587,9 @@ void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
 
   // SPD TPC
   //const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
-  const char* fileName[] = { "LHC09b9_0.9TeV_0T/mb1/spd/multiplicity.root", "LHC09b8_0.9TeV_0.5T/mb1/tpc/multiplicity.root" };
+  //const char* fileName[] = { "LHC09b9_0.9TeV_0T/mb1/spd/multiplicity.root", "LHC09b8_0.9TeV_0.5T/mb1/tpc/multiplicity.root" };
   //const char* fileName[] = { "spd/multiplicity.root", "tpc/multiplicity.root" };
-  //const char* fileName[] = { "multiplicity.root", "multiplicity.root" };
+  const char* fileName[] = { "multiplicityparticle-efficiency.root", "multiplicity.root" };
   Float_t etaRangeArr[] = {0.49, 0.9};
   const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
 
@@ -1493,7 +1604,7 @@ void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
   
   TLegend* legends[2];
   
-  for (Int_t loop=0; loop<2; ++loop)
+  for (Int_t loop=0; loop<1; ++loop)
   {
     Printf("%s", fileName[loop]);
 
@@ -1575,12 +1686,12 @@ void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
       TH1* genePt = gene->Project3D(Form("z_%d", i));
       TH1* measPt = meas->Project3D(Form("z_%d", i));
       
-      if (i == 2)
-      {
-        Printf("WARNING: Rebinning for protons!");
-        genePt->Rebin(2);
-        measPt->Rebin(2);
-      }
+//       if (i == 2)
+//       {
+//         Printf("WARNING: Rebinning for protons!");
+//         genePt->Rebin(2);
+//         measPt->Rebin(2);
+//       }
 
       genePt->Sumw2();
       measPt->Sumw2();
@@ -1702,6 +1813,8 @@ void EfficiencySpecies(Bool_t addDecayStopped = kFALSE)
     legend->Draw();
     canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));  
   }
+  
+  return;
 
   canvas->cd();
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
@@ -1751,50 +1864,19 @@ void ParticleSpeciesComparison()
   TH1* mc = 0;
   
   // loop over cases (normal, enhanced/reduced ratios)
-  Int_t nMax = 9;
+  Int_t nMax = 7;
   for (Int_t i = 0; i<nMax; ++i)
   {
-    AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2_species_%d.root", i), Form("Multiplicity_%d", i));
-    if (i == 0)
-      mc = mult->GetMultiplicityVtx(etaRange)->ProjectionY("mymchist", 1, 1);
+    AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(Form("chi2a_inel_species_%d.root", i), "Multiplicity");
     results[i] = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone(Form("result_%d", i));
   }
 
-  DrawResultRatio(mc, results[0], "ParticleSpeciesComparison1_1.eps");
-
   for (Int_t i=1; i<=results[0]->GetNbinsX(); i++)
-  {
     results[0]->SetBinError(i, 0);
-    mc->SetBinError(i, 0);
-  }
 
-  const char* legendStrings[] = { "K #times 0.5", "K #times 1.5", "p #times 0.5", "p #times 1.5", "K #times 0.5, p #times 0.5", "K #times 1.5, p #times 1.5", "K #times 0.5, p #times 1.5", "K #times 1.5, p #times 0.5" };
+  const char* legendStrings[] = { "K #times 0.7", "K #times 1.3", "p #times 0.7", "p #times 1.3", "others #times 0.7", "others #times 1.3",  };
 
   DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
-
-  //not valid: draw chi2 uncertainty on top!
-  /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
-  TH1* errorHist = (TH1*) gFile->Get("errorBoth");
-  errorHist->SetLineColor(1);
-  errorHist->SetLineWidth(2);
-  TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
-  for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
-  {
-    errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
-    errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
-  }
-
-  errorHist->DrawCopy("SAME");
-  errorHist2->DrawCopy("SAME");*/
-
-  //canvas->SaveAs(canvas->GetName());
-
-  DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
-
-  //errorHist->DrawCopy("SAME");
-  //errorHist2->DrawCopy("SAME");
-
-  //canvas2->SaveAs(canvas2->GetName());
 }
 
 /*void ParticleSpeciesComparison2()
@@ -1874,6 +1956,30 @@ TH1* Invert(TH1* eff)
   return corr;
 }
 
+void CompareVertexRecoEfficiencies()
+{
+  loadlibs();
+  
+  const char* file1 = "LHC10a12_run10482X";
+  const char* file2 = "LHC10a14_run10482X_Phojet";
+  
+  const char* suffix = "/all/spd/multiplicityMC_xsection.root";
+  
+  hist1 = AliMultiplicityCorrection::Open(Form("%s%s", file1, suffix))->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB);
+  hist2 = AliMultiplicityCorrection::Open(Form("%s%s", file2, suffix))->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB);
+  
+  ratio = (TH1*) hist1->Clone("ratio");
+  ratio->Divide(hist2);
+  
+  new TCanvas;
+  hist1->Draw();
+  hist2->SetLineColor(2);
+  hist2->Draw("SAME");
+  
+  new TCanvas;
+  ratio->Draw();
+}
+
 void TriggerVertexCorrection()
 {
   //
@@ -1889,6 +1995,8 @@ void TriggerVertexCorrection()
   TH1* corrINEL = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kINEL));
   TH1* corrNSD = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
   TH1* corrMB   = Invert(mult->GetEfficiency(etaRange, AliMultiplicityCorrection::kMB));
+  
+  TH1* triggerEff = Invert(mult->GetTriggerEfficiency(etaRange, AliMultiplicityCorrection::kNSD));
 
   TCanvas* canvas = new TCanvas("TriggerVertexCorrection", "TriggerVertexCorrection", 800, 500);
   gPad->SetGridx();
@@ -1914,9 +2022,14 @@ void TriggerVertexCorrection()
   corrNSD->SetMarkerColor(4);
   corrNSD->Draw("SAME PE");
   
-  Printf("       MB  INEL  NSD");
-  Printf("bin 0: %f %f %f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1));
-  Printf("bin 1: %f %f %f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2));
+  triggerEff->SetLineColor(6);
+  triggerEff->SetMarkerStyle(25);
+  triggerEff->SetMarkerColor(6);
+  //triggerEff->Draw("SAME PE");
+  
+  Printf("       MB  INEL  NSD  TRIGINEL");
+  Printf("bin 0: %.2f %.2f %.2f %.2f", corrMB->GetBinContent(1), corrINEL->GetBinContent(1), corrNSD->GetBinContent(1), triggerEff->GetBinContent(1));
+  Printf("bin 1: %.2f %.2f %.2f %.2f", corrMB->GetBinContent(2), corrINEL->GetBinContent(2), corrNSD->GetBinContent(2), triggerEff->GetBinContent(2));
 
   TLegend* legend = new TLegend(0.3, 0.6, 0.85, 0.85);
   legend->SetFillColor(0);
@@ -1930,61 +2043,9 @@ void TriggerVertexCorrection()
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 }
 
-void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
+void DrawStatisticalUncertainty(const char* fileName = "StatisticalUncertainty.root")
 {
-  loadlibs();
-
-  TFile::Open(correctionFile);
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  mult->LoadHistograms("Multiplicity");
-
-  TFile::Open(measuredFile);
-  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
-  mult2->LoadHistograms("Multiplicity");
-
-  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
-
-  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc", 1, 1);
-
-  AliUnfolding::SetNbins(70, 70);
-  //AliUnfolding::SetNbins(35, 35);
-  AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 1e3);
-  AliUnfolding::SetBayesianParameters(1, 10);
-  
-  TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
-  
-  new TCanvas; errorMeasured->Draw();
-  new TCanvas; 
-  
-  mult->GetMultiplicityESDCorrected(etaRange)->Scale(1.0 / mult->GetMultiplicityESDCorrected(etaRange)->Integral());
-  mult->GetMultiplicityESDCorrected(etaRange)->Draw();
-  mcHist->Scale(1.0 / mcHist->Integral()); 
-  mcHist->SetLineColor(2);
-  mcHist->Draw("SAME");
-  gPad->SetLogy();
-  
-  return;
-  
-  TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
-
-  TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
-
-  if (!mc)
-  {
-    TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-    DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
-  }
-
-  TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
-  errorResponse->Write();
-  errorMeasured->Write();
-  errorBoth->Write();
-  file->Close();
-}
-
-void DrawStatisticalUncertainty()
-{
-  TFile::Open("StatisticalUncertainty.root");
+  TFile::Open(fileName);
   
   errorResponse = (TH1*) gFile->Get("errorResponse");
   errorMeasured = (TH1*) gFile->Get("errorMeasured");
@@ -1998,7 +2059,7 @@ void DrawStatisticalUncertainty()
 
   errorResponse->SetLineColor(1);
   errorResponse->GetXaxis()->SetRangeUser(0, longDisplayRange);
-  errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
+  errorResponse->GetYaxis()->SetRangeUser(0, 1);
   errorResponse->SetStats(kFALSE);
   errorResponse->SetTitle(";true multiplicity;Uncertainty");
 
@@ -2010,6 +2071,10 @@ void DrawStatisticalUncertainty()
   errorBoth->SetLineColor(4);
   errorBoth->Draw("SAME");
 
+  errorResponse->Scale(1.0 / sqrt(2));
+  errorMeasured->Scale(1.0 / sqrt(2));
+  errorBoth->Scale(1.0 / sqrt(2));
+  
   Printf("Average errorResponse: %f", errorResponse->Integral(2, displayRange) / (displayRange - 1));
   Printf("Average errorMeasured: %f", errorMeasured->Integral(2, displayRange) /  (displayRange - 1));
   Printf("Average errorBoth: %f", errorBoth->Integral(2, displayRange) /  (displayRange - 1));
@@ -3751,12 +3816,14 @@ void FindUnfoldedLimit()
   new TCanvas; esd_proj->DrawCopy();
   
   TH1* percentage = (TH1*) (esd_proj->Clone("percentage"));
+  percentage->SetTitle("percentage");
   percentage->Reset();
   
   for (Int_t i=1; i<=esd_proj->GetNbinsX(); i++)
     if (esd_proj->GetBinContent(i) > 0)
       esd_proj->SetBinContent(i, 1);
   
+  Int_t limit = -1;
   for (Int_t i=1; i<=percentage->GetNbinsX(); i++)
   {
     TH1* binResponse = corr->ProjectionY("proj", i, i);
@@ -3765,17 +3832,30 @@ void FindUnfoldedLimit()
     binResponse->Scale(1.0 / binResponse->Integral());
     binResponse->Multiply(esd_proj);
     //new TCanvas; binResponse->Draw();
-    percentage->SetBinContent(i, binResponse->Integral());
+    Float_t value = binResponse->Integral();
+    percentage->SetBinContent(i, value);
+    if (limit == -1 && value < 0.9)
+      limit = percentage->GetXaxis()->GetBinCenter(i-1);
     //return;
   }
   
+  Printf("Limit is %d", limit);
+  
   new TCanvas; percentage->Draw();
+  
   new TCanvas;
   mc = esd->GetMultiplicityVtx(etaRange)->ProjectionY("mc", 1, esd->GetMultiplicityVtx(etaRange)->GetNbinsX());
   mc->SetLineColor(2);
   mc->Draw("");
 }
    
+void FindUnfoldedLimitAll()
+{
+  loadlibs();
+  for (etaRange = 0; etaRange <= 2; etaRange++)
+    FindUnfoldedLimit();
+}
+   
 void CompareUnfoldedWithMC()
 {
   loadlibs();
@@ -3792,3 +3872,390 @@ void CompareUnfoldedWithMC()
   mcHist->Draw("SAME");
   gPad->SetLogy();
 }
+
+void PaperNumbers()
+{
+  const char* label[] = { "SD", "DD", "ND" };
+  const char* files[] = { "multiplicity_syst_sd.root", "multiplicity_syst_dd.root", "multiplicity_syst_nd.root" };
+  
+  loadlibs();
+  
+  Printf("vertex reco");
+  
+  Float_t oneParticle[3];
+  for (Int_t i=0; i<3; i++)
+  {
+    AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(files[i]);
+    
+    eff = mult->GetEfficiency(2, AliMultiplicityCorrection::kMB);
+    //eff = mult->GetTriggerEfficiency(2);
+    
+    oneParticle[i] = 100.0 * eff->GetBinContent(2);
+    
+    Printf("%s: %.2f", label[i], oneParticle[i]);
+  }
+  
+  Float_t ua5_SD = 0.153;
+  Float_t ua5_DD = 0.080;
+  Float_t ua5_ND = 0.767;
+  
+  Float_t vtxINELUA5 = ua5_SD * oneParticle[0] + ua5_DD * oneParticle[1] + ua5_ND * oneParticle[2];
+  Float_t vtxNSDUA5  = (ua5_DD * oneParticle[1] + ua5_ND * oneParticle[2]) / (ua5_DD + ua5_ND);
+  
+  Printf("INEL (UA5)  = %.1f", vtxINELUA5);
+  Printf("NSD (UA5)  = %.1f", vtxNSDUA5);
+  
+  Printf("total for >= 2 charged tracks");
+  
+  AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicity.root");
+  
+  vtx = mult->GetMultiplicityVtx(2)->ProjectionY("vtx", 1, mult->GetMultiplicityVtx(2)->GetNbinsX(), "e");
+  all = mult->GetMultiplicityINEL(2)->ProjectionY("all", 1, mult->GetMultiplicityINEL(2)->GetNbinsX(), "e"); 
+  
+  Int_t begin = vtx->GetXaxis()->FindBin(2);
+  Int_t end = vtx->GetNbinsX();
+  
+  Float_t above2 = vtx->Integral(begin, end) / all->Integral(begin, end);
+  
+  Printf("%.1f", 100.0 * above2);
+}
+
+void DrawRawEta()
+{
+  TFile::Open(measuredFile);
+  
+  Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
+  
+  for (Int_t i=0; i<3; i++)
+  {
+    hist = dynamic_cast<TH1*> (gFile->Get(Form("fEta_%d", i)));
+    
+    hist->GetYaxis()->SetRangeUser(0, 8);
+    hist->SetLineColor(colors[i]);
+    hist->SetStats(kFALSE);
+    hist->Draw((i == 0) ? "HIST" : "HIST SAME");
+  }   
+  
+  gPad->SetGridx();
+  gPad->SetGridy();
+}
+
+void CrossSectionUncertainties(Int_t xsectionID, Int_t eventType)
+{
+  const char* files[] = { "multiplicitySD.root", "multiplicityDD.root", "multiplicityND.root" };
+
+  loadlibs();
+
+  AliMultiplicityCorrection* data[3];
+
+  Float_t ref_SD = -1;
+  Float_t ref_DD = -1;
+  Float_t ref_ND = -1;
+
+  Float_t error_SD = -1;
+  Float_t error_DD = -1;
+  Float_t error_ND = -1;
+
+  gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawSystematics.C"));
+  GetRelativeFractions(xsectionID, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
+  
+  TH1* unc[9*3];
+  TH1* unc2[9*3];
+  const char* legendStrings[8];
+  
+  for (Int_t i=0; i<9; i++)
+  {
+    Float_t factorSD = 0;
+    Float_t factorDD = 0;
+    
+    if (i > 0 && i < 3)
+      factorSD = (i % 2 == 0) ? 1 : -1;
+    else if (i >= 3 && i < 5)
+      factorDD = (i % 2 == 0) ? 1 : -1;
+    else if (i >= 5 && i < 9)
+    {
+      factorSD = ((i % 2 == 0) ? 1.0 : -1.0) / TMath::Sqrt(2);
+      if (i == 5 || i == 6)
+        factorDD = factorSD;
+      else
+        factorDD = -factorSD;
+    }
+    
+    Float_t scalesSD = ref_SD + factorSD * error_SD;
+    Float_t scalesDD = ref_DD + factorDD * error_DD;
+    Float_t scalesND = 1.0 - scalesDD[i] - scalesSD[i];
+    
+    str = new TString;
+    str->Form("Case %d: SD: %.2f DD: %.2f ND: %.2f", i, scalesSD, scalesDD, scalesND);
+    Printf("%s", str->Data());
+    if (i > 0)
+      legendStrings[i-1] = str->Data();
+  
+    for (Int_t j=0; j<3; ++j)
+    {
+      TString name;
+      name.Form("Multiplicity_%d", j);
+  
+      TFile::Open(files[j]);
+      data[j] = new AliMultiplicityCorrection(name, name);
+      data[j]->LoadHistograms("Multiplicity");
+    }
+    
+    // calculating relative
+    Float_t sd = data[0]->GetMultiplicityINEL(0)->Integral(0, data[0]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+    Float_t dd = data[1]->GetMultiplicityINEL(0)->Integral(0, data[1]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+    Float_t nd = data[2]->GetMultiplicityINEL(0)->Integral(0, data[2]->GetMultiplicityINEL(0)->GetNbinsX()+1);
+    Float_t total = nd + dd + sd;
+    
+    nd /= total;
+    sd /= total;
+    dd /= total;
+    
+    Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
+    
+    Float_t ratio[3];
+    ratio[0] = scalesSD / sd;
+    ratio[1] = scalesDD / dd;
+    ratio[2] = scalesND / nd;
+    
+    Printf("SD=%.2f, DD=%.2f, ND=%.2f", ratio[0], ratio[1], ratio[2]);
+    
+    TList list;
+    for (Int_t k=0; k<3; ++k)
+    {
+      // modify x-section
+      for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
+      {
+        data[k]->GetMultiplicityVtx(j)->Scale(ratio[k]);
+        data[k]->GetMultiplicityMB(j)->Scale(ratio[k]);
+        data[k]->GetMultiplicityINEL(j)->Scale(ratio[k]);
+        data[k]->GetMultiplicityNSD(j)->Scale(ratio[k]);
+      }
+  
+      if (k > 0)
+        list.Add(data[k]);
+    }
+  
+    printf("Case %d, %s: Entries in response matrix: SD: %.2f DD: %.2f ND: %.2f", xsectionID, data[0]->GetName(), data[0]->GetCorrelation(0)->Integral(), data[1]->GetCorrelation(0)->Integral(), data[2]->GetCorrelation(0)->Integral());
+  
+    data[0]->Merge(&list);
+    data[0]->FixTriggerEfficiencies(20);
+  
+    Printf(" Total: %.2f", data[0]->GetCorrelation(0)->Integral());
+  
+    for (Int_t etaR = 0; etaR < 3; etaR++)
+    {
+      unc[i+(etaR*9)] = (TH1*) data[0]->GetTriggerEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+      unc2[i+(etaR*9)] = (TH1*) data[0]->GetEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+    }
+  }
+      
+  file = TFile::Open(Form("systunc_fractions_%s.root", (eventType == 2) ? "inel" : "nsd"), "RECREATE");
+  
+  for (Int_t etaR = 0; etaR < 3; etaR++)
+  {
+    largestError1 = (TH1*) unc[0]->Clone(Form("fractions_trig_%d", etaR));
+    largestError1->Reset();
+    largestError2 = (TH1*) unc[0]->Clone(Form("fractions_trigvtx_%d", etaR));
+    largestError2->Reset();
+  
+    etaRange = etaR; // correct labels in DrawRatio
+    DrawRatio(unc[etaR*9],  8, unc+1+etaR*9,  Form("xsections_trig%d.png", etaR),    kFALSE, legendStrings, kFALSE, largestError1);
+    DrawRatio(unc2[etaR*9], 8, unc2+1+etaR*9, Form("xsections_trigvtx%d.png", etaR), kFALSE, legendStrings, kFALSE, largestError2);
+    
+    largestError2->SetBinContent(1, largestError1->GetBinContent(1));
+    
+    largestError2->GetYaxis()->UnZoom();
+    largestError2->Write();
+    
+    new TCanvas;
+    largestError2->DrawCopy();
+  }
+  file->Close();
+}
+
+void PythiaPhojetUncertainties()
+{
+  PythiaPhojetUncertainties(900, 2);
+  PythiaPhojetUncertainties(900, 3);
+  PythiaPhojetUncertainties(2360, 2);
+  PythiaPhojetUncertainties(2360, 3);
+}
+
+void PythiaPhojetUncertainties(Int_t energy, Int_t eventType)
+{
+  loadlibs();
+  
+  AliMultiplicityCorrection* mult[2];
+
+  const char* trigger = "all";
+  if (eventType == 3)
+    trigger = "v0and";
+    
+  if (energy == 2360)
+    trigger = "spdgfobits";
+
+  if (energy == 7000)
+    trigger = "all-onepart";
+  
+  if (energy == 900)
+    TFile::Open(Form("LHC10a12_run10482X/%s/spd/multiplicityMC_xsection.root", trigger));
+  else if (energy == 2360)
+    TFile::Open(Form("LHC10a10_run105054_7/%s/spd/multiplicityMC_xsection.root", trigger));
+  else
+    TFile::Open(Form("LHC10b2_run114931/%s/spd/multiplicity.root", trigger));
+  mult[0] = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  mult[0]->LoadHistograms("Multiplicity");
+  
+  if (energy == 900)
+    TFile::Open(Form("LHC10a14_run10482X_Phojet/%s/spd/multiplicityMC_xsection.root", trigger));
+  else if (energy == 2360)
+    TFile::Open(Form("LHC10a15_run10505X_Phojet/%s/spd/multiplicityMC_xsection.root", trigger));
+  else
+    TFile::Open(Form("LHC10b1_run114931/%s/spd/multiplicity.root", trigger));
+  mult[1] = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  mult[1]->LoadHistograms("Multiplicity");
+  
+  TH1* unc[6];
+  TH1* unc2[6];
+  
+  for (Int_t i=0; i<2; i++)
+  {
+    mult[i]->FixTriggerEfficiencies(20);
+    for (Int_t etaR = 0; etaR < 3; etaR++)
+    {
+      unc[i+(etaR*2)] = (TH1*) mult[i]->GetTriggerEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+      unc2[i+(etaR*2)] = (TH1*) mult[i]->GetEfficiency(etaR, eventType)->Clone(Form("eff_%d_%d", etaR, i));
+    }
+  }
+  
+  file = TFile::Open(Form("systunc_pythiaphojet_%d_%s.root", energy, (eventType == 2) ? "inel" : "nsd"), "RECREATE");
+  
+  for (Int_t etaR = 0; etaR < 3; etaR++)
+  {
+    largestError1Low = (TH1*) unc[0]->Clone(Form("pythiaphojet_trig_%d_low", etaR));
+    largestError1Low->Reset();
+    largestError1High = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trig_%d_high", etaR));
+    largestError2Low  = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trigvtx_%d_low", etaR));
+    largestError2High = (TH1*) largestError1Low->Clone(Form("pythiaphojet_trigvtx_%d_high", etaR));
+  
+    DrawRatio(unc[etaR*2], 1, unc+1+etaR*2, Form("pythiaphojet_trig%d.png", etaR), kFALSE, 0, kFALSE, largestError1Low, largestError1High);
+    DrawRatio(unc2[etaR*2], 1, unc2+1+etaR*2, Form("pythiaphojet_trigvtx%d.png", etaR), kFALSE, 0, kFALSE, largestError2Low, largestError2High);
+    
+    largestError2Low->SetBinContent(1, largestError1Low->GetBinContent(1));
+    largestError2High->SetBinContent(1, largestError1High->GetBinContent(1));
+    
+    largestError2Low->GetYaxis()->UnZoom();
+    largestError2High->GetYaxis()->UnZoom();
+    largestError2Low->Write();
+    largestError2High->Write();
+    
+    new TCanvas;
+    largestError2Low->DrawCopy();
+    
+    new TCanvas;
+    largestError2High->DrawCopy();
+  }
+  file->Close();
+}
+
+void CombinatoricsAsFunctionOfN(Int_t etaR)
+{
+  loadlibs();
+  
+  mult = AliMultiplicityCorrection::Open("multiplicityESD.root");
+  
+  //Float_t integratedCombinatorics[] = { 0.00062, 0.00074, 0.001 }; // 900 GeV
+  Float_t integratedCombinatorics[] = { 7.5e-4, 1.1e-3, 1.2e-3 };   // 2.36 TeV
+    
+  multHist = mult->GetMultiplicityESD(etaR)->ProjectionY();
+  multHist->Scale(1.0 / multHist->Integral());
+  
+  Float_t integral = 0;
+  Float_t tracks = 0;
+  for (Int_t i=1; i<= multHist->GetNbinsX(); i++)
+  {
+    Float_t x = multHist->GetXaxis()->GetBinCenter(i);
+    integral += x*x * multHist->GetBinContent(i+1);
+    tracks += x * multHist->GetBinContent(i+1);
+    Printf("%f %f %f", x, multHist->GetBinContent(i+1), integral);
+  }
+    
+  Printf("Integral is %f; %f", integral, integral / tracks);
+  Printf("alpha is %e", integratedCombinatorics[etaR] / (integral / tracks));
+}
+
+void TryCombinatoricsCorrection(Float_t etaR, Float_t alpha)
+{
+  loadlibs();
+  
+  mult = AliMultiplicityCorrection::Open("multiplicity.root");
+  
+  multHist = mult->GetMultiplicityESD(etaR)->ProjectionY();
+  
+  target = (TH1*) multHist->Clone("target");
+  target->Reset();
+  
+  Int_t totalTracks1 = 0;
+  Int_t totalTracks2 = 0;
+  
+  for (Int_t i=1; i<=multHist->GetNbinsX(); i++)
+  {
+    for (Int_t j=0; j<multHist->GetBinContent(i); j++)
+    {
+      Int_t origTracks = (Int_t) multHist->GetXaxis()->GetBinCenter(i);
+      tracks = origTracks;
+      //tracks -= gRandom->Poisson(1.3e-4 * origTracks * origTracks);
+      if (gRandom->Uniform() < alpha * origTracks * origTracks)
+        tracks--;
+      target->Fill(tracks);
+      
+      totalTracks1 += origTracks;
+      totalTracks2 += tracks;
+    }
+  }
+  
+  new TCanvas; multHist->Draw(); target->Draw("SAME"); target->SetLineColor(2);
+
+  Printf("%d %d %f %f", totalTracks1, totalTracks2, 1. - (Float_t) totalTracks2 / totalTracks1, 1. - target->GetMean() / multHist->GetMean());
+}
+
+void ApplyCombinatoricsCorrection()
+{
+  loadlibs();
+  
+  mult = AliMultiplicityCorrection::Open("multiplicity.root");
+  
+  fractionMC = new TF1("fractionMC", "[0]+[1]*x", 0, 200);
+  fractionMC->SetParameters(0.0004832, 5.82e-5); // mariella's presentation 16.04.10
+  
+  fractionData = new TF1("fractionData", "[0]+[1]*x", 0, 200);
+  fractionData->SetParameters(0.00052, 0.0001241); // mariella's presentation 16.04.10
+  
+  Int_t etaR = 1;
+  esd = mult->GetMultiplicityESD(etaR);
+  
+  for (Int_t x=1; x<=esd->GetNbinsX(); x++)
+  {
+    for (Int_t y=2; y<=esd->GetNbinsY(); y++)
+    {
+      printf("Before: %f  %f ", esd->GetBinContent(x, y-1), esd->GetBinContent(x, y));
+    
+      Float_t n = esd->GetYaxis()->GetBinCenter(y);
+      Int_t addComb = (fractionData->Eval(n) - fractionMC->Eval(n)) * n * esd->GetBinContent(x, y);
+      
+      // shift addComb down one bin
+      esd->SetBinContent(x, y-1, esd->GetBinContent(x, y-1) + addComb);
+      esd->SetBinContent(x, y,   esd->GetBinContent(x, y) - addComb);
+      
+      // recalc errors
+      esd->SetBinError(x, y-1, TMath::Sqrt(esd->GetBinContent(x, y-1)));
+      esd->SetBinError(x, y, TMath::Sqrt(esd->GetBinContent(x, y)));
+      
+      Printf("comb: %d; after: %f +- %f  %f +- %f", addComb, esd->GetBinContent(x, y-1), esd->GetBinError(x, y-1), esd->GetBinContent(x, y), esd->GetBinError(x, y));
+    }
+  }
+  
+  TFile::Open("out.root", "RECREATE");
+  mult->SaveHistograms();
+}
index 87f9b5a43749934194e344df83c702ca35640f23..096d75f5eb397053a2545f6498a3c304e0a2b007 100644 (file)
@@ -1,8 +1,13 @@
-void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug = kFALSE, Int_t aProof = 0, Bool_t mc = kTRUE, const char* option = "")
+void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug = kFALSE, Int_t aProof = 0, Int_t requiredData = 1, const char* option = "")
 {
   // aProof option: 0 no proof
   //                1 proof with chain
   //                2 proof with dataset
+  //
+  // requiredData option: 0 = only ESD
+  //                      1 = ESD+MC
+  //                      2 = RAW (ESD+check on event type)
+  //
 
   if (nRuns < 0)
     nRuns = 1234567890;
@@ -11,8 +16,7 @@ void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug =
   {
     gEnv->SetValue("XSec.GSI.DelegProxy", "2");
     TProof::Open("alicecaf");
-    //gProof->SetParallel(1);
-
+    
     // Enable the needed package
     if (1)
     {
@@ -51,10 +55,45 @@ void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug =
   // Create the analysis manager
   mgr = new AliAnalysisManager;
 
-  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn;
-  AliTriggerAnalysis::Trigger trigger      = AliTriggerAnalysis::kMB1;
+  // Add ESD handler
+  AliESDInputHandler* esdH = new AliESDInputHandler;
+  esdH->SetInactiveBranches("AliESDACORDE FMD ALIESDTZERO ALIESDZDC AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks AliESDTZERO ALIESDACORDE MuonTracks TrdTracks");
+  mgr->SetInputEventHandler(esdH);
+
+  // physics selection
+  gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  physicsSelectionTask = AddTaskPhysicsSelection((requiredData == 2) ? kFALSE : kTRUE);
+
+  // FO efficiency (for MC)
+  if (1 && requiredData != 2)
+  {
+    //const char* fastORFile = "../dNdEta/spdFOEff_run104824_52.root";
+    //const char* fastORFile = "../dNdEta/spdFOEff_run104867_92.root";
+    //const char* fastORFile = "../dNdEta/spdFOEff_run105054_7.root";
+    const char* fastORFile = "../dNdEta/spdFOEff_run114931.root";
+  
+    Printf("NOTE: Simulating FAST-OR efficiency on the analysis level using file %s", fastORFile);
+    TFile::Open(fastORFile);
+    
+    spdFOEff = (TH1F*) gFile->Get("spdFOEff");
+    physicsSelectionTask->GetPhysicsSelection()->Initialize(114931);
+    physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis()->SetSPDGFOEfficiency(spdFOEff);
+  }
+  
+  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn;
+  //AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPCITS | AliPWG0Helper::kFieldOn;
+  
+  //AliTriggerAnalysis::Trigger trigger      = AliTriggerAnalysis::kAcceptAll | AliTriggerAnalysis::kOfflineFlag;
+  AliTriggerAnalysis::Trigger trigger      = AliTriggerAnalysis::kAcceptAll | AliTriggerAnalysis::kOfflineFlag | AliTriggerAnalysis::kOneParticle; 
+  
+  //AliTriggerAnalysis::Trigger trigger      = AliTriggerAnalysis::kMB1Prime | AliTriggerAnalysis::kOfflineFlag;
+  //AliTriggerAnalysis::Trigger trigger      = AliTriggerAnalysis::kSPDGFOBits | AliTriggerAnalysis::kOfflineFlag;
+  //AliTriggerAnalysis::Trigger trigger      = AliTriggerAnalysis::kV0AND | AliTriggerAnalysis::kOfflineFlag; 
 
-  AliPWG0Helper::PrintConf(analysisMode, trigger);
+  AliPWG0Helper::DiffTreatment diffTreatment = AliPWG0Helper::kMCFlags;
+  //AliPWG0Helper::DiffTreatment diffTreatment = AliPWG0Helper::kE710Cuts;
+  
+  AliPWG0Helper::PrintConf(analysisMode, trigger, diffTreatment);
 
   TString taskName("AliMultiplicityTask.cxx+");
   if (aDebug)
@@ -66,7 +105,39 @@ void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug =
   } else
     gROOT->Macro(taskName);
 
-  task = new AliMultiplicityTask(option);
+  // 0 bin calculation
+  if (0)
+  {
+  }
+  
+  // V0 syst. study
+  if (0)
+  {
+    Printf("NOTE: Systematic study for VZERO enabled!");
+    //physicsSelectionTask->GetPhysicsSelection()->Initialize(104867);
+    for (Int_t i=0; i<1; i++)
+    {
+      // for MC and data
+      physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0HwPars(15, 61.5, 86.5);
+      physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0AdcThr(15);
+      // only for MC
+      //physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0HwPars(0, 0, 125);
+      //physicsSelectionTask->GetPhysicsSelection()->GetTriggerAnalysis(i)->SetV0AdcThr(0);
+    }
+  }
+
+  TString optionStr(option);
+  
+  // remove SAVE option if set
+  Bool_t save = kFALSE;
+  TString optionStr(option);
+  if (optionStr.Contains("SAVE"))
+  {
+    optionStr = optionStr(0,optionStr.Index("SAVE")) + optionStr(optionStr.Index("SAVE")+4, optionStr.Length());
+    save = kTRUE;
+  }
+  
+  task = new AliMultiplicityTask(optionStr);
 
   if (!(analysisMode & AliPWG0Helper::kSPD))
   {
@@ -81,22 +152,24 @@ void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug =
 
     task->SetTrackCuts(esdTrackCuts);
   }
-  else
-    task->SetDeltaPhiCut(0.05);
+  //else
+  //  task->SetDeltaPhiCut(0.05);
 
   task->SetAnalysisMode(analysisMode);
   task->SetTrigger(trigger);
+  task->SetDiffTreatment(diffTreatment);
 
-  if (mc)
+  if (requiredData == 1)
     task->SetReadMC();
 
   //task->SetUseMCVertex();
+  
+  if (requiredData != 2)
+    task->SetSkipParticles();
 
   mgr->AddTask(task);
 
-  TString optionStr(option);
-  
-  if (mc) {
+  if (requiredData == 1) {
     // Enable MC event handler
     AliMCEventHandler* handler = new AliMCEventHandler;
     if (!optionStr.Contains("particle-efficiency"))
@@ -119,18 +192,12 @@ void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug =
     task->SetPtSpectrum((TH1D*) hist->Clone("pt-spectrum"));
   }
 
-  // Add ESD handler
-  AliESDInputHandler* esdH = new AliESDInputHandler;
-  esdH->SetInactiveBranches("AliESDACORDE FMD ALIESDTZERO ALIESDZDC AliRawDataErrorLogs CaloClusters Cascades EMCALCells EMCALTrigger ESDfriend Kinks AliESDTZERO ALIESDACORDE MuonTracks TrdTracks");
-  mgr->SetInputEventHandler(esdH);
-
   // Attach input
   cInput  = mgr->GetCommonInputContainer();
   mgr->ConnectInput(task, 0, cInput);
 
   // Attach output
   cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
-  //cOutput->SetDataOwned(kTRUE);
   mgr->ConnectOutput(task, 0, cOutput);
 
   // Enable debug printouts
@@ -146,12 +213,61 @@ void run(Char_t* data, Long64_t nRuns = -1, Long64_t offset = 0, Bool_t aDebug =
     // process dataset
 
     mgr->StartAnalysis("proof", data, nRuns, offset);
+  
+    if (save)
+    {
+      TString path("maps/");
+      path += TString(data).Tokenize("/")->Last()->GetName();
+      
+      UInt_t triggerNoFlags = (UInt_t) trigger % (UInt_t) AliTriggerAnalysis::kStartOfFlags;
+      switch (triggerNoFlags)
+      {
+        case AliTriggerAnalysis::kAcceptAll: path += "/all"; break;
+        case AliTriggerAnalysis::kMB1: path += "/mb1"; break;
+        case AliTriggerAnalysis::kMB2: path += "/mb2"; break;
+        case AliTriggerAnalysis::kMB3: path += "/mb3"; break;
+        case AliTriggerAnalysis::kSPDGFO: path += "/spdgfo"; break;
+        case AliTriggerAnalysis::kSPDGFOBits: path += "/spdgfobits"; break;
+        case AliTriggerAnalysis::kV0AND: path += "/v0and"; break;
+        case AliTriggerAnalysis::kNSD1: path += "/nsd1"; break;
+        case AliTriggerAnalysis::kMB1Prime: path += "/mb1prime"; break;
+        default: Printf("ERROR: Trigger undefined for path to files"); return;
+      }
+      
+      if (trigger & AliTriggerAnalysis::kOneParticle)
+        path += "-onepart";
+      
+      if (analysisMode & AliPWG0Helper::kSPD)
+        path += "/spd";
+      
+      if (analysisMode & AliPWG0Helper::kTPC)
+        path += "/tpc";
+        
+      if (analysisMode & AliPWG0Helper::kTPCITS)
+        path += "/tpcits";
+
+      gSystem->mkdir(path, kTRUE);
+      
+      TString fileName("multiplicity");
+      if (optionStr.Contains("only-process-type-nd"))
+        fileName += "ND";
+      if (optionStr.Contains("only-process-type-sd"))
+        fileName += "SD";
+      if (optionStr.Contains("only-process-type-dd"))
+        fileName += "DD";
+      fileName += ".root";
+      
+      gSystem->Rename(fileName, path + "/" + fileName);
+      gSystem->Rename("event_stat.root", path + "/event_stat.root");
+      
+      Printf(">>>>> Moved files to %s", path.Data());
+    }  
   }
   else if (aProof == 3)
   {
     gROOT->ProcessLine(".L CreateChainFromDataSet.C");
     ds = gProof->GetDataSet(data)->GetStagedSubset();
-    chain = CreateChainFromDataSet(ds);
+    chain = CreateChainFromDataSet(ds, "esdTree", nRuns);
     mgr->StartAnalysis("local", chain, nRuns, offset);
   }
   else