changed histogram limits in multiplicity analysis
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 May 2009 07:36:17 +0000 (07:36 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 6 May 2009 07:36:17 +0000 (07:36 +0000)
new function to calculate multiplicity reach for high multiplicity trigger

PWG0/highMultiplicity/AliHighMultiplicitySelector.cxx
PWG0/highMultiplicity/AliHighMultiplicitySelector.h
PWG0/highMultiplicity/extrapolateXSection.C
PWG0/multiplicity/AliMultiplicityCorrection.cxx
PWG0/multiplicity/AliMultiplicityCorrection.h
PWG0/multiplicity/AliMultiplicityTask.cxx
PWG0/multiplicity/correct.C
PWG0/multiplicity/plots.C

index 49a1bfa..76543d0 100644 (file)
@@ -17,6 +17,8 @@
 #include <TLegend.h>
 #include <TLine.h>
 #include <TMath.h>
+#include <TLegend.h>
+#include <TText.h>
 
 #include <AliLog.h>
 #include <AliESD.h>
@@ -530,7 +532,7 @@ void AliHighMultiplicitySelector::ReadHistograms(const char* filename)
   }
 }
 
-TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut)
+TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t cut, Int_t upperCut)
 {
   //
   // returns the trigger efficiency as function of multiplicity with a given cut
@@ -541,7 +543,7 @@ TH1* AliHighMultiplicitySelector::GetTriggerEfficiency(TH2* multVsLayer, Int_t c
   //allEvents->Sumw2();
 
   //cut and multiply with x-section
-  TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, 1001);
+  TH1* proj = multVsLayer->ProjectionX(Form("%s_x", multVsLayer->GetName()), cut, upperCut);
   //proj->Sumw2();
 
   //new TCanvas; allEvents->DrawCopy(); gPad->SetLogy();
@@ -1310,6 +1312,7 @@ void AliHighMultiplicitySelector::Contamination2()
 void AliHighMultiplicitySelector::Contamination3()
 {
   //
+  // draws the contamination as function of treshold depending on a number a set of input MC and rate parameters
   //
 
   /*
@@ -1323,100 +1326,304 @@ void AliHighMultiplicitySelector::Contamination3()
 
   */
 
+  // output file
+  TFile* output = TFile::Open("contamination3.root", "RECREATE");
+
   // get x-sections
-  TFile* file = TFile::Open("crosssectionEx.root");
+  TFile* file = TFile::Open("crosssectionEx_10TeV.root");
   if (!file)
     return;
-
-  TH1* xSections[2];
-  xSections[0] = dynamic_cast<TH1*> (gFile->Get("xSection2Ex"));
-  xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
-
-  // prob for a collision in a bunch crossing
-  Int_t nRates = 3;
-  Float_t rates[] = {0.07, 0.1, 0.2};
-
-  new TCanvas;
-
-  Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
-  Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
-
-  // put to 2 for second layer
-  for (Int_t i=0; i<1; ++i)
+    
+  TCanvas* c = new TCanvas;
+  c->SetGridx();
+  c->SetGridy();
+  
+  TLegend* legend = new TLegend(0.7, 0.2, 1, 0.5);
+  legend->SetNColumns(2);
+  
+  TH2* dummy = new TH2F("dummy", ";Layer 1 Threshold;Contamination", 100, 95, 255, 100, 0, 1);
+  dummy->SetStats(kFALSE);
+  dummy->Draw();
+    
+  for (Int_t mc = 0; mc < 6; mc++)
   {
-    if (!xSections[i])
-      continue;
-
-    // relative x-section, once we have a collision
-    xSections[i]->Scale(1.0 / xSections[i]->Integral());
-
-    Int_t max = xSections[i]->GetNbinsX();
-    max = 500;
-
-    Float_t* xSection = new Float_t[max];
-    for (Int_t mult = 0; mult < max; mult++)
-      xSection[mult] = xSections[i]->GetBinContent(mult+1);
-
-    TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
-
-    for (Int_t currentRate = 0; currentRate<nRates; ++currentRate)
+    TH1* xSections[2];
+    TString str;
+    str.Form("xSection2Ex_%d_%d", mc/3, mc%3);
+    Printf("%s", str.Data());
+    file->cd();
+    xSections[0] = dynamic_cast<TH1*> (gFile->Get(str));
+    //xSections[1] = dynamic_cast<TH1*> (gFile->Get("xSection15Ex"));
+  
+    // prob for a collision in a bunch crossing
+    Int_t nRates = 1;
+    //Float_t rates[] = {0.02, 0.05, 0.1, 0.15, 0.2};
+    Float_t rates[] = {0.0636};
+    
+    // bunch crossing rate in Hz
+    Float_t bunchCrossingRate = 24. * 11245.5;
+  
+    Int_t colors[] = { 2, 3, 4, 6, 7, 8 };
+    Int_t markers[] = { 7, 2, 4, 5, 6, 27 };
+  
+    // put to 2 for second layer
+    for (Int_t i=0; i<1; ++i)
     {
-      TGraph* graph = new TGraph;
-      graph->SetMarkerColor(colors[currentRate]);
-      graph->SetMarkerStyle(markers[currentRate]);
-
-      Float_t rate = rates[currentRate];
-
-      for (Int_t cut = 100; cut <= 201; cut += 10)
+      if (!xSections[i])
+        continue;
+  
+      // relative x-section, once we have a collision
+      xSections[i]->Scale(1.0 / xSections[i]->Integral());
+  
+      Int_t max = xSections[i]->GetNbinsX();
+      max = 500;
+  
+      Float_t* xSection = new Float_t[max];
+      for (Int_t mult = 0; mult < max; mult++)
+        xSection[mult] = xSections[i]->GetBinContent(mult+1);
+  
+      TH2* fMvsL = (i == 0) ? fMvsL1: fMvsL2;
+  
+      for (Int_t currentRate = 0; currentRate<nRates; ++currentRate)
       {
-        Printf("cut at %d", cut);
+        TGraph* graph = new TGraph;
+        graph->SetMarkerColor(colors[currentRate]);
+        graph->SetMarkerStyle(markers[currentRate]);
+       TGraph* graph2 = new TGraph;
+        Float_t rate = rates[currentRate];
+  
+        Double_t singleRate = TMath::Poisson(1, rate);
+        Double_t doubleRate = TMath::Poisson(2, rate);
+        Double_t tripleRate = TMath::Poisson(3, rate);
+  
+        Printf("single = %f, double = %f, triple = %f", singleRate, doubleRate, tripleRate);
+    
+        for (Int_t cut = 100; cut <= 251; cut += 10)
+        {
+          Printf("Cut at %d", cut);
+  
+          TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+          Float_t* triggerEff = new Float_t[max];
+          for (Int_t mult = 0; mult < max; mult++)
+            triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
+    
+          Double_t triggerRate = 0;
+          for (Int_t mult = 0; mult < max; mult++)
+            triggerRate += xSection[mult] * triggerEff[mult];
+    
+          //Printf("  Raw value for 1  collision is %e; Rate: %.1f Hz", triggerRate, triggerRate * singleRate * bunchCrossingRate);
+    
+          Double_t triggerRate2 = 0;
+          for (Int_t mult = 0; mult < max; mult++)
+            for (Int_t mult2 = mult; mult2 < max; mult2++)
+              if (mult+mult2 < max)
+                triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
+    
+          //Printf("  Raw value for 2 collisions is %e; Rate: %.1f Hz", triggerRate2, triggerRate2 * doubleRate * bunchCrossingRate);
+    
+          Double_t triggerRate3 = 0;
+          for (Int_t mult = 0; mult < max; mult++)
+            for (Int_t mult2 = 0; mult2 < max; mult2++)
+              for (Int_t mult3 = 0; mult3 < max; mult3++)
+                if (mult+mult2+mult3 < max)
+                  triggerRate3 += xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3];
+          
+          //Printf("  Raw value for 3 collisions is %e; Rate: %.1f Hz", triggerRate3, triggerRate3 * tripleRate * bunchCrossingRate);
+          
+          Printf("  Rates: %.1f Hz; %.1f Hz; %.1f Hz", triggerRate * singleRate * bunchCrossingRate, triggerRate2 * doubleRate * bunchCrossingRate, triggerRate3 * tripleRate * bunchCrossingRate);
+          
+          Float_t totalTrigger = (triggerRate * singleRate + triggerRate2 * doubleRate + triggerRate3 * tripleRate);
+          
+          Printf("  Total trigger rate: %.1f Hz", totalTrigger * bunchCrossingRate);
+          
+          //if (totalTrigger * bunchCrossingRate > 200)
+          //  continue;
+    
+          Float_t totalContamination = (triggerRate2 * doubleRate + triggerRate3 * tripleRate) / totalTrigger;
+          //if (totalContamination > 0.99)
+         //  break;
+    
+          Printf("  Total contamination is %.1f%%", totalContamination * 100);
+    
+          graph->SetPoint(graph->GetN(), cut, totalContamination);
+         graph2->SetPoint(graph->GetN(), cut, totalTrigger * bunchCrossingRate);
+        }
+  
+        graph->SetMarkerStyle(mc+20);
+        graph->SetMarkerColor(currentRate+1);
+        graph->Draw("P SAME");
+        graph->GetXaxis()->SetTitle("Layer 1 threshold");
+        graph->GetYaxis()->SetTitle("Contamination");
+        graph->GetYaxis()->SetRangeUser(0, 1);
+        
+        if (currentRate == 0)
+        {
+          const char* legendLabel[] = { "Pythia Slope 1", "Pythia Slope 2", "Pythia Slope 3", "Phojet Slope 1", "Phojet Slope 2", "Phojet Slope 3" };
+          legend->AddEntry(graph, legendLabel[mc], "P");
+        }
 
-      TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
-      Float_t* triggerEff = new Float_t[max];
-      for (Int_t mult = 0; mult < max; mult++)
-        triggerEff[mult] = triggerEffHist->GetBinContent(mult+1);
+       output->cd();
+       graph->Write(Form("%s_%d_cont", str.Data(), currentRate));
+       graph2->Write(Form("%s_%d_rate", str.Data(), currentRate));
+      }
+    }
+  }
+  
+  output->Close();
+  
+  legend->Draw();
+}
 
