removing PWG0depHelper
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 Dec 2007 18:35:26 +0000 (18:35 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 Dec 2007 18:35:26 +0000 (18:35 +0000)
moved dNdEta analysis to task

19 files changed:
PWG0/dNdEta/AliMultiplicityCorrection.cxx
PWG0/dNdEta/AliMultiplicityCorrection.h
PWG0/dNdEta/AliMultiplicityMCSelector.cxx
PWG0/dNdEta/AlidNdEtaCorrection.cxx
PWG0/dNdEta/AlidNdEtaCorrection.h
PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx
PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx [new file with mode: 0644]
PWG0/dNdEta/AlidNdEtaCorrectionTask.h [new file with mode: 0644]
PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx
PWG0/dNdEta/AlidNdEtaTask.cxx [new file with mode: 0644]
PWG0/dNdEta/AlidNdEtaTask.h [new file with mode: 0644]
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/dNdEtaAnalysis.h
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/plotsMultiplicity.C
PWG0/dNdEta/run.C [new file with mode: 0644]
PWG0/dNdEta/runCorrection.C [new file with mode: 0644]
PWG0/dNdEta/runMultiplicitySelector.C
PWG0/dNdEta/rundNdEtaAnalysis.C

index 59ebee5..0041897 100644 (file)
@@ -108,7 +108,7 @@ void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
 
 //____________________________________________________________________
 AliMultiplicityCorrection::AliMultiplicityCorrection() :
-  TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
+  TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
 {
   //
   // default constructor
@@ -198,6 +198,8 @@ AliMultiplicityCorrection::~AliMultiplicityCorrection()
   // Destructor
   //
 
+  Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called");
+
   for (Int_t i = 0; i < kESDHists; ++i)
   {
     if (fMultiplicityESD[i])
@@ -594,7 +596,7 @@ Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
       chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
   }
 
-  return 10.0 + chi2;
+  return 100.0 + chi2;
 }
 
 //____________________________________________________________________
@@ -668,11 +670,14 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c
   // normal way is like this:
   // measGuessVector *= (*fgCorrelationCovarianceMatrix);
   // optimized way like this:
-  for (Int_t i=0; i<fgNParamsMinuit; ++i)
+  for (Int_t i=0; i<fgkMaxInput; ++i)
     measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 
   //measGuessVector.Print();
 
+  // reduce influence of first bin
+  copy(0) *= 0.01;
+
   // (Ad - m) W (Ad - m)
   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
   // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
@@ -686,6 +691,73 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c
 }
 
 //____________________________________________________________________
+void AliMultiplicityCorrection::DrawGuess(Double_t *params)
+{
+  //
+  // draws residuals of solution suggested by params and effect of regularization
+  //
+
+  // d
+  static TVectorD paramsVector(fgNParamsMinuit);
+  for (Int_t i=0; i<fgNParamsMinuit; ++i)
+    paramsVector[i] = params[i] * params[i];
+
+  // Ad
+  TVectorD measGuessVector(fgkMaxInput);
+  measGuessVector = (*fgCorrelationMatrix) * paramsVector;
+
+  // Ad - m
+  measGuessVector -= (*fgCurrentESDVector);
+
+  TVectorD copy(measGuessVector);
+
+  // (Ad - m) W
+  // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
+  // normal way is like this:
+  // measGuessVector *= (*fgCorrelationCovarianceMatrix);
+  // optimized way like this:
+  for (Int_t i=0; i<fgkMaxInput; ++i)
+    measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
+
+  // (Ad - m) W (Ad - m)
+  // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
+  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
+  //Double_t chi2FromFit = measGuessVector * copy * 1e6;
+
+  // draw residuals
+  TH1* residuals = new TH1F("residuals", "residuals", fgkMaxInput, -0.5, fgkMaxInput - 0.5);
+  for (Int_t i=0; i<fgkMaxInput; ++i)
+    residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
+  new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
+
+  // draw penalty
+  if (fgRegularizationType != kPol1) {
+    Printf("Drawing not supported");
+    return;
+  }
+
+  TH1* penalty = new TH1F("penalty", "penalty", fgNParamsMinuit, -0.5, fgNParamsMinuit - 0.5);
+
+  for (Int_t i=2; i<fgNParamsMinuit; ++i)
+  {
+    if (params[i-1] == 0)
+      continue;
+
+    Double_t right  = params[i];
+    Double_t middle = params[i-1];
+    Double_t left   = params[i-2];
+
+    Double_t der1 = (right - middle);
+    Double_t der2 = (middle - left);
+
+    Double_t diff = (der1 - der2) / middle;
+
+    penalty->SetBinContent(i, diff * diff);
+  }
+  new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
+}
+
+//____________________________________________________________________
 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
 {
   //
@@ -723,17 +795,17 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
 
   if (createBigBin)
   {
-    Int_t maxBin = 0;
+    fLastBinLimit = 0;
     for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
     {
       if (fCurrentESD->GetBinContent(i) <= 5)
       {
-        maxBin = i;
+        fLastBinLimit = i;
         break;
       }
     }
 
-    if (maxBin > 0)
+    if (fLastBinLimit > 0)
     {
       TCanvas* canvas = 0;
 
@@ -757,25 +829,25 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
         fCurrentCorrelation->DrawCopy("COLZ");
       }
 
-      printf("Bin limit in measured spectrum is %d.\n", maxBin);
-      fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
-      for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i)
+      printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit);
+      fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX()));
+      for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i)
       {
         fCurrentESD->SetBinContent(i, 0);
         fCurrentESD->SetBinError(i, 0);
       }
       // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
-      fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin)));
+      fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit)));
 
-      printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin));
+      printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit));
 
       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
       {
-        fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY()));
+        fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY()));
         // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
-        fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin)));
+        fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit)));
 
-        for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+        for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
         {
           fCurrentCorrelation->SetBinContent(i, j, 0);
           fCurrentCorrelation->SetBinError(i, j, 0);
@@ -945,18 +1017,17 @@ Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEffici
   for (Int_t i=0; i<fgNParamsMinuit; i++)
     (*fgEntropyAPriori)[i] = 1;
 
-  if (initialConditions)
-  {
+  if (!initialConditions) {
+    initialConditions = measured;
+  } else {
     printf("Using different starting conditions...\n");
     //new TCanvas; initialConditions->DrawCopy();
-
     initialConditions->Scale(1.0 / initialConditions->Integral());
-    for (Int_t i=0; i<fgNParamsMinuit; i++)
-      if (initialConditions->GetBinContent(i+1) > 0)
-        (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
   }
-  else
-    initialConditions = measured;
+
+  for (Int_t i=0; i<fgNParamsMinuit; i++)
+    if (initialConditions->GetBinContent(i+1) > 0)
+      (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
 
   Double_t* results = new Double_t[fgNParamsMinuit+1];
   for (Int_t i=0; i<fgNParamsMinuit; ++i)
@@ -991,15 +1062,17 @@ Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEffici
       (*fgCorrelationCovarianceMatrix)(i, i) = 0;
   }
 
-  Int_t dummy;
-  Double_t chi2;
+  Int_t dummy = 0;
+  Double_t chi2 = 0;
   MinuitFitFunction(dummy, 0, chi2, results, 0);
+  DrawGuess(results);
   printf("Chi2 of initial parameters is = %f\n", chi2);
 
-  delete[] results;
-
   if (check)
+  {
+    delete[] results;
     return -1;
+  }
 
   // first param is number of iterations, second is precision....
   arglist[0] = 1e6;
@@ -1016,9 +1089,10 @@ Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEffici
   {
     if (aEfficiency->GetBinContent(i+1) > 0)
     {
-      result->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / aEfficiency->GetBinContent(i+1));
+      results[i] = minuit->GetParameter(i);
+      result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1));
       // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
-      result->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) /  aEfficiency->GetBinContent(i+1));
+      result->SetBinError(i+1, minuit->GetParError(i) * results[i] /  aEfficiency->GetBinContent(i+1));
     }
     else
     {
@@ -1026,6 +1100,9 @@ Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEffici
       result->SetBinError(i+1, 0);
     }
   }
+  DrawGuess(results);
+
+  delete[] results;
 
   return status;
 }
@@ -1233,7 +1310,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   delete ratio;
   mcHist->Scale(scalingFactor);*/
 
-  // find bin with <= 100 entries. this is used as limit for the chi2 calculation
+  // 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)
   {
@@ -1280,8 +1357,9 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     if (fCurrentESD->GetBinError(i) > 0)
     {
       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
-      if (i > 1)
+      if (i > 1 && i <= fLastBinLimit)
         chi2 += value * value;
+      Printf("%d --> %f (%f)", i, value * value, chi2);
       residual->SetBinContent(i, value);
       residualHist->Fill(value);
     }
@@ -1292,17 +1370,15 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     }
     convolutedProj->SetBinError(i, 0);
     residual->SetBinError(i, 0);
-
-    if (i == 200)
-      fLastChi2Residuals = chi2;
   }
+  fLastChi2Residuals = chi2;
 
   new TCanvas; residualHist->DrawCopy();
 
   //residualHist->Fit("gaus", "N");
   //delete residualHist;
 
-  printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
+  printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit);
 
   TCanvas* canvas1 = 0;
   if (simple)
@@ -1322,12 +1398,17 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   proj->GetXaxis()->SetRangeUser(0, 200);
   proj->SetTitle(";true multiplicity;Entries");
   proj->SetStats(kFALSE);
-  proj->DrawCopy("HIST");
-  gPad->SetLogy();
 
-  //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
-  fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
+
+  if (proj->GetEntries() > 0) {
+    proj->DrawCopy("HIST");
+    fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
+  }
+  else
+    fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
+
+  gPad->SetLogy();
 
   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
   legend->AddEntry(proj, "true distribution");
@@ -1506,10 +1587,10 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     fftResult2->GetYaxis()->UnZoom();
     fftResult3->GetYaxis()->UnZoom();
 
-    Printf("AFTER FFT");
+    //Printf("AFTER FFT");
     for (Int_t i=0; i<kFFT; ++i)
     {
-      Printf("%d: %f", i, fftReal[i]);
+      //Printf("%d: %f", i, fftReal[i]);
       fftResult->SetBinContent(i+1, fftReal[i]);
       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
       {
@@ -1522,10 +1603,10 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     fftResult->DrawCopy();
 
 
-    Printf("IMAG");
+    //Printf("IMAG");
     for (Int_t i=0; i<kFFT; ++i)
     {
-      Printf("%d: %f", i, fftImag[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]));
@@ -1549,7 +1630,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
         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]);
+        //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
       }
     }
 
@@ -1558,10 +1639,10 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
     fftResult4->Reset();
 
-    Printf("AFTER BACK-TRAFO");
+    //Printf("AFTER BACK-TRAFO");
     for (Int_t i=0; i<kFFT; ++i)
     {
-      Printf("%d: %f", i, fftReal[i]);
+      //Printf("%d: %f", i, fftReal[i]);
       fftResult4->SetBinContent(i+1+10, fftReal[i]);
     }
 