-      Double_t triggerRate = 0;
-      for (Int_t mult = 0; mult < max; mult++)
-        triggerRate += xSection[mult] * triggerEff[mult];
+void AliHighMultiplicitySelector::Contamination_Reach()
+{
+  // plot the multiplicity reach based on the output from Contamination3()
+  // for each rate case and each MC, a certain number of events is required starting counting from the highest multiplicity
+  // note that the reach of different MC cannot be compared with each other
 
-      Printf("Raw value for 1 collision is %e", triggerRate);
+  /*
 
-      Double_t triggerRate2 = 0;
-      for (Int_t mult = 0; mult < max; mult++)
-        for (Int_t mult2 = mult; mult2 < max; mult2++)
-          if (mult+mult2 < max)
-            triggerRate2 += ((mult2 > mult) ? 2. : 1.) * xSection[mult] * xSection[mult2] * triggerEff[mult+mult2];
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libPWG0base");
+  .L AliHighMultiplicitySelector.cxx+g
+  x = new AliHighMultiplicitySelector();
+  x->ReadHistograms("highmult_hijing100k.root");
+  x->Contamination_Reach();
 
-      Printf("Raw value for 2 collisions is %e", triggerRate2);
+  */
 
-      Double_t triggerRate3 = 0;
-      for (Int_t mult = 0; mult < max; mult++)
-        for (Int_t mult2 = 0; mult2 < max; mult2++)
-          for (Int_t mult3 = 0; mult3 < max; mult3++)
-            if (mult+mult2+mult3 < max)
-              triggerRate3 += xSection[mult] * xSection[mult2] * xSection[mult3] * triggerEff[mult+mult2+mult3];
+  TCanvas* c = new TCanvas("c", "c", 800, 600);
+  c->Divide(2, 3);
+  
+  // prob for a collision in a bunch crossing
+  Int_t nRates = 1;
+  //Float_t rates[] = {0.02, 0.05, 0.1, 0.15, 0.2};
+  Float_t rates[] = {0.0636};
 
-      Printf("Raw value for 3 collisions is %e", triggerRate3);
+  // bunch crossing rate in Hz
+  Float_t bunchCrossingRate = 24. * 11245.5;
 
-      Double_t singleRate = TMath::Poisson(1, rate);
-      Double_t doubleRate = TMath::Poisson(2, rate);
-      Double_t tripleRate = TMath::Poisson(3, rate);
+  TH2* dummy = new TH2F("dummy", ";Coll/bunch crossing;multiplicity reach", 100, 0, 0.3, 100, 50, 350);
+  //TH2* dummy = new TH2F("dummy", ";Coll/bunch crossing;fractional cross-section", 100, 0, 0.3, 1000, 1e-6, 0.1);
+  dummy->SetStats(kFALSE);
 
-      Printf("single = %f, double = %f, triple = %f", singleRate, doubleRate, tripleRate);
+  const char* legendLabel[] = { "Pythia Slope 1", "Pythia Slope 2", "Pythia Slope 3", "Phojet Slope 1", "Phojet Slope 2", "Phojet Slope 3" };
 
-      Float_t totalContamination = (triggerRate2 * doubleRate + triggerRate3 * tripleRate) / (triggerRate * singleRate + triggerRate2 * doubleRate + triggerRate3 * tripleRate);
+  TFile* mcFile = TFile::Open("crosssectionEx_10TeV.root");
+  TFile* contFile = TFile::Open("contamination3.root");
 
-        //Printf("Total contamination is %.1f%%", totalContamination * 100);
+  // for comparison: how many MB events can one take at the same time
+  Int_t mbEvents = 2e6 * 500;
 
-      graph->SetPoint(graph->GetN(), cut, totalContamination);
+  for (Int_t mc=0; mc<6; mc++)
+  {
+    mcFile->cd();
+    TH1* mcHist = (TH1*) gFile->Get(Form("xSection2Ex_%d_%d", mc/3, mc%3));
+    mcHist->Scale(1.0 / mcHist->Integral());
+    
+    c->cd(mc+1);//->SetLogy();
+    c->SetGridx();
+    c->SetGridy();
+    dummy->Draw();
+    
+    Int_t color = 0;
+    for (Int_t requiredEvents = 300; requiredEvents <= 3000000; requiredEvents *= 10)
+    {
+      TGraph* reach = new TGraph;
+
+      color++;
+      if (color == 5)
+        color++;
+      
+      Float_t requiredRate = (Float_t) requiredEvents / 1e6;
+      Printf("Required rate is %f", requiredRate);
+    
+      // find reach without trigger
+      Int_t mbReach = 1000;
+      while (mcHist->Integral(mcHist->FindBin(mbReach), mcHist->GetNbinsX()) < (Float_t) requiredEvents / mbEvents && mbReach > 1)
+        mbReach--;
+      Printf("MB reach is %d with %f events", mbReach, mcHist->Integral(mcHist->FindBin(mbReach), mcHist->GetNbinsX()) * mbEvents);
+      
+      for (Int_t rate=0; rate<nRates; rate++)
+      {
+        contFile->cd();
+        TGraph* cont = (TGraph*) gFile->Get(Form("xSection2Ex_%d_%d_%d_cont", mc/3, mc%3, rate));
+        TGraph* rateh = (TGraph*) gFile->Get(Form("xSection2Ex_%d_%d_%d_rate", mc/3, mc%3, rate));
+  
+        Double_t singleRate = TMath::Poisson(1, rates[rate]);
+        Double_t totalCollRate = singleRate * bunchCrossingRate;
+        Printf("collisions/bc: %f; coll. rate: %f", singleRate, totalCollRate);
+  
+        // find 200 Hz limit
+        Int_t low = 100;
+        while (rateh->Eval(low) > 200)
+          low++;
+  
+        // find contamination limit
+        Int_t high = 100;
+        while (cont->Eval(high) < 0.9 && high < 250)
+          high++;
+  
+        Printf("MC %d Rate %d: Acceptable threshold region is %d to %d", mc, rate, low, high);
+        // find reachable multiplicity; include contamination in rate calculation
+        // TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL, cut)->Clone("triggerEff");
+  
+        // trigger efficiency as function of multiplicity in range mult <= ... <= high
+        //new TCanvas; fMvsL1->Draw("COLZ");
+        TH1* triggerEffHist = (TH1*) GetTriggerEfficiency(fMvsL1, low, high);
+        
+        //new TCanvas; triggerEffHist->DrawCopy();
+        triggerEffHist->Multiply(mcHist);
+        
+        Float_t fractionXSection = triggerEffHist->Integral();
+        Printf("The fraction of the cross-section is %f", fractionXSection);
+        
+        //new TCanvas; triggerEffHist->DrawCopy();
+        triggerEffHist->Scale(totalCollRate);
+        //new TCanvas; triggerEffHist->DrawCopy(); gPad->SetLogy();
+        
+        Float_t achievedRate = 0;
+        Int_t mult = 1000;
+        while (1)
+        {
+          achievedRate = triggerEffHist->Integral(triggerEffHist->FindBin(mult), triggerEffHist->GetNbinsX());
+  
+          if (achievedRate >= requiredRate)
+            break;
+      
+          if (mult == 1)
+            break;
+      
+          mult--;
+        } 
+  
+        Printf("Achieved rate %f above multiplicity %d", achievedRate, mult);
+        
+        if (achievedRate < requiredRate)
+        {
+          Printf("Achieved rate too low");
+          continue;
+        }
+  
+        reach->SetPoint(reach->GetN(), rates[rate], mult);
+        //reach->SetPoint(reach->GetN(), rates[rate], fractionXSection);
+        
+      
+        //return;
       }
-
-      graph->Draw((currentRate == 0) ? "A*" : "* SAME");
-      graph->GetXaxis()->SetTitle("layer 1 threshold");
-      graph->GetYaxis()->SetTitle("contamination");
-      graph->GetYaxis()->SetRangeUser(0, 1);
-    }
+  
+      reach->SetMarkerColor(color);
+      reach->Draw("SAME*");
+      
+      TLine* line = new TLine(0, mbReach, 0.3, mbReach);
+      line->SetLineColor(color);
+      //line->Draw();
+    }    
+  
+    TText* text = new TText;
+    text->DrawText(0.2, 325, legendLabel[mc]);
+    //return;
   }
 }
 
index 110d857..118b75f 100644 (file)
@@ -38,6 +38,7 @@ class AliHighMultiplicitySelector : public AliSelectorRL {
     void Contamination();
     void Contamination2();
     void Contamination3();
+    void Contamination_Reach();
     void MBComparison();
 
  protected:
@@ -45,7 +46,7 @@ class AliHighMultiplicitySelector : public AliSelectorRL {
     void MakeGraphs2(const char* title, TH1* xSection, TH2* fMvsL);
 
     TH1* GetXSectionCut(TH1* xSection, TH2* multVsLayer, Int_t cut);
-    TH1* GetTriggerEfficiency(TH2* multVsLayer, Int_t cut);
+    TH1* GetTriggerEfficiency(TH2* multVsLayer, Int_t cut, Int_t upperCut = 1001);
 
     TH1F* fChipsLayer1;   // fired chips in layer 1
     TH1F* fChipsLayer2;   // fired chips in layer 2
index 312b772..39219a4 100644 (file)
@@ -1,63 +1,80 @@
 void extrapolateXSection()
 {
-  TFile* file = TFile::Open("crosssection.root");
-  TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("xSection2"));
-  TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15"));
-
-  gROOT->cd();
-
-  TH1F* xSection2Ex = new TH1F("xSection2Ex", ";Npart", 1001, -0.5, 1000.5);
-  TH1F* xSection15Ex = new TH1F("xSection15Ex", ";Npart", 1001, -0.5, 1000.5);
-
-  new TCanvas;
-  xSection2->Draw();
-  xSection2->Fit("expo", "", "", 200, 250);
-  gPad->SetLogy();
-
-  for (Int_t i=1; i<=1000; ++i)
+  TFile* file2 = TFile::Open("crosssectionEx_10TeV.root", "RECREATE");
+    
+  c = new TCanvas;
+  TH1* base = 0;
+  
+  for (Int_t fileId = 0; fileId < 2; fileId++)
   {
-    if (i < 250)
+    if (fileId == 1)
     {
-      xSection2Ex->SetBinContent(i, xSection2->GetBinContent(i));
-      xSection2Ex->SetBinError(i, xSection2->GetBinError(i));
+      TFile* file = TFile::Open("out_phojet.root");
     }
     else
-      xSection2Ex->SetBinContent(i, xSection2->GetFunction("expo")->Eval(i));
-  }
-
-  new TCanvas;
-  xSection2Ex->Draw();
-  xSection2->SetLineColor(2);
-  xSection2->Draw("SAME");
-  gPad->SetLogy();
-
-  new TCanvas;
-  xSection15->Draw();
-  xSection15->Fit("expo", "", "", 145, 250);
-  gPad->SetLogy();
-
-  for (Int_t i=1; i<=1000; ++i)
-  {
-    if (i < 145)
+      TFile* file = TFile::Open("out_pythia.root");
+    
+    TH1* xSection2 = dynamic_cast<TH1*> (gFile->Get("fMult"));
+    if (!base)
+      base = xSection2;
+    xSection2->Sumw2();
+    xSection2->Scale(1.0 / xSection2->Integral());
+    //TH1* xSection15 = dynamic_cast<TH1*> (gFile->Get("xSection15"));
+  
+    //TH1F* xSection15Ex = new TH1F("xSection15Ex", ";Npart", 1001, -0.5, 1000.5);
+  
+    TF1* func[3];
+    Float_t lowLimit[] = { 150, 175, 200 };
+    if (fileId == 1)
     {
-      xSection15Ex->SetBinContent(i,xSection15->GetBinContent(i));
-      xSection15Ex->SetBinError(i, xSection15->GetBinError(i));
+      lowLimit[0] = 50;
+      lowLimit[1] = 75;
+      lowLimit[2] = 100;
+    }
+    
+    c->cd();
+    xSection2->Draw((fileId == 0) ? "" : "SAME");
+    
+    for (Int_t i=0; i<3; i++)
+    {
+      func[i] = new TF1("func", "[0]*exp([1]*x)", 0, 1000);
+      func[i]->SetParameters(1, -1e-4);
+    
+      xSection2->Fit(func[i], "0", "", lowLimit[i], 250);
+      func[i]->SetRange(lowLimit[i], 500);
+      func[i]->SetLineColor(i+1);
+      func[i]->Draw("SAME");
+      gPad->SetLogy();
+    }
+    
+    base->GetXaxis()->SetRangeUser(0, 500);
+    base->GetYaxis()->SetRangeUser(func[2]->Eval(500), xSection2->GetMaximum());
+    
+    for (Int_t j=0; j<3; j++)
+    {
+      gROOT->cd();
+      TH1F* xSection2Ex = new TH1F(Form("xSection2Ex_%d_%d", fileId, j), ";Npart", 1001, -0.5, 1000.5);
+      
+      for (Int_t i=1; i<=1000; ++i)
+      {
+        if (i < lowLimit[j])
+        {
+          xSection2Ex->SetBinContent(i, xSection2->GetBinContent(i));
+          xSection2Ex->SetBinError(i, xSection2->GetBinError(i));
+        }
+        else
+          xSection2Ex->SetBinContent(i, func[j]->Eval(i));
+      }
+      
+      new TCanvas;
+      xSection2Ex->Draw();
+      gPad->SetLogy();
+      
+      file2->cd();
+      xSection2Ex->Write();
     }
-    else
-      xSection15Ex->SetBinContent(i, xSection15->GetFunction("expo")->Eval(i));
   }
-
-  new TCanvas;
-  xSection15Ex->Draw();
-  xSection15->SetLineColor(2);
-  xSection15->Draw("SAME");
-  gPad->SetLogy();
-
-  TFile* file2 = TFile::Open("crosssectionEx.root", "RECREATE");
-
-  xSection2Ex->Write();
-  xSection15Ex->Write();
-
+  
   file2->Close();
 }
 
index 45d4217..187ba41 100644 (file)
@@ -175,29 +175,30 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C
   #define NBINNING fgkMaxParams, binLimitsN*/
 
   #define NBINNING 201, -0.5, 200.5
-  #define VTXBINNING 1, -6, 6
+  
+  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", VTXBINNING, NBINNING);
+    fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 1, -vtxRange[i], vtxRange[i], NBINNING);
 
   for (Int_t i = 0; i < kMCHists; ++i)
   {
-    fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
+    fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[i%3]->Clone(Form("fMultiplicityVtx%d", i)));
     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
 