@@ -1765,6 +1846,41 @@ TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
 }
 
 //____________________________________________________________________
+TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
+{
+  // calculate standard deviation of (results[0] - results[k]) k=1...max-1
+  // per bin one gets: sigma(r[0] - r[n]) / r[0]
+
+  TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
+  standardDeviation->Reset();
+
+  for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
+  {
+    if (results[0]->GetBinContent(x) > 0)
+    {
+      Double_t average = 0;
+      for (Int_t n=1; n<max; ++n)
+        average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
+      average /= max-1;
+
+      Double_t variance = 0;
+      for (Int_t n=1; n<max; ++n)
+      {
+        Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
+        variance += value * value;
+      }
+      variance /= max-1;
+
+      Double_t standardDev = TMath::Sqrt(variance);
+      standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
+      //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
+    }
+  }
+
+  return standardDeviation;
+}
+
+//____________________________________________________________________
 TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
 {
   //
@@ -1950,32 +2066,7 @@ TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, In
 //     }
 //   }
 
-  // calculate standard deviation of (results[0] - results[n])
-  TH1* standardDeviation = (TH1*) maxError->Clone("standardDeviation");
-  standardDeviation->Reset();
-
-  for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
-  {
-    if (results[0]->GetBinContent(x) > 0)
-    {
-      Double_t average = 0;
-      for (Int_t n=1; n<kErrorIterations; ++n)
-        average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
-      average /= kErrorIterations-1;
-
-      Double_t variance = 0;
-      for (Int_t n=1; n<kErrorIterations; ++n)
-      {
-        Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
-        variance += value * value;
-      }
-      variance /= kErrorIterations-1;
-
-      Double_t standardDev = TMath::Sqrt(variance);
-      standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
-      Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
-    }
-  }
+  TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
 
   // compare maxError to RMS of profile (variable name: average)
   // first: calculate rms in percent of value
@@ -2073,21 +2164,15 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
-  // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
+  // spectrum calculated. This is performed N times and the sigma is taken as the statistical
   // error of the unfolding method itself.
 
-  TH1* maxError = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("maxError");
-  maxError->Reset();
-
-  TH1* normalizedResult = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("normalizedResult");
-  normalizedResult->Scale(1.0 / normalizedResult->Integral());
-
   const Int_t kErrorIterations = 20;
 
   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
 
   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
-  TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
+  TH1* resultArray[kErrorIterations+1];
   for (Int_t n=0; n<kErrorIterations; ++n)
   {
     // randomize the content of clone following a poisson with the mean = the value of that bin
@@ -2099,29 +2184,27 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
       randomized->SetBinError(x, TMath::Sqrt(randomValue));
     }
 
+    TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
     result2->Reset();
     if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
       return;
 
     result2->Scale(1.0 / result2->Integral());
 
-    // calculate ratio
-    result2->Divide(normalizedResult);
-
-    //new TCanvas; ratio->DrawCopy("HIST");
-
-    // find max. deviation
-    for (Int_t x=1; x<=result2->GetNbinsX(); x++)
-      maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - result2->GetBinContent(x))));
+    resultArray[n+1] = result2;
   }
   delete randomized;
-  delete result2;
+
+  resultArray[0] = fMultiplicityESDCorrected[correlationID];
+  TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
+
+  for (Int_t n=0; n<kErrorIterations; ++n)
+    delete resultArray[n+1];
 
   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, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
 
-  delete maxError;
-  delete normalizedResult;
+  delete error;
 }
 
 //____________________________________________________________________
@@ -2139,11 +2222,15 @@ Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEffi
 
   const Int_t kStartBin = 0;
 
-  const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
+  const Int_t kMaxM = fgkMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
   const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
 
+  // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
+  const Double_t kConvergenceLimit = kMaxT * 1e-6;
+
   // store information in arrays, to increase processing speed (~ factor 5)
   Double_t measuredCopy[kMaxM];
+  Double_t measuredError[kMaxM];
   Double_t prior[kMaxT];
   Double_t result[kMaxT];
   Double_t efficiency[kMaxT];
@@ -2156,6 +2243,7 @@ Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEffi
   for (Int_t m=0; m<kMaxM; m++)
   {
     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
+    measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
 
     for (Int_t t=0; t<kMaxT; t++)
     {
@@ -2181,22 +2269,30 @@ Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEffi
       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
   }
 
-  //TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
+  //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
 
   // unfold...
-  for (Int_t i=0; i<nIterations; i++)
+  for (Int_t i=0; i<nIterations || nIterations < 0; i++)
   {
     //printf(" iteration %i \n", i);
 
     // calculate IR from Beyes theorem
     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
 
+    Double_t chi2Measured = 0;
     for (Int_t m=0; m<kMaxM; m++)
     {
       Float_t norm = 0;
       for (Int_t t = kStartBin; t<kMaxT; t++)
         norm += response[t][m] * prior[t];
 
+      // calc. chi2: (measured - response * prior) / error
+      if (measuredError[m] > 0)
+      {
+        Double_t value = (measuredCopy[m] - norm) / measuredError[m];
+        chi2Measured += value * value;
+      }
+
       if (norm > 0)
       {
         for (Int_t t = kStartBin; t<kMaxT; t++)
@@ -2208,6 +2304,7 @@ Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEffi
           inverseResponse[t][m] = 0;
       }
     }
+    //Printf("chi2Measured of the last prior is %e", chi2Measured);
 
     for (Int_t t = kStartBin; t<kMaxT; t++)
     {
@@ -2221,6 +2318,7 @@ Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEffi
         result[t] = 0;
     }
 
+    Double_t chi2LastIter = 0;
     // regularization (simple smoothing)
     for (Int_t t=kStartBin; t<kMaxT; t++)
     {
@@ -2236,32 +2334,29 @@ Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEffi
       else
         newValue = result[t];
 
+      // calculate chi2 (change from last iteration)
+      if (prior[t] > 1e-5)
+      {
+        Double_t diff = (prior[t] - newValue) / prior[t];
+        chi2LastIter += diff * diff;
+      }
+
       prior[t] = newValue;
     }
+    //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
+    //convergence->Fill(i+1, chi2LastIter);
 
-    // calculate chi2 (change from last iteration)
-    //Double_t chi2 = 0;
-    /*for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
+    if (nIterations < 0 && chi2LastIter < kConvergenceLimit)
     {
-      Float_t newValue = hTrueSmoothed->GetBinContent(t);
-      //Float_t diff = hPrior->GetBinContent(t) - newValue;
-      //chi2 += (Double_t) diff * diff;
-
-      hPrior->SetBinContent(t, newValue);
+      Printf("Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
+      break;
     }
-    //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
-    //convergence->Fill(i+1, chi2);
 
     //if (i % 5 == 0)
     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
-
-    delete hTrueSmoothed;*/
   } // end of iterations
 
-  //new TCanvas;
-  //convergence->DrawCopy();
-  //gPad->SetLogy();
-
+  //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
   //delete convergence;
 
   for (Int_t t=0; t<kMaxT; t++)
index fee2576..95df174 100644 (file)
@@ -60,6 +60,7 @@ class AliMultiplicityCorrection : public TNamed {
 
     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);
 
+    static TH1* CalculateStdDev(TH1** results, Int_t max);
     TH1* StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0);
 
     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
@@ -109,6 +110,8 @@ class AliMultiplicityCorrection : public TNamed {
     static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
     static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
 
+    static void DrawGuess(Double_t *params);
+
     void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
 
     Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
@@ -147,6 +150,7 @@ class AliMultiplicityCorrection : public TNamed {
 
     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)
index 04226a0..a07ac82 100644 (file)
@@ -25,7 +25,6 @@
 
 #include "esdTrackCuts/AliESDtrackCuts.h"
 #include "AliPWG0Helper.h"
-#include "AliPWG0depHelper.h"
 #include "dNdEta/AliMultiplicityCorrection.h"
 #include "AliCorrection.h"
 #include "AliCorrectionMatrix3D.h"
@@ -246,7 +245,7 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
   if (fSelectProcessType > 0)
   {
     // getting process information; NB: this only works for Pythia
-    Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
+    Int_t processtype = AliPWG0Helper::GetPythiaEventProcessType(header);
 
     Bool_t processEvent = kFALSE;
 
@@ -483,7 +482,7 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
       // preserve label for later
       Int_t labelCopy = label;
       if (labelCopy >= 0)
-        labelCopy = AliPWG0depHelper::FindPrimaryMotherLabel(stack, labelCopy);
+        labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
       if (labelCopy >= 0)
         mother = stack->Particle(labelCopy);
 
@@ -533,7 +532,7 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
 
       // find mother
       if (label >= 0)
-        motherLabel = AliPWG0depHelper::FindPrimaryMotherLabel(stack, label);
+        motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
       if (motherLabel >= 0)
         mother = stack->Particle(motherLabel);
 
index 3079edd..48def25 100644 (file)
@@ -42,7 +42,7 @@ AlidNdEtaCorrection::AlidNdEtaCorrection()
 }
 
 //____________________________________________________________________
-AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title)
+AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title, const char* analysis)
   : TNamed(name, title),
   fTrack2ParticleCorrection(0),
   fVertexRecoCorrection(0),
@@ -54,12 +54,12 @@ AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title
   // constructor
   //
 
-  fTrack2ParticleCorrection = new AliCorrection("Track2Particle", "Track2Particle");
-  fVertexRecoCorrection     = new AliCorrection("VertexReconstruction", "VertexReconstruction");
+  fTrack2ParticleCorrection = new AliCorrection("Track2Particle", "Track2Particle", analysis);
+  fVertexRecoCorrection     = new AliCorrection("VertexReconstruction", "VertexReconstruction", analysis);
 
-  fTriggerBiasCorrectionMBToINEL = new AliCorrection("TriggerBias_MBToINEL", "TriggerBias_MBToINEL");
-  fTriggerBiasCorrectionMBToNSD  = new AliCorrection("TriggerBias_MBToNSD", "TriggerBias_MBToNSD");
-  fTriggerBiasCorrectionMBToND   = new AliCorrection("TriggerBias_MBToND", "TriggerBias_MBToND");
+  fTriggerBiasCorrectionMBToINEL = new AliCorrection("TriggerBias_MBToINEL", "TriggerBias_MBToINEL", analysis);
+  fTriggerBiasCorrectionMBToNSD  = new AliCorrection("TriggerBias_MBToNSD", "TriggerBias_MBToNSD", analysis);
+  fTriggerBiasCorrectionMBToND   = new AliCorrection("TriggerBias_MBToND", "TriggerBias_MBToND", analysis);
 }
 
 //____________________________________________________________________
index ea53cdf..7d54dc4 100644 (file)
@@ -36,7 +36,7 @@ public:
   };
 
   AlidNdEtaCorrection();