-    fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
+    fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityMB%d", i)));
     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
 
-    fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
+    fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityINEL%d", i)));
     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
     
-    fMultiplicityNSD[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityNSD%d", i)));
+    fMultiplicityNSD[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityNSD%d", i)));
     fMultiplicityNSD[i]->SetTitle("fMultiplicityNSD");
   }
 
   for (Int_t i = 0; i < kCorrHists; ++i)
   {
-    fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
+    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);
   }
 
@@ -455,7 +456,7 @@ void AliMultiplicityCorrection::SaveHistograms(const char* dir)
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, AliPWG0Helper::MCProcessType processType, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
+void AliMultiplicityCorrection::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)
 {
   //
   // Fills an event from MC
@@ -465,38 +466,34 @@ void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Boo
   {
     fMultiplicityMB[0]->Fill(vtx, generated05);
     fMultiplicityMB[1]->Fill(vtx, generated10);
-    fMultiplicityMB[2]->Fill(vtx, generated15);
-    fMultiplicityMB[3]->Fill(vtx, generated20);
-    fMultiplicityMB[4]->Fill(vtx, generatedAll);
+    fMultiplicityMB[2]->Fill(vtx, generated14);
+    fMultiplicityMB[3]->Fill(vtx, generatedAll);
 
     if (vertex)
     {
       fMultiplicityVtx[0]->Fill(vtx, generated05);
       fMultiplicityVtx[1]->Fill(vtx, generated10);
-      fMultiplicityVtx[2]->Fill(vtx, generated15);
-      fMultiplicityVtx[3]->Fill(vtx, generated20);
-      fMultiplicityVtx[4]->Fill(vtx, generatedAll);
+      fMultiplicityVtx[2]->Fill(vtx, generated14);
+      fMultiplicityVtx[3]->Fill(vtx, generatedAll);
     }
   }
 
   fMultiplicityINEL[0]->Fill(vtx, generated05);
   fMultiplicityINEL[1]->Fill(vtx, generated10);
-  fMultiplicityINEL[2]->Fill(vtx, generated15);
-  fMultiplicityINEL[3]->Fill(vtx, generated20);
-  fMultiplicityINEL[4]->Fill(vtx, generatedAll);
+  fMultiplicityINEL[2]->Fill(vtx, generated14);
+  fMultiplicityINEL[3]->Fill(vtx, generatedAll);
   
   if (processType != AliPWG0Helper::kSD)
   {
     fMultiplicityNSD[0]->Fill(vtx, generated05);
     fMultiplicityNSD[1]->Fill(vtx, generated10);
-    fMultiplicityNSD[2]->Fill(vtx, generated15);
-    fMultiplicityNSD[3]->Fill(vtx, generated20);
-    fMultiplicityNSD[4]->Fill(vtx, generatedAll);
+    fMultiplicityNSD[2]->Fill(vtx, generated14);
+    fMultiplicityNSD[3]->Fill(vtx, generatedAll);
   }
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
+void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured14)
 {
   //
   // Fills an event from ESD
@@ -504,12 +501,11 @@ void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_
 
   fMultiplicityESD[0]->Fill(vtx, measured05);
   fMultiplicityESD[1]->Fill(vtx, measured10);
-  fMultiplicityESD[2]->Fill(vtx, measured15);
-  fMultiplicityESD[3]->Fill(vtx, measured20);
+  fMultiplicityESD[2]->Fill(vtx, measured14);
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
+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)
 {
   //
   // Fills an event into the correlation map with the information from MC and ESD
@@ -517,13 +513,11 @@ void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, I
 
   fCorrelation[0]->Fill(vtx, generated05, measured05);
   fCorrelation[1]->Fill(vtx, generated10, measured10);
-  fCorrelation[2]->Fill(vtx, generated15, measured15);
-  fCorrelation[3]->Fill(vtx, generated20, measured20);
+  fCorrelation[2]->Fill(vtx, generated14, measured14);
 
-  fCorrelation[4]->Fill(vtx, generatedAll, measured05);
-  fCorrelation[5]->Fill(vtx, generatedAll, measured10);
-  fCorrelation[6]->Fill(vtx, generatedAll, measured15);
-  fCorrelation[7]->Fill(vtx, generatedAll, measured20);
+  fCorrelation[3]->Fill(vtx, generatedAll, measured05);
+  fCorrelation[4]->Fill(vtx, generatedAll, measured10);
+  fCorrelation[5]->Fill(vtx, generatedAll, measured14);
 }
 
 //____________________________________________________________________
index 63204d3..0ff4a34 100644 (file)
@@ -32,7 +32,7 @@ class TCollection;
 class AliMultiplicityCorrection : public TNamed {
   public:
     enum EventType { kTrVtx = 0, kMB, kINEL, kNSD };
-    enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8, kQualityRegions = 3 };
+    enum { kESDHists = 3, kMCHists = 4, kCorrHists = 6, kQualityRegions = 3 };
 
     AliMultiplicityCorrection();
     AliMultiplicityCorrection(const Char_t* name, const Char_t* title);
@@ -42,10 +42,10 @@ class AliMultiplicityCorrection : public TNamed {
 
     virtual Long64_t Merge(TCollection* list);
 
-    void FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20);
-    void FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, AliPWG0Helper::MCProcessType processType, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll);
+    void FillMeasured(Float_t vtx, 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 generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20);
+    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);
 
     Bool_t LoadHistograms(const Char_t* dir = 0);
     void SaveHistograms(const char* dir = 0);
@@ -104,14 +104,14 @@ class AliMultiplicityCorrection : public TNamed {
     TH2* fCurrentCorrelation; //! current correlation
     TH1* fCurrentEfficiency;  //! current efficiency
 
-    TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2 (0..3)
+    TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; 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.5, 2, inf (0..4)
-    TH2F* fMultiplicityMB[kMCHists];   // multiplicity histogram of triggered events                        : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2, inf (0..4)
-    TH2F* fMultiplicityINEL[kMCHists]; // multiplicity histogram of all (inelastic) events                  : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2, inf (0..4)
-    TH2F* fMultiplicityNSD[kMCHists]; // multiplicity histogram of NSD events                  : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.5, 2, inf (0..4)
+    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)
+    TH2F* fMultiplicityINEL[kMCHists]; // multiplicity histogram of all (inelastic) events                  : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4, inf (0..3)
+    TH2F* fMultiplicityNSD[kMCHists]; // multiplicity histogram of NSD events                  : vtx vs multiplicity; array: |eta| < 0.5, 1.0, 1.4, inf (0..3)
 