-  AlidNdEtaCorrection(const Char_t* name, const Char_t* title);
+  AlidNdEtaCorrection(const Char_t* name, const Char_t* title, const char* analysis = "TPC");
 
   virtual Long64_t Merge(TCollection* list);
 
index a62e5f9..41b3b30 100644 (file)
@@ -29,7 +29,6 @@
 #include "esdTrackCuts/AliESDtrackCuts.h"
 #include "dNdEta/AlidNdEtaCorrection.h"
 #include "AliPWG0Helper.h"
-#include "AliPWG0depHelper.h"
 
 #include "dNdEta/dNdEtaAnalysis.h"
 
@@ -242,7 +241,7 @@ Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
   genHeader->PrimaryVertex(vtxMC);
 
   // get process type
-  Int_t processType = AliPWG0depHelper::GetPythiaEventProcessType(header);
+  Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
   AliDebug(AliLog::kDebug+1, Form("Found pythia procces type %d", processType));
 
   if (processType<0)
diff --git a/PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx b/PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
new file mode 100644 (file)
index 0000000..67dd58e
--- /dev/null
@@ -0,0 +1,436 @@
+/* $Id$ */
+
+#include "AlidNdEtaCorrectionTask.h"
+
+#include <TCanvas.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TParticle.h>
+
+#include <AliLog.h>
+#include <AliESDVertex.h>
+#include <AliESDEvent.h>
+#include <AliStack.h>
+#include <AliHeader.h>
+#include <AliGenEventHeader.h>
+#include <AliMultiplicity.h>
+#include <AliAnalysisManager.h>
+#include <AliMCEventHandler.h>
+#include <AliMCEvent.h>
+#include <AliESDInputHandler.h>
+
+#include "esdTrackCuts/AliESDtrackCuts.h"
+#include "AliPWG0Helper.h"
+//#include "AliCorrection.h"
+//#include "AliCorrectionMatrix3D.h"
+#include "dNdEta/dNdEtaAnalysis.h"
+#include "dNdEta/AlidNdEtaCorrection.h"
+
+ClassImp(AlidNdEtaCorrectionTask)
+
+AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
+  AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
+  fESD(0),
+  fOutput(0),
+  fOption(opt),
+  fAnalysisMode(kTPC),
+  fSignMode(0),
+  fEsdTrackCuts(0),
+  fdNdEtaCorrection(0),
+  fdNdEtaAnalysisMC(0),
+  fdNdEtaAnalysisESD(0),
+  fPIDParticles(0),
+  fPIDTracks(0)
+{
+  //
+  // Constructor. Initialization of pointers
+  //
+
+  // Define input and output slots here
+  DefineInput(0, TChain::Class());
+  DefineOutput(0, TList::Class());
+}
+
+AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
+{
+  //
+  // Destructor
+  //
+
+  // histograms are in the output list and deleted when the output
+  // list is deleted by the TSelector dtor
+
+  if (fOutput) {
+    delete fOutput;
+    fOutput = 0;
+  }
+}
+
+//________________________________________________________________________
+void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
+{
+  // Connect ESD
+  // Called once
+
+  Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
+
+  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+  if (!tree) {
+    Printf("ERROR: Could not read tree from input slot 0");
+  } else {
+    // Disable all branches and enable only the needed ones
+    tree->SetBranchStatus("*", 0);
+
+    tree->SetBranchStatus("fTriggerMask", 1);
+    tree->SetBranchStatus("fSPDVertex*", 1);
+
+    if (fAnalysisMode == kSPD)
+      tree->SetBranchStatus("fSPDMult*", 1);
+
+    if (fAnalysisMode == kTPC) {
+      AliESDtrackCuts::EnableNeededBranches(tree);
+      tree->SetBranchStatus("fTracks.fLabel", 1);
+    }
+
+    tree->SetCacheSize(0);
+
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+    if (!esdH) {
+      Printf("ERROR: Could not get ESDInputHandler");
+    } else
+      fESD = esdH->GetEvent();
+  }
+}
+
+void AlidNdEtaCorrectionTask::CreateOutputObjects()
+{
+  // create result objects and add to output list
+
+  if (fOption.Contains("only-positive"))
+  {
+    Printf("INFO: Processing only positive particles.");
+    fSignMode = 1;
+  }
+  else if (fOption.Contains("only-negative"))
+  {
+    Printf("INFO: Processing only negative particles.");
+    fSignMode = -1;
+  }
+
+  TString detector;
+  if (fAnalysisMode == kTPC)
+  {
+    detector = "TPC";
+  }
+  else if (fAnalysisMode == kSPD)
+    detector = "SPD";
+
+  fOutput = new TList;
+  fOutput->SetOwner();
+
+  fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", detector);
+  fOutput->Add(fdNdEtaCorrection);
+
+  fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
+  fOutput->Add(fPIDParticles);
+
+  fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
+  fOutput->Add(fPIDTracks);
+
+  fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", detector);
+  fOutput->Add(fdNdEtaAnalysisMC);
+
+  fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", detector);
+  fOutput->Add(fdNdEtaAnalysisESD);
+}
+
+void AlidNdEtaCorrectionTask::Exec(Option_t*)
+{
+  // process the event
+
+  // Check prerequisites
+  if (!fESD)
+  {
+    AliDebug(AliLog::kError, "ESD branch not available");
+    return;
+  }
+
+  // trigger definition
+  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
+
+  Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD->GetVertex());
+
+  // post the data already here
+  PostData(0, fOutput);
+
+  // create list of (label, eta, pt) tuples
+  Int_t inputCount = 0;
+  Int_t* labelArr = 0;
+  Float_t* etaArr = 0;
+  Float_t* ptArr = 0;
+  if (fAnalysisMode == 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()];
+    ptArr = new Float_t[mult->GetNumberOfTracklets()];
+
+    // 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));
+
+      // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
+      if (mult->GetDeltaPhi(i) < -1000)
+        continue;
+
+      etaArr[inputCount] = mult->GetEta(i);
+      labelArr[inputCount] = mult->GetLabel(i);
+      ptArr[inputCount] = 0; // no pt for tracklets
+      ++inputCount;
+    }
+  }
+  else if (fAnalysisMode == kTPC)
+  {
+    if (!fEsdTrackCuts)
+    {
+      AliDebug(AliLog::kError, "fESDTrackCuts not available");
+      return;
+    }
+
+    // get multiplicity from ESD tracks
+    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+    Int_t nGoodTracks = list->GetEntries();
+
+    labelArr = new Int_t[nGoodTracks];
+    etaArr = new Float_t[nGoodTracks];
+    ptArr = new Float_t[nGoodTracks];
+
+    // loop over esd tracks
+    for (Int_t i=0; i<nGoodTracks; i++)
+    {
+      AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
+      if (!esdTrack)
+      {
+        AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
+        continue;
+      }
+
+      etaArr[inputCount] = esdTrack->Eta();
+      labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
+      ptArr[inputCount] = esdTrack->Pt();
+      ++inputCount;
+    }
+  }
+  else
+    return;
+
+  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+  if (!eventHandler) {
+    Printf("ERROR: Could not retrieve MC event handler");
+    return;
+  }
+
+  AliMCEvent* mcEvent = eventHandler->MCEvent();
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }
+
+  AliStack* stack = mcEvent->Stack();
+  if (!stack)
+  {
+    AliDebug(AliLog::kError, "Stack not available");
+    return;
+  }
+
+  AliHeader* header = mcEvent->Header();
+  if (!header)
+  {
+    AliDebug(AliLog::kError, "Header not available");
+    return;
+  }
+
+  // get process type
+  Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
+  AliDebug(AliLog::kDebug+1, Form("Found pythia procces type %d", processType));
+
+  if (processType<0)
+    AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
+
+  // get the MC vertex
+  AliGenEventHeader* genHeader = header->GenEventHeader();
+  TArrayF vtxMC(3);
+  genHeader->PrimaryVertex(vtxMC);
+
+  // loop over mc particles
+  Int_t nPrim  = stack->GetNprimary();
+
+  for (Int_t iMc = 0; iMc < nPrim; ++iMc)
+  {
+    //Printf("Getting particle %d", iMc);
+    TParticle* particle = stack->Particle(iMc);
+
+    if (!particle)
+      continue;
+
+    if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+      continue;
+
+    if (SignOK(particle->GetPDG()) == kFALSE)
+      continue;
+
+    Float_t eta = particle->Eta();
+    Float_t pt = particle->Pt();
+
+    fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
+
+    if (eventTriggered)
+      if (eventVertex)
+        fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
+  }
+
+  for (Int_t i=0; i<inputCount; ++i)
+  {
+    Int_t label = labelArr[i];
+
+    if (label < 0)
+    {
+      AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", label));
+      continue;
+    }
+
+    TParticle* particle = stack->Particle(label);
+    if (!particle)
+    {
+      AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
+      continue;
+    }
+
+    if (SignOK(particle->GetPDG()) == kFALSE)
+        continue;
+
+    if (eventTriggered && eventVertex)
+    {
+      fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
+      fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+      if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
+      {
+        fPIDTracks->Fill(particle->GetPdgCode());
+      }
+    }
+  }
+
+  if (eventTriggered && eventVertex)
+    fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
+
+  // stuff regarding the vertex reco correction and trigger bias correction
+  fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
+  if (eventTriggered) {
+    if (eventVertex)
+    {
+      fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
+    }
+  }
+
+  delete[] etaArr;
+  delete[] labelArr;
+  delete[] ptArr;
+}
+
+void AlidNdEtaCorrectionTask::Terminate(Option_t *)
+{
+  // The Terminate() function is the last function to be called during
+  // a query. It always runs on the client, it can be used to present
+  // the results graphically or save the results to file.
+
+  fOutput = dynamic_cast<TList*> (GetOutputData(0));
+  if (!fOutput) {
+    Printf("ERROR: fOutput not available");
+    return;
+  }
+
+  fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
+  fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
+  fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
+  if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
+  {
+    AliDebug(AliLog::kError, "Could not read object from output list");
+    return;
+  }
+
+  fdNdEtaCorrection->Finish();
+
+  TString fileName;
+  fileName.Form("correction_map%s.root", fOption.Data());
+  TFile* fout = new TFile(fileName, "RECREATE");
+
+  if (fEsdTrackCuts)
+    fEsdTrackCuts->SaveHistograms("esd_track_cuts");
+  fdNdEtaCorrection->SaveHistograms();
+  fdNdEtaAnalysisMC->SaveHistograms();
+  fdNdEtaAnalysisESD->SaveHistograms();
+
+  fout->Write();
+  fout->Close();
+
+  fdNdEtaCorrection->DrawHistograms();
+
+  if (fPIDParticles && fPIDTracks)
+  {
+    new TCanvas("pidcanvas", "pidcanvas", 500, 500);
+
+    fPIDParticles->Draw();
+    fPIDTracks->SetLineColor(2);
+    fPIDTracks->Draw("SAME");
+
+    TDatabasePDG* pdgDB = new TDatabasePDG;
+
+    for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
+      if (fPIDParticles->GetBinContent(i) > 0)
+        printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
+
+    delete pdgDB;
+    pdgDB = 0;
+  }
+
+  Printf("Writting result to %s", fileName.Data());
+}
+
+Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
+{
+  // returns if a particle with this sign should be counted
+  // this is determined by the value of fSignMode, which should have the same sign
+  // as the charge
+  // if fSignMode is 0 all particles are counted
+
+  if (fSignMode == 0)
+    return kTRUE;
+
+  if (!particle)
+  {
+    printf("WARNING: not counting a particle that does not have a pdg particle\n");
+    return kFALSE;
+  }
+
+  Double_t charge = particle->Charge();
+
+  if (fSignMode > 0)
+    if (charge < 0)
+      return kFALSE;
+
+  if (fSignMode < 0)
+    if (charge > 0)
+      return kFALSE;
+
+  return kTRUE;
+}
diff --git a/PWG0/dNdEta/AlidNdEtaCorrectionTask.h b/PWG0/dNdEta/AlidNdEtaCorrectionTask.h
new file mode 100644 (file)
index 0000000..36d7519
--- /dev/null
@@ -0,0 +1,59 @@
+/* $Id$ */
+
+#ifndef AlidNdEtaCorrectionTask_H
+#define AlidNdEtaCorrectionTask_H
+
+#include "AliAnalysisTask.h"
+#include <TString.h>
+
+class AliESDtrackCuts;
+class dNdEtaAnalysis;
+class AlidNdEtaCorrection;
+class TH1F;
+class AliESDEvent;
+class TParticlePDG;
+
+class AlidNdEtaCorrectionTask : public AliAnalysisTask {
+  public:
+    enum AnalysisMethod { kSPD = 0, kTPC };
+
+    AlidNdEtaCorrectionTask(const char* opt = "");
+    virtual ~AlidNdEtaCorrectionTask();
+
+    virtual void   ConnectInputData(Option_t *);
+    virtual void   CreateOutputObjects();
+    virtual void   Exec(Option_t*);
+    virtual void   Terminate(Option_t *);
+
+    void SetTrackCuts(AliESDtrackCuts* cuts) { fEsdTrackCuts = cuts; }
+    void SetAnalysisMode(AnalysisMethod mode) { fAnalysisMode = mode; }
+
+ protected:
+    Bool_t SignOK(TParticlePDG* particle);
+
+    AliESDEvent *fESD;               //! ESD object
+    TList* fOutput;                  //! list send on output slot 0
+
+    TString fOption;                 // option string
+    AnalysisMethod fAnalysisMode;    // detector that is used for analysis
+    Int_t fSignMode;                 // if 0 process all particles, if +-1 process only particles with that sign
+
+    AliESDtrackCuts*  fEsdTrackCuts;             // Object containing the parameters of the esd track cuts
+
+    AlidNdEtaCorrection* fdNdEtaCorrection;      //! contains the intermediate histograms (on each slave)
+
+    dNdEtaAnalysis* fdNdEtaAnalysisMC;           //! analysis from MC (only triggered, vertex events)
+    dNdEtaAnalysis* fdNdEtaAnalysisESD;          //! analysis from ESD (not yet corrected!)
+
+    // control histograms
+    TH1F* fPIDParticles;                         //! pid of primary particles
+    TH1F* fPIDTracks;                            //! pid of reconstructed tracks
+
+ private:
+    AlidNdEtaCorrectionTask(const AlidNdEtaCorrectionTask&);
+    AlidNdEtaCorrectionTask& operator=(const AlidNdEtaCorrectionTask&);
+
+  ClassDef(AlidNdEtaCorrectionTask, 1);
+};
+
+#endif
index 20ca14a..2761919 100644 (file)
@@ -23,7 +23,6 @@
 
 #include "esdTrackCuts/AliESDtrackCuts.h"
 #include "AliPWG0Helper.h"