-    TH3F* fCorrelation[kCorrHists];              // vtx vs. (gene multiplicity (trig+vtx)) vs. (meas multiplicity); array: |eta| < 0.5, 1, 1.5, 2 (0..3 and 4..7), the first corrects to the eta range itself, the second to full phase space
+    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
 
@@ -129,7 +129,7 @@ class AliMultiplicityCorrection : public TNamed {
     AliMultiplicityCorrection(const AliMultiplicityCorrection&);
     AliMultiplicityCorrection& operator=(const AliMultiplicityCorrection&);
 
-  ClassDef(AliMultiplicityCorrection, 4);
+  ClassDef(AliMultiplicityCorrection, 5);
 };
 
 #endif
index 621b30c..7cf4de0 100644 (file)
@@ -318,8 +318,7 @@ void AliMultiplicityTask::Exec(Option_t*)
     {
       Int_t nESDTracks05 = 0;
       Int_t nESDTracks10 = 0;
-      Int_t nESDTracks15 = 0;
-      Int_t nESDTracks20 = 0;
+      Int_t nESDTracks14 = 0;
 
       for (Int_t i=0; i<inputCount; ++i)
       {
@@ -331,14 +330,11 @@ void AliMultiplicityTask::Exec(Option_t*)
         if (TMath::Abs(eta) < 1.0)
           nESDTracks10++;
 
-        if (TMath::Abs(eta) < 1.5)
-          nESDTracks15++;
-
-        if (TMath::Abs(eta) < 2.0)
-          nESDTracks20++;
+        if (TMath::Abs(eta) < 1.4)
+          nESDTracks14++;
       }
 
-      fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
+      fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
     }
   }
   else if (fReadMC)   // Processing of MC information
@@ -419,8 +415,7 @@ void AliMultiplicityTask::Exec(Option_t*)
       // tracks in different eta ranges
       Int_t nMCTracks05 = 0;
       Int_t nMCTracks10 = 0;
-      Int_t nMCTracks15 = 0;
-      Int_t nMCTracks20 = 0;
+      Int_t nMCTracks14 = 0;
       Int_t nMCTracksAll = 0;
 
       // tracks per particle species, in |eta| < 2 (systematic study)
@@ -489,11 +484,8 @@ void AliMultiplicityTask::Exec(Option_t*)
         if (TMath::Abs(particle->Eta()) < 1.0)
           nMCTracks10 += particleWeight;
 
-        if (TMath::Abs(particle->Eta()) < 1.5)
-          nMCTracks15 += particleWeight;
-
-        if (TMath::Abs(particle->Eta()) < 2.0)
-          nMCTracks20 += particleWeight;
+        if (TMath::Abs(particle->Eta()) < 1.4)
+          nMCTracks14 += particleWeight;
 
         nMCTracksAll += particleWeight;
         
@@ -519,14 +511,13 @@ void AliMultiplicityTask::Exec(Option_t*)
         }
       } // end of mc particle
 
-      fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
+      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 nESDTracks15 = 0;
-        Int_t nESDTracks20 = 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)
@@ -600,12 +591,8 @@ void AliMultiplicityTask::Exec(Option_t*)
           if (TMath::Abs(eta) < 1.0)
             nESDTracks10 += particleWeight;
 