-#include "AliPWG0depHelper.h"
 #include "dNdEta/AlidNdEtaCorrection.h"
 
 ClassImp(AlidNdEtaSystematicsSelector)
@@ -223,7 +222,7 @@ Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
   }
   
   // getting process information NB: this only works for Pythia !!!
-  Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
+  Int_t processtype = AliPWG0Helper::GetPythiaEventProcessType(header);
 
   // can only read pythia headers, either directly or from cocktalil header
   AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader());
@@ -392,7 +391,7 @@ Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
     }
 
     // find primary particle that created this particle
-    TParticle* mother = AliPWG0depHelper::FindPrimaryMother(stack, label);
+    TParticle* mother = AliPWG0Helper::FindPrimaryMother(stack, label);
     if (!mother)
       continue;
 
diff --git a/PWG0/dNdEta/AlidNdEtaTask.cxx b/PWG0/dNdEta/AlidNdEtaTask.cxx
new file mode 100644 (file)
index 0000000..65c9533
--- /dev/null
@@ -0,0 +1,544 @@
+/* $Id$ */
+
+#include "AlidNdEtaTask.h"
+
+#include <TStyle.h>
+#include <TSystem.h>
+#include <TCanvas.h>
+#include <TVector3.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TParticle.h>
+#include <TRandom.h>
+#include <TNtuple.h>
+#include <TObjString.h>
+#include <TF1.h>
+
+#include <AliLog.h>
+#include <AliESDVertex.h>
+#include <AliESDEvent.h>
+#include <AliStack.h>
+#include <AliHeader.h>
+#include <AliGenEventHeader.h>
+#include <AliMultiplicity.h>
+#include <AliAnalysisManager.h>
+#include <AliMCEventHandler.h>
+#include <AliMCEvent.h>
+#include <AliESDInputHandler.h>
+
+#include "esdTrackCuts/AliESDtrackCuts.h"
+#include "AliPWG0Helper.h"
+#include "AliCorrection.h"
+#include "AliCorrectionMatrix3D.h"
+#include "dNdEta/dNdEtaAnalysis.h"
+
+ClassImp(AlidNdEtaTask)
+
+AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
+  AliAnalysisTask("AlidNdEtaTask", ""),
+  fESD(0),
+  fOutput(0),
+  fOption(opt),
+  fAnalysisMode(kTPC),
+  fReadMC(kFALSE),
+  fEsdTrackCuts(0),
+  fdNdEtaAnalysisESD(0),
+  fMult(0),
+  fMultVtx(0),
+  fEvents(0),
+  fdNdEtaAnalysis(0),
+  fdNdEtaAnalysisTr(0),
+  fdNdEtaAnalysisTrVtx(0),
+  fdNdEtaAnalysisTracks(0),
+  fVertex(0),
+  fPartPt(0)
+{
+  //
+  // Constructor. Initialization of pointers
+  //
+
+  // Define input and output slots here
+  DefineInput(0, TChain::Class());
+  DefineOutput(0, TList::Class());
+}
+
+AlidNdEtaTask::~AlidNdEtaTask()
+{
+  //
+  // Destructor
+  //
+
+  // histograms are in the output list and deleted when the output
+  // list is deleted by the TSelector dtor
+
+  if (fOutput) {
+    delete fOutput;
+    fOutput = 0;
+  }
+}
+
+//________________________________________________________________________
+void AlidNdEtaTask::ConnectInputData(Option_t *)
+{
+  // Connect ESD
+  // Called once
+
+  Printf("AlidNdEtaTask::ConnectInputData called");
+
+  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+  if (!tree) {
+    Printf("ERROR: Could not read tree from input slot 0");
+  } else {
+    // Disable all branches and enable only the needed ones
+    //tree->SetBranchStatus("*", 0);
+
+    tree->SetBranchStatus("fTriggerMask", 1);
+    tree->SetBranchStatus("fSPDVertex*", 1);
+
+    if (fAnalysisMode == kSPD)
+      tree->SetBranchStatus("fSPDMult*", 1);
+
+    if (fAnalysisMode == kTPC) {
+      AliESDtrackCuts::EnableNeededBranches(tree);
+      tree->SetBranchStatus("fTracks.fLabel", 1);
+    }
+
+    tree->SetCacheSize(0);
+
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+    if (!esdH) {
+      Printf("ERROR: Could not get ESDInputHandler");
+    } else
+      fESD = esdH->GetEvent();
+  }
+}
+
+void AlidNdEtaTask::CreateOutputObjects()
+{
+  // create result objects and add to output list
+
+  fOutput = new TList;
+  fOutput->SetOwner();
+
+  TString detector;
+  if (fAnalysisMode == kTPC)
+  {
+    detector = "TPC";
+  }
+  else if (fAnalysisMode == kSPD)
+    detector = "SPD";
+
+  fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", detector);
+  fOutput->Add(fdNdEtaAnalysisESD);
+
+  fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
+  fOutput->Add(fMult);
+
+  fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
+  fOutput->Add(fMultVtx);
+
+  for (Int_t i=0; i<3; ++i)
+  {
+    fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
+    fPartEta[i]->Sumw2();
+    fOutput->Add(fPartEta[i]);
+  }
+
+  fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
+  fOutput->Add(fEvents);
+
+  if (fReadMC)
+  {
+    fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", detector);
+    fOutput->Add(fdNdEtaAnalysis);
+
+    fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", detector);
+    fOutput->Add(fdNdEtaAnalysisTr);
+
+    fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", detector);
+    fOutput->Add(fdNdEtaAnalysisTrVtx);
+
+    fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", detector);
+    fOutput->Add(fdNdEtaAnalysisTracks);
+
+    fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
+    fOutput->Add(fVertex);
+
+    fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
+    fPartPt->Sumw2();
+    fOutput->Add(fPartPt);
+  }
+}
+
+void AlidNdEtaTask::Exec(Option_t*)
+{
+  // process the event
+
+  // Check prerequisites
+  if (!fESD)
+  {
+    AliDebug(AliLog::kError, "ESD branch not available");
+    return;
+  }
+
+  // trigger definition
+  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
+
+  Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD->GetVertex());
+
+  // get the ESD vertex
+  const AliESDVertex* vtxESD = fESD->GetVertex();
+  Double_t vtx[3];
+  vtxESD->GetXYZ(vtx);
+
+  // post the data already here
+  PostData(0, fOutput);
+
+  // create list of (label, eta, pt) tuples
+  Int_t inputCount = 0;
+  Int_t* labelArr = 0;
+  Float_t* etaArr = 0;
+  Float_t* ptArr = 0;
+  if (fAnalysisMode == 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()];
+    ptArr = new Float_t[mult->GetNumberOfTracklets()];
+
+    // 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));
+
+      // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
+      if (mult->GetDeltaPhi(i) < -1000)
+        continue;
+
+      etaArr[inputCount] = mult->GetEta(i);
+      labelArr[inputCount] = mult->GetLabel(i);
+      ptArr[inputCount] = 0; // no pt for tracklets
+      ++inputCount;
+    }
+  }
+  else if (fAnalysisMode == kTPC)
+  {
+    if (!fEsdTrackCuts)
+    {
+      AliDebug(AliLog::kError, "fESDTrackCuts not available");
+      return;
+    }
+
+    // get multiplicity from ESD tracks
+    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
+    Int_t nGoodTracks = list->GetEntries();
+
+    labelArr = new Int_t[nGoodTracks];
+    etaArr = new Float_t[nGoodTracks];
+    ptArr = new Float_t[nGoodTracks];
+
+    // loop over esd tracks
+    for (Int_t i=0; i<nGoodTracks; i++)
+    {
+      AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
+      if (!esdTrack)
+      {
+        AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
+        continue;
+      }
+
+      etaArr[inputCount] = esdTrack->Eta();
+      labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
+      ptArr[inputCount] = esdTrack->Pt();
+      ++inputCount;
+    }
+  }
+  else
+    return;
+
+  // Processing of ESD information (always)
+  if (eventTriggered)
+  {
+    // control hists
+    fMult->Fill(inputCount);
+    if (eventVertex)
+      fMultVtx->Fill(inputCount);
+
+    // count all events
+    if (!eventVertex)
+      vtx[2] = 0;
+
+    for (Int_t i=0; i<inputCount; ++i)
+    {
+      Float_t eta = etaArr[i];
+      Float_t pt  = ptArr[i];
+
+      fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
+
+      if (TMath::Abs(vtx[2]) < 20)
+      {
+        fPartEta[0]->Fill(eta);
+
+        if (vtx[2] < 0)
+          fPartEta[1]->Fill(eta);
+        else
+          fPartEta[2]->Fill(eta);
+      }
+    }
+
+    // for event count per vertex
+    fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
+
+    // control hists
+    fEvents->Fill(vtx[2]);
+  }
+
+  if (fReadMC)   // Processing of MC information (optional)
+  {
+    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+    if (!eventHandler) {
+      Printf("ERROR: Could not retrieve MC event handler");
+      return;
+    }
+
+    AliMCEvent* mcEvent = eventHandler->MCEvent();
+    if (!mcEvent) {
+      Printf("ERROR: Could not retrieve MC event");
+      return;
+    }
+
+    AliStack* stack = mcEvent->Stack();
+    if (!stack)
+    {
+      AliDebug(AliLog::kError, "Stack not available");
+      return;
+    }
+
+    AliHeader* header = mcEvent->Header();
+    if (!header)
+    {
+      AliDebug(AliLog::kError, "Header not available");
+      return;
+    }
+
+    // get the MC vertex
+    AliGenEventHeader* genHeader = header->GenEventHeader();
+    TArrayF vtxMC(3);
+    genHeader->PrimaryVertex(vtxMC);
+
+    // loop over mc particles
+    Int_t nPrim  = stack->GetNprimary();
+
+    Int_t nAcceptedParticles = 0;
+
+    for (Int_t iMc = 0; iMc < nPrim; ++iMc)
+    {
+      //Printf("Getting particle %d", iMc);
+      TParticle* particle = stack->Particle(iMc);
+
+      if (!particle)
+        continue;
+
+      if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+        continue;
+
+      AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
+      ++nAcceptedParticles;
+
+      Float_t eta = particle->Eta();
+      Float_t pt = particle->Pt();
+
+      fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
+      fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
+
+      if (eventTriggered)
+      {
+        fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
+        if (eventVertex)
+          fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
+      }
+
+      if (TMath::Abs(eta) < 0.8)
+        fPartPt->Fill(pt);
+    }
+
+    fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
+    if (eventTriggered)
+    {
+      fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
+      if (eventVertex)
+        fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
+    }
+
+    if (eventTriggered && eventVertex)
+    {
+      // from tracks is only done for triggered and vertex reconstructed events
+
+      for (Int_t i=0; i<inputCount; ++i)
+      {
+        Int_t label = labelArr[i];
+
+        if (label < 0)
+          continue;
+
+        //Printf("Getting particle of track %d", label);
+        TParticle* particle = stack->Particle(label);
+
+        if (!particle)
+        {
+          AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
+          continue;
+        }
+
+        fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
+      } // end of track loop
+
+      // for event count per vertex
+      fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
+    }
+  }
+
+  delete[] etaArr;
+  delete[] labelArr;
+  delete[] ptArr;
+}
+
+void AlidNdEtaTask::Terminate(Option_t *)
+{
+  // The Terminate() function is the last function to be called during
+  // a query. It always runs on the client, it can be used to present
+  // the results graphically or save the results to file.
+
+  fOutput = dynamic_cast<TList*> (GetOutputData(0));
+  if (!fOutput) {
+    Printf("ERROR: fOutput not available");
+    return;
+  }
+
+  fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
+  fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
+  fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
+  fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
+  for (Int_t i=0; i<3; ++i)
+    fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
+
+  if (!fdNdEtaAnalysisESD)
+  {
+    AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
+    return;
+  }
+
+  if (fMult && fMultVtx)
+  {
+    new TCanvas;
+    fMult->Draw();
+    fMultVtx->SetLineColor(2);
+    fMultVtx->DrawCopy("SAME");
+  }
+
+  if (fPartEta[0])
+  {
+    Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
+    Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
+
+    fPartEta[0]->Scale(1.0 / (events1 + events2));
+    fPartEta[1]->Scale(1.0 / events1);
+    fPartEta[2]->Scale(1.0 / events2);
+
+    for (Int_t i=0; i<3; ++i)
+      fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
+
+    new TCanvas("control", "control", 500, 500);
+    for (Int_t i=0; i<3; ++i)
+    {
+      fPartEta[i]->SetLineColor(i+1);
+      fPartEta[i]->Draw((i==0) ? "" : "SAME");
+    }
+  }
+
+  if (fEvents)
+  {
+    new TCanvas("control3", "control3", 500, 500);
+    fEvents->DrawCopy();
+  }
+
+  TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
+
+  if (fdNdEtaAnalysisESD)
+    fdNdEtaAnalysisESD->SaveHistograms();
+
+  if (fEsdTrackCuts)
+    fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
+
+  if (fMult)
+    fMult->Write();
+
+  if (fMultVtx)
+    fMultVtx->Write();
+
+  for (Int_t i=0; i<3; ++i)
+    if (fPartEta[i])
+      fPartEta[i]->Write();
+
+  if (fEvents)
+    fEvents->Write();
+
+  fout->Write();
+  fout->Close();
+
+  Printf("Writting result to analysis_esd_raw.root");
+
+  if (fReadMC)
+  {
+    fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
+    fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
+    fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
+    fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
+    fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
+
+    if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
+    {
+      AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
+      return;
+    }
+
+    fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
+    fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
+    fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
+    fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
+
+    Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
+    fPartPt->Scale(1.0/events);
+    fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
+
+    TFile* fout = new TFile("analysis_mc.root","RECREATE");
+
+    fdNdEtaAnalysis->SaveHistograms();
+    fdNdEtaAnalysisTr->SaveHistograms();
+    fdNdEtaAnalysisTrVtx->SaveHistograms();
+    fdNdEtaAnalysisTracks->SaveHistograms();
+    fPartPt->Write();
+
+    fout->Write();
+    fout->Close();
+
+    Printf("Writting result to analysis_mc.root");
+
+    if (fPartPt)
+    {
+      new TCanvas("control2", "control2", 500, 500);
+      fPartPt->DrawCopy();
+    }
+  }
+}
diff --git a/PWG0/dNdEta/AlidNdEtaTask.h b/PWG0/dNdEta/AlidNdEtaTask.h
new file mode 100644 (file)
index 0000000..6489d9e
--- /dev/null
@@ -0,0 +1,66 @@
+/* $Id$ */
+
+#ifndef AlidNdEtaTask_H
+#define AlidNdEtaTask_H
+
+#include "AliAnalysisTask.h"
+#include <TString.h>
+
+class AliESDtrackCuts;
+class dNdEtaAnalysis;
+class TH1F;
+class TH3F;
+class AliESDEvent;
+
+class AlidNdEtaTask : public AliAnalysisTask {
+  public:
+    enum AnalysisMethod { kSPD = 0, kTPC };
+
+    AlidNdEtaTask(const char* opt = "");
+    virtual ~AlidNdEtaTask();
+
+    virtual void   ConnectInputData(Option_t *);
+    virtual void   CreateOutputObjects();
+    virtual void   Exec(Option_t*);
+    virtual void   Terminate(Option_t*);
+
+    void SetTrackCuts(AliESDtrackCuts* cuts) { fEsdTrackCuts = cuts; }
+    void SetAnalysisMode(AnalysisMethod mode) { fAnalysisMode = mode; }
+    void SetReadMC(Bool_t flag = kTRUE) { fReadMC = flag; }
+
+ protected:
+    AliESDEvent *fESD;    //! ESD object
+    TList* fOutput;                  //! list send on output slot 0
+
+    TString fOption;      // option string
+    AnalysisMethod fAnalysisMode; // detector that is used for analysis
+
+    Bool_t  fReadMC;       // if true reads MC data (to build correlation maps)
+
+    AliESDtrackCuts* fEsdTrackCuts;         // Object containing the parameters of the esd track cuts
+
+    // Gathered from ESD
+    dNdEtaAnalysis* fdNdEtaAnalysisESD;     //! contains the dndeta from the ESD
+    // control hists
+    TH1F* fMult;                            //! raw multiplicity histogram (control histogram)
+    TH1F* fMultVtx;                            //! raw multiplicity histogram of evts with vtx (control histogram)
+    TH1F* fPartEta[3];            //! counted particles as function of eta (full vertex range, below 0 range, above 0 range)
+    TH1F* fEvents;                //! events counted as function of vtx
+
+    // Gathered from MC (when fReadMC is set)
+    dNdEtaAnalysis* fdNdEtaAnalysis;        //! contains the dndeta from the full sample
+    dNdEtaAnalysis* fdNdEtaAnalysisTr;      //! contains the dndeta from the triggered events
+    dNdEtaAnalysis* fdNdEtaAnalysisTrVtx;   //! contains the dndeta from the triggered events with vertex
+    dNdEtaAnalysis* fdNdEtaAnalysisTracks;  //! contains the dndeta from the triggered events with vertex counted from the mc particles associated to the tracks (comparing this to the raw values from the esd shows the effect of the detector resolution)
+    // the following are control histograms to check the dNdEtaAnalysis class
+    TH3F* fVertex;                //! vertex of counted particles
+    TH1F* fPartPt;                //! counted particles as function of pt
+
+ private:
+    AlidNdEtaTask(const AlidNdEtaTask&);
+    AlidNdEtaTask& operator=(const AlidNdEtaTask&);
+
+  ClassDef(AlidNdEtaTask, 1);
+};
+
+#endif
index 7f2966a..b4341dd 100644 (file)
@@ -39,16 +39,17 @@ dNdEtaAnalysis::dNdEtaAnalysis() :
 }
 
 //____________________________________________________________________
-dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
+dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title, const char* analysis) :
   TNamed(name, title),
   fData(0),
   fPtDist(0)
 {
   // constructor
 
+  // TODO this binning has to be the same than in AliCorrection, somehow passed?!
   Float_t binLimitsPt[] = {0.0, 0.05, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
 
-  fData = new AliCorrection("Analysis", Form("%s Analysis", title));
+  fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysis);
 
   // do not add this hists to the directory
   Bool_t oldStatus = TH1::AddDirectoryStatus();
@@ -184,11 +185,11 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
     if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
       trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
 
-    if (correctionType >= AlidNdEtaCorrection::kVertexReco)
+    /*if (correctionType >= AlidNdEtaCorrection::kVertexReco)
     {
       trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
       eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
-    }
+    }*/
 
     switch (correctionType)
     {
index c8568d5..ad8e70f 100644 (file)
@@ -32,7 +32,7 @@ public:
   enum { kVertexBinning = 1+2 }; // the first is for the whole vertex range, the others divide the vertex range
 
   dNdEtaAnalysis();
-  dNdEtaAnalysis(Char_t* name, Char_t* title);
+  dNdEtaAnalysis(Char_t* name, Char_t* title, const char* analysis = "tpc");
   virtual ~dNdEtaAnalysis();
 
   dNdEtaAnalysis(const dNdEtaAnalysis &c);
index 0258de6..3723788 100644 (file)
@@ -248,7 +248,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMBVtxNoPt->SetMarkerStyle(22);
   histESDMBTracksNoPt->SetMarkerStyle(23);
 
-  TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 3.1, histESDMBVtx->GetMaximum() * 1.1);
+  TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 0, histESDMBVtx->GetMaximum() * 1.1);
   Prepare1DPlot(dummy);
   dummy->SetStats(kFALSE);
   dummy->SetXTitle("#eta");