-          if (TMath::Abs(eta) < 1.5)
-            nESDTracks15 += particleWeight;
-
-          if (TMath::Abs(eta) < 2.0)
-            nESDTracks20 += particleWeight;
-
+          if (TMath::Abs(eta) < 1.4)
+            nESDTracks14 += particleWeight;
 
           if (fParticleSpecies)
           {
@@ -788,16 +775,16 @@ void AliMultiplicityTask::Exec(Option_t*)
         delete[] foundPrimaries;
         delete[] foundPrimaries2;
 
-        if ((Int_t) nMCTracks15 > 10 && nESDTracks15 <= 3)
+        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(), nMCTracks15, nESDTracks15);
+            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);
         }
 
         // fill response matrix using vtxMC (best guess)
-        fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks15,  nMCTracks20,  nMCTracksAll,  nESDTracks05,  nESDTracks10,  nESDTracks15,  nESDTracks20);
+        fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks14,  nMCTracksAll,  nESDTracks05,  nESDTracks10, nESDTracks14);
 
-        fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
+        fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 
         if (fParticleSpecies)
           fParticleSpecies->Fill(vtxMC[2], nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
index 25fd74f..df19334 100644 (file)
@@ -66,47 +66,51 @@ void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder
 
   AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
   AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
+  
+  AliUnfolding::SetNbins(150, 150);
 
-  TH2F* hist = esd->GetMultiplicityESD(histID);
-  TH2F* hist2 = esd->GetMultiplicityMC(histID, eventType);
-  TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
-
-  mult->SetMultiplicityESD(histID, hist);
-
-  // small hack to get around charge conservation for full phase space ;-)
-  if (fullPhaseSpace)
-  {
-    TH1* corr = mult->GetCorrelation(histID + 4);
-
-    for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
-      for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
-      {
-        corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
-        corr->SetBinError(i, j, corr->GetBinError(i-1, j));
-      }
-  }
-
-  if (chi2)
+  for (hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
   {
-    AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, beta);
-    //mult->SetCreateBigBin(kFALSE);
-    //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
-    //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
-    //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e4);
-    //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
-    
-    //mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kTRUE, mcCompare);
-    //mult->SetMultiplicityESDCorrected(histID, (TH1F*) mcCompare);
-    
-    mult->ApplyMinuitFit(histID, fullPhaseSpace, eventType, kFALSE); //hist2->ProjectionY("mymchist"));
-    //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType);
-  }
-  else
-  {
-    mult->ApplyBayesianMethod(histID, fullPhaseSpace, eventType, 1, 10, 0, kTRUE);
+    TH2F* hist = esd->GetMultiplicityESD(hID);
+    TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
+  
+    mult->SetMultiplicityESD(histID, hist);
+  
+    // small hack to get around charge conservation for full phase space ;-)
+    if (fullPhaseSpace)
+    {
+      TH1* corr = mult->GetCorrelation(histID + 4);
+  
+      for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
+        for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
+        {
+          corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
+          corr->SetBinError(i, j, corr->GetBinError(i-1, j));
+        }
+    }
+  
+    if (chi2)
+    {
+      AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, beta);
+      //mult->SetCreateBigBin(kFALSE);
+      //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
+      //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
+      //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e4);
+      //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
+      
+      //mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kTRUE, mcCompare);
+      //mult->SetMultiplicityESDCorrected(histID, (TH1F*) mcCompare);
+      
+      mult->ApplyMinuitFit(hID, fullPhaseSpace, eventType, kFALSE); //hist2->ProjectionY("mymchist"));
+      //mult->ApplyNBDFit(histID, fullPhaseSpace, eventType);
+    }
+    else
+    {
+      mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, 10, 0, kTRUE);
+    }
+  
+    mult->SetMultiplicityMC(hID, eventType, hist2);
   }
-
-  mult->SetMultiplicityMC(histID, eventType, hist2);
   
   Printf("Writing result to %s", targetFile);
   TFile* file = TFile::Open(targetFile, "RECREATE");
@@ -114,7 +118,12 @@ void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder
   file->Write();
   file->Close();
 
-  mult->DrawComparison((chi2) ? "MinuitChi2" : "Bayesian", histID, fullPhaseSpace, kTRUE, mcCompare);
+  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);
+  }
 }
 
 void DrawBayesianIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 1, Int_t eventType = 0 /* AliMultiplicityCorrection::kTrVtx */)
index 1005a0d..3927aa6 100644 (file)
@@ -929,7 +929,9 @@ void DrawResiduals(const char* fileName, const char* epsName)
 
   // projection
   TH1* residualHist = new TH1F("residualHist", ";", 11, -3, 3);
+  residualHist->Sumw2();
 
+  Float_t chi2 = 0;
   for (Int_t i=1; i<=displayRange+1; ++i)
   {
     if (measured->GetBinError(i) > 0)
@@ -938,6 +940,7 @@ void DrawResiduals(const char* fileName, const char* epsName)
       residual->SetBinError(i, 1);
 
       residualHist->Fill(residual->GetBinContent(i));
+      chi2 += residual->GetBinContent(i) * residual->GetBinContent(i);
     }
     else
     {
@@ -945,6 +948,8 @@ void DrawResiduals(const char* fileName, const char* epsName)
       residual->SetBinError(i, 0);
     }
   }
+  
+  Printf("chi2 / ndf = %f / %d = %f", chi2, displayRange+1, chi2 / (displayRange+1));
 
   residual->GetYaxis()->SetTitle("Residuals:   (1/e) (M - R  #otimes U)");
   residual->GetYaxis()->SetRangeUser(-4.5, 4.5);
@@ -958,7 +963,8 @@ void DrawResiduals(const char* fileName, const char* epsName)
   residualHist->SetStats(kFALSE);
   residualHist->SetLabelSize(0.08, "xy");
   residualHist->Fit("gaus");