@@ -266,7 +266,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMBVtxNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
   histESDMBTracksNoPt->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
 
-  /*TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
+  TCanvas* canvas = new TCanvas("dNdEta1", "dNdEta1", 500, 500);
 
   dummy->DrawCopy();
   histESDMBVtx->Draw("SAME");
@@ -274,7 +274,7 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESD->Draw("SAME");
 
   canvas->SaveAs("dNdEta1.gif");
-  canvas->SaveAs("dNdEta1.eps");*/
+  canvas->SaveAs("dNdEta1.eps");
 
   if (onlyESD)
     return;
index c452855..439eea7 100644 (file)
@@ -31,7 +31,7 @@
 const char* correctionFile = "multiplicityMC_2M.root";
 const char* measuredFile   = "multiplicityMC_1M_3.root";
 Int_t etaRange = 3;
-Int_t displayRange = 150; // axis range
+Int_t displayRange = 200; // axis range
 Int_t ratioRange = 151;   // range to calculate difference
 Int_t longDisplayRange = 200;
 
@@ -729,6 +729,9 @@ void bayesianExample()
   //mult->DrawComparison("bayesianExample", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
   DrawResultRatio(mcHist, result, "bayesianExample.eps");
 
+  //Printf("KolmogorovTest says PROB = %f", mcHist->KolmogorovTest(result, "D"));
+  //Printf("Chi2Test says PROB = %f", mcHist->Chi2Test(result));
+
   // draw residual plot
 
   // TODO take out efficiency correction if other than AliMultiplicityCorrection::kTrVtx
@@ -1597,7 +1600,8 @@ void StatisticalUncertaintyCompare(const char* det = "SPD")
   errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
   errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
   errorResponse->SetStats(kFALSE);
-  errorResponse->SetTitle(";true multiplicity;Uncertainty");
+  errorResponse->GetYaxis()->SetTitleOffset(1.2);
+  errorResponse->SetTitle(";true multiplicity;#sigma(U-T)/T");
 
   errorResponse->Draw();
 
@@ -1625,7 +1629,7 @@ void StatisticalUncertaintyCompare(const char* det = "SPD")
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 }
 
-void EfficiencyComparison(Int_t eventType = 2)
+void EfficiencyComparison(Int_t eventType = 2, Bool_t uncertainty = kTRUE)
 {
  const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
 
@@ -1699,17 +1703,19 @@ void EfficiencyComparison(Int_t eventType = 2)
       for (Int_t bin=1; bin<=20; bin++)
         if (eff->GetBinContent(bin) > 0)
           Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
-
-      effError = (TH1*) eff->Clone("effError");
-      effError->Reset();
-
-      for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
-        if (eff->GetBinContent(bin) > 0)
-          effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
-
-      effError->SetLineColor(1);
-      effError->SetMarkerStyle(1);
-      effError->DrawCopy("SAME HIST");
+      
+      if (uncertainty) {
+       effError = (TH1*) eff->Clone("effError");
+       effError->Reset();
+
+       for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
+         if (eff->GetBinContent(bin) > 0)
+           effError->SetBinContent(bin, 10.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
+
+       effError->SetLineColor(1);
+       effError->SetMarkerStyle(1);
+       effError->DrawCopy("SAME HIST");
+      }
     }
 
     eff->SetBinContent(1, 0);
@@ -1726,7 +1732,8 @@ void EfficiencyComparison(Int_t eventType = 2)
     legend->AddEntry(eff, (((i == 0) ? "non diffractive" : ((i == 1) ? "single diffractive" : ((i == 2) ? "double diffractive" : "Pythia combined")))));
   }
 
-  legend->AddEntry(effError, "relative syst. uncertainty #times 10");
+  if (uncertainty)
+    legend->AddEntry(effError, "relative syst. uncertainty #times 10");
 
   legend->Draw();
 
@@ -2532,7 +2539,7 @@ TH1* SystematicsSummary(Bool_t tpc = 1)
   return total;
 }
 
-void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
+void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE, Bool_t small = kFALSE)
 {
   gSystem->Load("libPWG0base");
 
@@ -2584,7 +2591,7 @@ void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
   // change error drawing style
   systError->SetFillColor(15);
 
-  TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", 800, 400);
+  TCanvas* canvas = new TCanvas("finalPlot.eps", "finalPlot.eps", (small) ? 600 : 800, 400);
   canvas->SetRightMargin(0.05);
   canvas->SetTopMargin(0.05);
 
@@ -2647,3 +2654,161 @@ void finalPlot(Bool_t tpc = kTRUE, Bool_t chi2 = kTRUE)
 
   canvas->SaveAs(canvas->GetName());
 }
+
+void BlobelUnfoldingExample()
+{
+  const Int_t kSize = 20;
+
+  TMatrixD matrix(kSize, kSize);
+  for (Int_t x=0; x<kSize; x++)
+  {
+    for (Int_t y=0; y<kSize; y++)
+    {
+      if (x == y)
+      {
+        if (x == 0 || x == kSize -1)
+        {
+          matrix(x, y) = 0.75;
+        }
+        else
+          matrix(x, y) = 0.5;
+      }
+      else if (TMath::Abs(x - y) == 1)
+      {
+        matrix(x, y) = 0.25;
+      }
+    }
+  }
+
+  //matrix.Print();
+
+  TMatrixD inverted(matrix);
+  inverted.Invert();
+
+  //inverted.Print();
+
+  TH1F* inputDist = new TH1F("inputDist", ";t;#tilde{T}(t)", kSize, -0.5, (Float_t) kSize - 0.5);
+  TVectorD inputDistVector(kSize);
+  TH1F* unfolded = inputDist->Clone("unfolded");
+  TH1F* measuredIdealDist = inputDist->Clone("measuredIdealDist");
+  measuredIdealDist->SetTitle(";m;#tilde{M}(m)");
+  TH1F* measuredDist = measuredIdealDist->Clone("measuredDist");
+
+  TF1* gaus = new TF1("func", "gaus(0)", -0.5, kSize);
+  // norm: 1/(sqrt(2pi)sigma)
+  gaus->SetParameters(10000 / sqrt(2 * TMath::Pi()) / ((Float_t) kSize / 8), (Float_t) kSize / 2, (Float_t) kSize / 8);
+  //gaus->Print();
+
+  for (Int_t x=1; x<=inputDist->GetNbinsX(); x++)
+  {
+    Float_t value = gaus->Eval(inputDist->GetBinCenter(x));
+    inputDist->SetBinContent(x, value);
+    inputDistVector(x-1) = value;
+  }
+
+  TVectorD measuredDistIdealVector = matrix * inputDistVector;
+  
+  for (Int_t x=1; x<=measuredIdealDist->GetNbinsX(); x++)
+    measuredIdealDist->SetBinContent(x, measuredDistIdealVector(x-1));
+
+  measuredDist->FillRandom(measuredIdealDist, 10000);
+
+  // fill error matrix before scaling
+  TMatrixD covarianceMatrixMeasured(kSize, kSize);
+  for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
+    covarianceMatrixMeasured(x-1, x-1) = TMath::Sqrt(measuredDist->GetBinContent(x));
+
+  TMatrixD covarianceMatrix = inverted * covarianceMatrixMeasured * inverted;
+  //covarianceMatrix.Print();
+
+  TVectorD measuredDistVector(kSize);
+  for (Int_t x=1; x<=measuredDist->GetNbinsX(); x++)
+    measuredDistVector(x-1) = measuredDist->GetBinContent(x);
+
+  TVectorD unfoldedVector = inverted * measuredDistVector;
+  for (Int_t x=1; x<=unfolded->GetNbinsX(); x++)
+    unfolded->SetBinContent(x, unfoldedVector(x-1));
+
+  TCanvas* canvas = new TCanvas("BlobelUnfoldingExample", "BlobelUnfoldingExample", 1000, 500);
+  canvas->SetTopMargin(0.05);
+  canvas->Divide(2, 1);
+
+  canvas->cd(1);
+  canvas->cd(1)->SetLeftMargin(0.15);
+  canvas->cd(1)->SetRightMargin(0.05);
+  measuredDist->GetYaxis()->SetTitleOffset(1.7);
+  measuredDist->SetStats(0);
+  measuredDist->DrawCopy();
+  gaus->Draw("SAME");
+
+  canvas->cd(2);
+  canvas->cd(2)->SetLeftMargin(0.15);
+  canvas->cd(2)->SetRightMargin(0.05);
+  unfolded->GetYaxis()->SetTitleOffset(1.7);
+  unfolded->SetStats(0);
+  unfolded->DrawCopy();
+  gaus->Draw("SAME");
+
+  canvas->SaveAs("BlobelUnfoldingExample.eps");
+}
+
+void E735Fit()
+{
+  TH1* fCurrentESD = new TH1F("mult", "mult", 501, -0.5, 500.5);
+  fCurrentESD->Sumw2();
+
+  // Open the input stream
+  ifstream in;
+  in.open("e735data.txt");
+
+  while(in.good())
+  {
+    Float_t x, y, ye;
+    in >> x >> y >> ye;
+
+    //Printf("%f %f %f", x, y, ye);
+    fCurrentESD->SetBinContent(fCurrentESD->FindBin(x), y);
+    fCurrentESD->SetBinError(fCurrentESD->FindBin(x), ye);
+  }
+
+  in.close();
+
+  //new TCanvas; fCurrentESD->DrawCopy(); gPad->SetLogy();
+
+  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+
+  TF1* func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])");
+  func->SetParNames("scaling", "averagen", "k");
+  func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
+  func->SetParLimits(1, 0.001, 1000);
+  func->SetParLimits(2, 0.001, 1000);
+  func->SetParameters(fCurrentESD->GetMaximum() * 100, 10, 2);
+
+  TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
+  lognormal->SetParNames("scaling", "mean", "sigma");
+  lognormal->SetParameters(1, 1, 1);
+  lognormal->SetParLimits(0, 0, 10);
+  lognormal->SetParLimits(1, 0, 100);
+  lognormal->SetParLimits(2, 1e-3, 10);
+
+  TCanvas* canvas = new TCanvas("c1", "c1", 700, 400);
+  fCurrentESD->SetStats(kFALSE);
+  fCurrentESD->GetYaxis()->SetTitleOffset(1.3);
+  fCurrentESD->SetTitle(";true multiplicity (N);P_{N}");
+  fCurrentESD->Draw("");
+  fCurrentESD->GetXaxis()->SetRangeUser(0, 250);
+  fCurrentESD->Fit(func, "0", "", 0, 150);
+  func->SetRange(0, 250);
+  func->Draw("SAME");
+  printf("chi2 = %f\n", func->GetChisquare());
+
+  fCurrentESD->Fit(lognormal, "0", "", 0.01, 150);
+  lognormal->SetLineColor(2);
+  lognormal->SetLineStyle(2);
+  lognormal->SetRange(0, 250);
+  lognormal->Draw("SAME");
+
+  gPad->SetLogy();
+
+  canvas->SaveAs("E735Fit.eps");
+}
diff --git a/PWG0/dNdEta/run.C b/PWG0/dNdEta/run.C
new file mode 100644 (file)
index 0000000..5c7e744
--- /dev/null
@@ -0,0 +1,90 @@
+void run(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, Bool_t mc = kTRUE, const char* option = "")
+{
+  if (aProof)
+  {
+    TProof::Open("lxb6046");
+
+    // Enable the needed package
+    gProof->UploadPackage("STEERBase");
+    gProof->EnablePackage("STEERBase");
+    gProof->UploadPackage("ESD");
+    gProof->EnablePackage("ESD");
+    gProof->UploadPackage("ANALYSIS");
+    gProof->EnablePackage("ANALYSIS");
+    gProof->UploadPackage("PWG0base");
+    gProof->EnablePackage("PWG0base");
+  }
+  else
+  {
+    gSystem->Load("libVMC");
+    gSystem->Load("libTree");
+    gSystem->Load("libSTEERBase");
+    gSystem->Load("libESD");
+    gSystem->Load("libANALYSIS");
+    gSystem->Load("libPWG0base");
+  }
+
+  // Create chain of input files
+  gROOT->LoadMacro("../CreateESDChain.C");
+  chain = CreateESDChain(data, nRuns, offset);
+
+  // Create the analysis manager
+  mgr = new AliAnalysisManager;
+
+  // selection of esd tracks
+  gROOT->ProcessLine(".L CreateCuts.C");
+  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
+  if (!esdTrackCuts)
+  {
+    printf("ERROR: esdTrackCuts could not be created\n");
+    return;
+  }
+
+  TString taskName("AlidNdEtaTask.cxx+");
+  if (aDebug)
+    taskName += "+g";
+
+  // Create, add task
+  if (aProof) {
+    gProof->Load(taskName);
+  } else
+    gROOT->Macro(taskName);
+
+  task = new AlidNdEtaTask(option);
+  task->SetTrackCuts(esdTrackCuts);
+  task->SetAnalysisMode(AlidNdEtaTask::kSPD);
+
+  if (mc)
+    task->SetReadMC();
+
+  mgr->AddTask(task);
+
+  if (mc) {
+    // Enable MC event handler
+    AliMCEventHandler* handler = new AliMCEventHandler;
+    handler->SetReadTR(kFALSE);
+    mgr->SetMCtruthEventHandler(handler);
+  }
+
+  // Add ESD handler
+  AliESDInputHandler* esdH = new AliESDInputHandler;
+  mgr->SetInputEventHandler(esdH);
+
+  // Attach input
+  cInput  = mgr->CreateContainer("cInput", TChain::Class(), AliAnalysisManager::kInputContainer);
+  mgr->ConnectInput(task, 0, cInput);
+
+  // Attach output
+  cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
+  mgr->ConnectOutput(task, 0, cOutput);
+
+  // Enable debug printouts
+  if (aDebug)
+    mgr->SetDebugLevel(2);
+
+  // Run analysis
+  mgr->InitAnalysis();
+  mgr->PrintStatus();
+
+  mgr->StartAnalysis((aProof) ? "proof" : "local", chain);
+}
diff --git a/PWG0/dNdEta/runCorrection.C b/PWG0/dNdEta/runCorrection.C
new file mode 100644 (file)
index 0000000..94f2cfa
--- /dev/null
@@ -0,0 +1,86 @@
+void runCorrection(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = ""){
+  if (aProof)
+  {
+    TProof::Open("lxb6046");
+
+    // Enable the needed package
+    gProof->UploadPackage("STEERBase");
+    gProof->EnablePackage("STEERBase");
+    gProof->UploadPackage("ESD");
+    gProof->EnablePackage("ESD");
+    gProof->UploadPackage("ANALYSIS");
+    gProof->EnablePackage("ANALYSIS");
+    gProof->UploadPackage("PWG0base");
+    gProof->EnablePackage("PWG0base");
+  }
+  else
+  {
+    gSystem->Load("libVMC");
+    gSystem->Load("libTree");
+    gSystem->Load("libSTEERBase");
+    gSystem->Load("libESD");
+    gSystem->Load("libANALYSIS");
+    gSystem->Load("libPWG0base");
+  }
+
+  // Create chain of input files
+  gROOT->LoadMacro("../CreateESDChain.C");
+  chain = CreateESDChain(data, nRuns, offset);
+
+  // Create the analysis manager
+  mgr = new AliAnalysisManager;
+
+  // selection of esd tracks
+  gROOT->ProcessLine(".L CreateCuts.C");
+  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
+  if (!esdTrackCuts)
+  {
+    printf("ERROR: esdTrackCuts could not be created\n");
+    return;
+  }
+
+  TString taskName("AlidNdEtaCorrectionTask.cxx+");
+  if (aDebug)
+    taskName += "+g";
+
+  // Create, add task
+  if (aProof) {
+    gProof->Load(taskName);
+  } else
+    gROOT->Macro(taskName);
+
+  task = new AlidNdEtaCorrectionTask(option);
+  task->SetTrackCuts(esdTrackCuts);
+  task->SetAnalysisMode(AlidNdEtaCorrectionTask::kSPD);
+
+  mgr->AddTask(task);
+
+  // Enable MC event handler
+  AliMCEventHandler* handler = new AliMCEventHandler;
+  handler->SetReadTR(kFALSE);
+  mgr->SetMCtruthEventHandler(handler);
+
+  // Add ESD handler
+  AliESDInputHandler* esdH = new AliESDInputHandler;
+  mgr->SetInputEventHandler(esdH);
+
+  // Attach input
+  cInput  = mgr->CreateContainer("cInput", TChain::Class(), AliAnalysisManager::kInputContainer);
+  mgr->ConnectInput(task, 0, cInput);
+
+  // Attach output
+  cOutput = mgr->CreateContainer("cOutput", TList::Class(), AliAnalysisManager::kOutputContainer);
+  mgr->ConnectOutput(task, 0, cOutput);
+
+  // Enable debug printouts
+  if (aDebug)
+    mgr->SetDebugLevel(2);
+
+  // Run analysis
+  mgr->InitAnalysis();
+  mgr->PrintStatus();
+
+  //aProof = kFALSE;
+
+  mgr->StartAnalysis((aProof) ? "proof" : "local", chain);
+}
index d9e6fed..a52d2b7 100644 (file)
@@ -125,10 +125,20 @@ void draw(const char* fileName = "multiplicityMC.root", const char* folder = "Mu
   canvas->SaveAs("Plot_Correlation.pdf");
 }
 
-void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
+void loadlibs()
 {
+  gSystem->Load("libTree");
+  gSystem->Load("libVMC");
+
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libANALYSIS");
   gSystem->Load("libPWG0base");
+}
 
+void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE, Float_t beta  = 1e4)
+{
+  loadlibs();
+  
   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
 
   TFile::Open(fileNameMC);
@@ -160,12 +170,14 @@ void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fol
 
   if (chi2)
   {
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
-    //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e3);
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 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, 1e5);
     //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
     //mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kTRUE, hist2->ProjectionY("mymchist"));
-    mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, hist2->ProjectionY("mymchist"));
+    mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); //hist2->ProjectionY("mymchist"));
     mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
   }
   else
@@ -247,168 +259,6 @@ const char* GetEventTypeName(Int_t type)
   return 0;
 }
 