-  residualHist->Draw();
+  residualHist->Draw("HIST");
+  residualHist->FindObject("gaus")->Draw("SAME");
 
   canvas->Modified();
   canvas->SaveAs(canvas->GetName());
@@ -2058,7 +2064,7 @@ void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
 
   loadlibs();
 
-  TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 500);
+  TCanvas* canvas = new TCanvas("EfficiencyComparison", "EfficiencyComparison", 800, 600);
   canvas->SetGridx();
   canvas->SetGridy();
   canvas->SetRightMargin(0.05);
@@ -2914,7 +2920,7 @@ void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE)
 
   //TH1* errorNSD = SystematicsSummary(tpc, 1);
 
-  TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 500);
+  TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 500); 
   canvas->SetRightMargin(0.05);
   canvas->SetTopMargin(0.05);
   canvas->SetGridx();
@@ -3026,6 +3032,122 @@ void finalPlot(Bool_t tpc = 0, Bool_t small = kFALSE)
   canvas->SaveAs(canvas->GetName());
 }
 
+void finalPlot2(Bool_t tpc = 0)
+{
+  loadlibs();
+
+  if (tpc)
+    SetTPC();
+
+  //TH1* errorNSD = SystematicsSummary(tpc, 1);
+
+  TCanvas* canvas = new TCanvas("finalPlot2.eps", "finalPlot2.eps", 800, 600); 
+  canvas->SetRightMargin(0.05);
+  canvas->SetTopMargin(0.05);
+  canvas->SetGridx();
+  canvas->SetGridy();
+  canvas->SetLogy();
+  
+  legend = new TLegend(0.5, 0.7, 0.9, 0.9);
+  legend->SetFillColor(0);
+  legend->SetTextSize(0.04);
+  
+  Int_t displayRanges[] = { 50, 80, 120 };
+  
+  TH2* dummy = new TH2F("dummy", ";True multiplicity N_{ch};P(N_{ch})", 100, -0.5, displayRanges[2]+10, 1000, 5e-5, 5);
+  dummy->SetStats(0);
+  dummy->Draw();
+  
+  TList list;
+  
+  Int_t count = 0;
+  for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kINEL; eventType <= AliMultiplicityCorrection::kNSD; eventType++)
+  {
+    AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open((eventType == AliMultiplicityCorrection::kINEL) ? "chi2_inel.root" : "chi2_nsd.root");
+    //TH1* mcHist = mult->GetMultiplicityMC(etaRange, eventType)->ProjectionY("mymc");
+    for (Int_t etaR = 0; etaR <= 2; etaR++)
+    {
+      TH1* result = mult->GetMultiplicityESDCorrected(etaR);
+    
+      //DrawResultRatio(mcHist, result, Form("finalPlotCheck_%d.eps", (Int_t) eventType));
+  
+      // normalize result
+      result->Scale(TMath::Power(5, etaR) / result->Integral(1, displayRanges[etaR]));
+    
+      result->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
+      //result->SetBinContent(1, 0); result->SetBinError(1, 0);
+      //result->SetTitle(Form(""));
+      result->SetMarkerStyle(0);
+      result->SetLineColor(1);
+      result->SetLineWidth(2);
+      //result->SetMarkerStyle(4);
+      //result->SetStats(kFALSE);
+    
+      // systematic error
+      TH1* error = SystematicsSummary(tpc, (eventType == AliMultiplicityCorrection::kNSD));
+      
+      TH1* systError = (TH1*) result->Clone("systError");
+      // small hack until we have syst. errors for all eta ranges
+      Float_t factor = 80.0 / displayRanges[etaR];
+      for (Int_t i=1; i<=systError->GetNbinsX(); ++i)
+      {
+        systError->SetBinError(i, systError->GetBinContent(i) * error->GetBinContent(factor * i));
+      }
+    
+      // change error drawing style
+      systError->SetFillColor(15);
+      
+      if (eventType == AliMultiplicityCorrection::kNSD)
+      {
+        result->SetLineColor(2);
+        result->SetMarkerColor(2);
+        result->SetMarkerStyle(5);
+      }
+      
+      canvas->cd();
+      systError->GetXaxis()->SetRangeUser(0, displayRanges[etaR]);
+      systError->DrawCopy("E2 ][ SAME");
+      list.Add(result);
+      
+      if (eventType == AliMultiplicityCorrection::kINEL)
+      {
+        TLatex* Tl = new TLatex;
+        Tl->SetTextSize(0.04);
+        //Tl->SetBit(TLatex::kTextNDC);
+        Tl->SetTextAlign(32);
+
+        Float_t etaRangeArr[] = { 0.5, 1.0, 1.4 };
+        TString tmpStr;
+        tmpStr.Form("|#eta| < %.1f (x %d)", etaRangeArr[etaR], (Int_t) TMath::Power(5, etaR));
+        if (etaR == 0)
+          Tl->DrawLatex(36, result->GetBinContent(41), tmpStr);
+        if (etaR == 1)
+        {
+          Tl->SetTextAlign(12);
+          Tl->DrawLatex(82, result->GetBinContent(81), tmpStr);
+        }
+        if (etaR == 2)
+        {
+          Tl->SetTextAlign(12);
+          Tl->DrawLatex(106, result->GetBinContent(101), tmpStr);
+        }
+      }
+      
+      if (etaR == 0)
+        legend->AddEntry(result, (eventType == AliMultiplicityCorrection::kINEL) ? "Inelastic events" : "NSD events", (eventType == AliMultiplicityCorrection::kINEL) ? "L" : "P");
+      
+      count++;
+    }
+  }
+  
+  for (Int_t i=0; i<list.GetEntries(); i++)
+    ((TH1*) list.At(i))->DrawCopy("SAME E ][");
+  
+  legend->Draw();
+
+  canvas->SaveAs(canvas->GetName());
+}
+
 TMatrixD* NonInvertable()
 {
   const Int_t kSize = 5;