-/*void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
-{
-  gSystem->mkdir(targetDir);
-
-  gSystem->Load("libPWG0base");
-
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  TFile::Open(fileNameMC);
-  mult->LoadHistograms("Multiplicity");
-
-  AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
-  TFile::Open(fileNameESD);
-  multESD->LoadHistograms("Multiplicity");
-  mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
-
-  TCanvas* canvas = new TCanvas("EvaluateBayesianMethod", "EvaluateBayesianMethod", 800, 600);
-  TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
-  legend->SetFillColor(0);
-
-  Float_t min = 1e10;
-  Float_t max = 0;
-
-  TGraph* first = 0;
-  Int_t count = 0; // just to order the saved images...
-
-  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
-  {
-    TGraph* fitResultsMC = new TGraph;
-    fitResultsMC->SetTitle(";Weight Parameter");
-    TGraph* fitResultsRes = new TGraph;
-    fitResultsRes->SetTitle(";Weight Parameter");
-
-    fitResultsMC->SetFillColor(0);
-    fitResultsRes->SetFillColor(0);
-    fitResultsMC->SetMarkerStyle(20+type);
-    fitResultsRes->SetMarkerStyle(24+type);
-    fitResultsRes->SetMarkerColor(kRed);
-    fitResultsRes->SetLineColor(kRed);
-
-    legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
-    legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
-
-    for (Float_t weight = 0; weight < 1.01; weight += 0.1)
-    {
-      Float_t chi2MC = 0;
-      Float_t residuals = 0;
-
-      mult->ApplyBayesianMethod(histID, kFALSE, type, weight, 100, 0, kFALSE);
-      mult->DrawComparison(Form("%s/Bayesian_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
-      mult->GetComparisonResults(&chi2MC, 0, &residuals);
-
-      fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
-      fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
-
-      min = TMath::Min(min, TMath::Min(chi2MC, residuals));
-      max = TMath::Max(max, TMath::Max(chi2MC, residuals));
-    }
-
-    fitResultsMC->Print();
-    fitResultsRes->Print();
-
-    canvas->cd();
-    fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME"));
-    fitResultsRes->Draw("SAME CP");
-
-    if (first == 0)
-      first = fitResultsMC;
-  }
-
-  gPad->SetLogy();
-  printf("min = %f, max = %f\n", min, max);
-  if (min <= 0)
-    min = 1e-5;
-  first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
-
-  legend->Draw();
-
-  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
-  canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
-}
-
-void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
-{
-  gSystem->mkdir(targetDir);
-
-  gSystem->Load("libPWG0base");
-
-  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
-  TFile::Open(fileNameMC);
-  mult->LoadHistograms("Multiplicity");
-
-  AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
-  TFile::Open(fileNameESD);
-  multESD->LoadHistograms("Multiplicity");
-  mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
-
-  TCanvas* canvas = new TCanvas("EvaluateBayesianMethodIterations", "EvaluateBayesianMethodIterations", 800, 600);
-  TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
-  legend->SetFillColor(0);
-
-  Float_t min = 1e10;
-  Float_t max = 0;
-
-  TGraph* first = 0;
-  Int_t count = 0; // just to order the saved images...
-
-  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
-  {
-    TGraph* fitResultsMC = new TGraph;
-    fitResultsMC->SetTitle(";Iterations");
-    TGraph* fitResultsRes = new TGraph;
-    fitResultsRes->SetTitle(";Iterations");
-
-    fitResultsMC->SetFillColor(0);
-    fitResultsRes->SetFillColor(0);
-    fitResultsMC->SetMarkerStyle(20+type);
-    fitResultsRes->SetMarkerStyle(24+type);
-    fitResultsRes->SetMarkerColor(kRed);
-    fitResultsRes->SetLineColor(kRed);
-
-    legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
-    legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
-
-    for (Int_t iter = 5; iter <= 50; iter += 5)
-    {
-      Float_t chi2MC = 0;
-      Float_t residuals = 0;
-
-      mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1, iter);
-      mult->DrawComparison(Form("%s/BayesianIter_%02d_%d_%d", targetDir, count++, type, iter), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
-      mult->GetComparisonResults(&chi2MC, 0, &residuals);
-
-      fitResultsMC->SetPoint(fitResultsMC->GetN(), iter, chi2MC);
-      fitResultsRes->SetPoint(fitResultsRes->GetN(), iter, residuals);
-
-      min = TMath::Min(min, TMath::Min(chi2MC, residuals));
-      max = TMath::Max(max, TMath::Max(chi2MC, residuals));
-    }
-
-    fitResultsMC->Print();
-    fitResultsRes->Print();
-
-    canvas->cd();
-    fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME"));
-    fitResultsRes->Draw("SAME CP");
-
-    if (first == 0)
-      first = fitResultsMC;
-  }
-
-  gPad->SetLogy();
-  printf("min = %f, max = %f\n", min, max);
-  if (min <= 0)
-    min = 1e-5;
-  first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
-
-  legend->Draw();
-
-  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
-  canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
-}*/
-
 void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
 {
   gSystem->mkdir(targetDir);
@@ -429,7 +279,7 @@ void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multipl
   TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE");
 
   Int_t colors[3] = {1, 2, 4};
-  Int_t markers[12] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3};
+  Int_t markers[20] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3, 4, 5, 6};
 
   for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
   //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
@@ -439,9 +289,9 @@ void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multipl
 
     TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
 
-    for (Int_t i = 1; i <= 4; i++)
+    for (Int_t i = 1; i <= 5; i++)
     {
-      Int_t iterArray[4] = {5, 20, 50, 100};
+      Int_t iterArray[5] = {5, 20, 50, 100, -1};
       //Int_t iter = i * 40 - 20;
       Int_t iter = iterArray[i-1];
 
@@ -454,7 +304,8 @@ void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multipl
         fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
         fitResultsMC[region]->SetName(Form("%s_MC_%d", tmp.Data(), i * AliMultiplicityCorrection::kQualityRegions + region - 2));
         fitResultsMC[region]->SetFillColor(0);
-        fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]);
+        //fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]);
+        fitResultsMC[region]->SetMarkerStyle(markers[(i-1)]);
         fitResultsMC[region]->SetLineColor(colors[region]);
       }
 
@@ -639,6 +490,154 @@ void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes =
   canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
 }
 
+void EvaluateDrawResultRegions(const char* targetDir, Int_t type = 0)
+{
+  gSystem->Load("libPWG0base");
+
+  TString name;
+  TFile* graphFile = 0;
+  if (type == -1)
+  {
+    name = "EvaluateChi2Method";
+    graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
+  }
+  else
+  {
+    name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
+    graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
+  }
+
+  TCanvas* canvas = new TCanvas(name, name, 800, 1200);
+  //canvas->Divide(1, AliMultiplicityCorrection::kQualityRegions, 0, 0);
+  canvas->Range(0, 0, 1, 1);
+
+  TPad* pad[3];
+  pad[0] = new TPad(Form("%s_pad1", name.Data()), "", 0.02, 0.05, 0.98, 0.35);
+  pad[1] = new TPad(Form("%s_pad2", name.Data()), "", 0.02, 0.35, 0.98, 0.65);
+  pad[2] = new TPad(Form("%s_pad3", name.Data()), "", 0.02, 0.65, 0.98, 0.95);
+
+  Float_t yMinRegion[3];
+  for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
+    yMinRegion[i] = 1e20;
+
+  for (Int_t region = 1; region <= AliMultiplicityCorrection::kQualityRegions; region++)
+  {
+    canvas->cd();
+    pad[region-1]->Draw();
+    pad[region-1]->cd();
+    pad[region-1]->SetRightMargin(0.05);
+
+    if (region != 1)
+      pad[region-1]->SetBottomMargin(0);
+    if (region != AliMultiplicityCorrection::kQualityRegions)
+      pad[region-1]->SetTopMargin(0);
+
+    if (type == -1)
+    {
+      pad[region-1]->SetLogx();
+      pad[region-1]->SetLogy();
+    }
+    pad[region-1]->SetTopMargin(0.05);
+    pad[region-1]->SetGridx();
+    pad[region-1]->SetGridy();
+
+    TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98);
+    legend->SetFillColor(0);
+
+    Int_t count = 0;
+
+    Float_t xMin = 1e20;
+    Float_t xMax = 0;
+
+    Float_t yMin = 1e20;
+    Float_t yMax = 0;
+
+    TString xaxis, yaxis;
+
+    while (1)
+    {
+      count++;
+
+      TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
+      if (!mc)
+        break;
+
+      if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
+        continue;
+
+      xaxis = mc->GetXaxis()->GetTitle();
+      yaxis = mc->GetYaxis()->GetTitle();
+
+      mc->Print();
+
+      xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
+      yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
+
+      xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
+      yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
+
+      for (Int_t i=0; i<mc->GetN(); ++i)
+        yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
+    }
+
+    if (type >= 0)
+    {
+      xaxis = "smoothing parameter";
+    }
+    else if (type == -1)
+    {
+      xaxis = "weight parameter";
+      xMax *= 5;
+    }
+    yaxis = "P_{1}"; // (2 <= t <= 150)";
+
+    printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
+
+    TGraph* dummy = new TGraph;
+    dummy->SetPoint(0, xMin, yMin);
+    dummy->SetPoint(1, xMax, yMax);
+    dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
+
+    dummy->SetMarkerColor(0);
+    dummy->Draw("AP");
+    dummy->GetYaxis()->SetMoreLogLabels(1);
+
+    count = 0;
+
+    while (1)
+    {
+      count++;
+
+      TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
+      if (!mc)
+        break;
+      if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
+        continue;
+
+      printf("Loaded %d sucessful.\n", count);
+
+      if (type == -1)
+      {
+        legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
+      }
+      else
+        legend->AddEntry(mc);
+
+      mc->Draw("SAME PC");
+
+     }
+
+     legend->Draw();
+  }
+
+  for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
+    Printf("Minimum for region %d is %f", i, yMinRegion[i]);
+
+  canvas->Modified();
+  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
+  canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
+}
+
 void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", const char* targetDir, Int_t histID = 3)
 {
   gSystem->mkdir(targetDir);
@@ -664,9 +663,9 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC_2M.root", const
 
   TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
 
-  for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kLog; ++type)
+  for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
 //  for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kEntropy; type <= AliMultiplicityCorrection::kEntropy; ++type)
-//  for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kPol0; ++type)
+  //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
   {
     TGraph* fitResultsMC[3];
     for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
index 10d23b7..b653112 100644 (file)
@@ -70,9 +70,19 @@ void rundNdEtaAnalysis(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC
   }
 }
 
-void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const char* dataOutput = "analysis_esd.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
+void loadlibs()
 {
+  gSystem->Load("libTree");
+  gSystem->Load("libVMC");
+
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libANALYSIS");
   gSystem->Load("libPWG0base");
+}
+
+void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const char* dataOutput = "analysis_esd.root", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
+{
+  loadlibs();
 
   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
   TFile::Open(correctionMapFile);
@@ -87,7 +97,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
   }
 
   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
-  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
   fdNdEtaAnalysis->DrawHistograms(kTRUE);
   TFile* file2 = TFile::Open(dataOutput, "RECREATE");
@@ -95,7 +105,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
 
   file->cd();
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
-  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco);
   fdNdEtaAnalysis->DrawHistograms(kTRUE);
   file2->cd();
@@ -103,7 +113,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
 
   file->cd();
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
-  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
   fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle);
   fdNdEtaAnalysis->DrawHistograms(kTRUE);
   file2->cd();
@@ -111,7 +121,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
 
   file->cd();
   fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
-  fdNdEtaAnalysis->LoadHistograms("dndeta");
+  fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
   fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
   fdNdEtaAnalysis->DrawHistograms(kTRUE);
   file2->cd();
@@ -120,7 +130,7 @@ void FinishAnalysisAll(const char* dataInput = "analysis_esd_raw.root", const ch
 
 void FinishAnalysis(const char* analysisFile = "analysis_esd.root", const char* analysisDir = "dndeta", const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", Bool_t useUncorrected = kFALSE, Bool_t simple = kFALSE)
 {
-  gSystem->Load("libPWG0base");
+  loadlibs();
 
   TFile* file = TFile::Open(analysisFile);
 
@@ -133,7 +143,9 @@ void FinishAnalysis(const char* analysisFile = "analysis_esd.root", const char*
     TFile::Open(correctionMapFile);
     dNdEtaCorrection->LoadHistograms();
 
-    fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3);
+    //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.2, AlidNdEtaCorrection::kINEL);
+    fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kINEL);
+    //fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
   }
 
   fdNdEtaAnalysis->DrawHistograms(simple);