PWGUD/base/dNdEtaAnalysis.cxx

   1 /* $Id$ */
   2
   3 #include "dNdEtaAnalysis.h"
   4
   5 #include <TFile.h>
   6 #include <TH3F.h>
   7 #include <TH2F.h>
   8 #include <TH1F.h>
   9 #include <TMath.h>
  10 #include <TCanvas.h>
  11 #include <TCollection.h>
  12 #include <TIterator.h>
  13 #include <TList.h>
  14 #include <TLegend.h>
  15 #include <TLine.h>
  16 #include <TParameter.h>
  17
  18 #include <AliLog.h>
  19 #include "AlidNdEtaCorrection.h"
  20 #include <AliCorrection.h>
  21 #include <AliPWG0Helper.h>
  22 #include <AliCorrectionMatrix2D.h>
  23 #include <AliCorrectionMatrix3D.h>
  24
  25 //____________________________________________________________________
  26 ClassImp(dNdEtaAnalysis)
  27
  28 //____________________________________________________________________
  29 dNdEtaAnalysis::dNdEtaAnalysis() :
  30   TNamed(),
  31   fData(0),
  32   fMult(0),
  33   fPtDist(0),
  34   fAnalysisMode(AliPWG0Helper::kInvalid),
  35   fTag(),
  36   fvtxMin(-9.99),
  37   fvtxMax(9.99)
  38 {
  39   // default constructor
  40
  41   for (Int_t i=0; i<kVertexBinning; ++i)
  42   {
  43     fdNdEta[i] = 0;
  44     fdNdEtaPtCutOffCorrected[i] = 0;
  45   }
  46 }
  47
  48 //____________________________________________________________________
  49 dNdEtaAnalysis::dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
  50   TNamed(name, title),
  51   fData(0),
  52   fMult(0),
  53   fPtDist(0),
  54   fAnalysisMode(analysisMode),
  55   fTag(),
  56   fvtxMin(-9.99),
  57   fvtxMax(9.99)
  58 {
  59   // constructor
  60
  61   fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode);
  62
  63   // do not add this hists to the directory
  64   Bool_t oldStatus = TH1::AddDirectoryStatus();
  65   TH1::AddDirectory(kFALSE);
  66
  67   fMult = new TH1F("TriggeredMultiplicity", "Triggered Events;raw multiplicity;entries", 1000, -0.5, 999.5);
  68
  69   TH1* histForBinning = fData->GetTrackCorrection()->GetGeneratedHistogram();
  70   fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", histForBinning->GetNbinsY(), histForBinning->GetYaxis()->GetXbins()->GetArray());
  71
  72   fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
  73
  74   for (Int_t i=1; i<kVertexBinning; ++i)
  75   {
  76     fdNdEta[i]    = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
  77     fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
  78   }
  79
  80   fPtDist = new TH1F("Pt", "p_{T} distribution;p_{T} [GeV/c];#frac{dN}{d#eta dp_{T}} [c/GeV]", histForBinning->GetNbinsZ(), histForBinning->GetZaxis()->GetXbins()->GetArray());
  81
  82   TH1::AddDirectory(oldStatus);
  83 }
  84
  85 //____________________________________________________________________
  86 dNdEtaAnalysis::~dNdEtaAnalysis()
  87 {
  88   // destructor
  89
  90   if (fData)
  91   {
  92     delete fData;
  93     fData = 0;
  94   }
  95
  96   if (fMult)
  97   {
  98     delete fMult;
  99     fMult = 0;
 100   }
 101
 102   for (Int_t i=0; i<kVertexBinning; ++i)
 103   {
 104     if (fdNdEta[i])
 105     {
 106       delete fdNdEta[i];
 107       fdNdEta[i] = 0;
 108     }
 109     if (fdNdEtaPtCutOffCorrected[i])
 110     {
 111       delete fdNdEtaPtCutOffCorrected[i];
 112       fdNdEtaPtCutOffCorrected[i] = 0;
 113     }
 114   }
 115
 116   if (fPtDist)
 117   {
 118     delete fPtDist;
 119     fPtDist = 0;
 120   }
 121 }
 122
 123 //_____________________________________________________________________________
 124 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
 125   TNamed(c),
 126   fData(0),
 127   fMult(0),
 128   fPtDist(0),
 129   fAnalysisMode(AliPWG0Helper::kInvalid),
 130   fTag(),
 131   fvtxMin(-9.99),
 132   fvtxMax(9.99)
 133 {
 134   //
 135   // dNdEtaAnalysis copy constructor
 136   //
 137
 138   ((dNdEtaAnalysis &) c).Copy(*this);
 139 }
 140
 141 //_____________________________________________________________________________
 142 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
 143 {
 144   //
 145   // Assignment operator
 146   //
 147
 148   if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
 149   return *this;
 150 }
 151
 152 //_____________________________________________________________________________
 153 void dNdEtaAnalysis::Copy(TObject &c) const
 154 {
 155   //
 156   // Copy function
 157   //
 158
 159   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
 160
 161   target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
 162   target.fMult = dynamic_cast<TH1F*> (fMult->Clone());
 163
 164   for (Int_t i=0; i<kVertexBinning; ++i)
 165   {
 166     target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
 167     target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
 168   }
 169
 170   target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
 171
 172   target.fAnalysisMode = fAnalysisMode;
 173   target.fTag = fTag;
 174
 175   TNamed::Copy((TNamed &) c);
 176 }
 177
 178 //____________________________________________________________________
 179 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
 180 {
 181   // fills a track into the histograms
 182
 183   fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
 184 }
 185
 186 //____________________________________________________________________
 187 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
 188 {
 189   // fills an event into the histograms
 190
 191   fData->GetEventCorrection()->FillMeas(vtx, n);
 192 }
 193
 194 //____________________________________________________________________
 195 void dNdEtaAnalysis::FillTriggeredEvent(Float_t n)
 196 {
 197   // fills a triggered event into the histograms
 198
 199   fMult->Fill(n);
 200 }
 201
 202 //____________________________________________________________________
 203 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag, Int_t backgroundSubtraction, TH1* combinatoricsCorrection)
 204 {
 205   //
 206   // correct with the given correction values and calculate dNdEta and pT distribution
 207   // the corrections that are applied can be steered by the flag correctionType
 208   // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
 209   //
 210
 211   fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut);
 212   Printf("\n\n%s", fTag.Data());
 213
 214   if (combinatoricsCorrection)
 215     Printf("Combinatorics subtraction active!");
 216
 217   // set corrections to 1
 218   fData->SetCorrectionToUnity();
 219
 220   if (correction && correctionType != AlidNdEtaCorrection::kNone)
 221   {
 222     TH3* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
 223     TH2* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
 224
 225     if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
 226       trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
 227
 228     if (correctionType >= AlidNdEtaCorrection::kVertexReco)
 229     {
 230       trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
 231       eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
 232
 233       // set bin with multiplicity 0 to unity (correction has no meaning in this bin)
 234       for (Int_t i=0; i<=eventCorr->GetNbinsX()+1; i++)
 235         eventCorr->SetBinContent(i, 1, 1);
 236     }
 237
 238     switch (correctionType)
 239     {
 240       case AlidNdEtaCorrection::kINEL :
 241       {
 242         trackCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetCorrectionHistogram());
 243         eventCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram());
 244         break;
 245       }
 246       case AlidNdEtaCorrection::kNSD :
 247       {
 248         trackCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetTrackCorrection()->GetCorrectionHistogram());
 249         eventCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetEventCorrection()->GetCorrectionHistogram());
 250         break;
 251       }
 252       case AlidNdEtaCorrection::kND :
 253       {
 254         trackCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetTrackCorrection()->GetCorrectionHistogram());
 255         eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram());
 256         break;
 257       }
 258       case AlidNdEtaCorrection::kOnePart :
 259       {
 260         trackCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetTrackCorrection()->GetCorrectionHistogram());
 261         eventCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->GetCorrectionHistogram());
 262         break;
 263       }
 264       default : break;
 265     }
 266   }
 267   else
 268     printf("INFO: No correction applied\n");
 269
 270   TH2F* rawMeasured = (TH2F*) fData->GetEventCorrection()->GetMeasuredHistogram()->Clone("rawMeasured");
 271
 272   fData->ResetErrorsOnCorrections();
 273   fData->Multiply();
 274
 275   if (correction && correctionType >= AlidNdEtaCorrection::kVertexReco)
 276   {
 277     // There are no events with vertex that have 0 multiplicity, therefore
 278     //   populate bin with 0 multiplicity with the following idea:
 279     //     alpha = triggered events with vertex at a given vertex position / all triggered events with vertex
 280     //     triggered events without vertex and 0 multiplicity at a given vertex position = alpha * all triggered events with 0 multiplicity
 281     //   afterwards we still correct for the trigger efficiency
 282     // at the same time calculate expectation from MC (not used, just to check the value)
 283
 284     //TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
 285     TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
 286
 287     TH2* eAll  =    correction->GetCorrection(correctionType)->GetEventCorrection()->GetGeneratedHistogram();
 288     TH2* eTrig =    correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
 289     TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
 290     TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsY()+1);
 291
 292     //new TCanvas; correctedEvents->DrawCopy("TEXT");
 293
 294     // start above 0 mult. bin with integration
 295     TH1* vertexDist = rawMeasured->ProjectionX("vertexdist_measured", 2, rawMeasured->GetNbinsY()+1);
 296     //new TCanvas; vertexDist->DrawCopy();
 297
 298     Int_t allEventsWithVertex = (Int_t) vertexDist->Integral(0, vertexDist->GetNbinsX()+1); // include under/overflow!
 299     Int_t triggeredEventsWith0Mult = (Int_t) fMult->GetBinContent(1);
 300
 301     Printf("%d triggered events with 0 mult. -- %d triggered events with vertex", triggeredEventsWith0Mult, allEventsWithVertex);
 302
 303     if (backgroundSubtraction > 0)
 304     {
 305       triggeredEventsWith0Mult -= backgroundSubtraction;
 306       Printf("Subtracted %d background events from 0 mult. bin", backgroundSubtraction);
 307     }
 308
 309     TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
 310     kineBias->Reset();
 311
 312     // loop over vertex bins
 313     for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
 314     {
 315       Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
 316       Double_t events = alpha * triggeredEventsWith0Mult;
 317
 318       if (eTrigVtx_projx->GetBinContent(i) == 0)
 319         continue;
 320
 321       // calculate how many events we would have got with a pure MC-based correction
 322       //   in the given bin: measured events with vertex (mult > 0) * triggered events with mult 0 (mc) / events with vertex and mult > 0 (mc) * trigger correction for bin 0
 323       Printf("+++ 0-Bin Correction for bin %d +++", i);
 324       Printf("  Events: %f", vertexDist->GetBinContent(i));
 325       Printf("  Ratio triggered N==0 / triggered vertex N>0: %f", eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i));
 326       Printf("  Ratio all N==0 / triggered vertex N>0: %f", eAll->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i));
 327       Printf("  Correction factor: %f", fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1));
 328
 329       //Double_t mcEvents = vertexDist->GetBinContent(i) * eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i) * fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
 330       Double_t mcEvents = vertexDist->GetBinContent(i) * eAll->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i);
 331       if (backgroundSubtraction == -1)
 332       {
 333         Printf("Using MC value for 0-bin correction!");
 334         events = mcEvents;
 335       }
 336       else
 337       {
 338         Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) *
 339           eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1);
 340
 341         kineBias->SetBinContent(i, fZ);
 342
 343         events *= fZ;
 344
 345         // multiply with trigger correction if set above
 346         events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
 347
 348         Printf("  Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, mcEvents);
 349       }
 350
 351       correctedEvents->SetBinContent(i, 1, events);
 352     }
 353
 354     Printf("In |vtx-z| < 10 cm: %d events have been added", (Int_t) correctedEvents->Integral(vertexDist->FindBin(-9.9), vertexDist->FindBin(9.9), 1, 1));
 355
 356     //new TCanvas; correctedEvents->DrawCopy("TEXT");
 357     //new TCanvas; kineBias->DrawCopy();
 358   }
 359
 360   fData->PrintInfo(ptCut);
 361
 362   TH3* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
 363
 364   // integrate multiplicity axis out (include under/overflow bins!!!)
 365   TH2* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
 366
 367   TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
 368
 369   // create pt hist
 370   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 371   {
 372     // reset all ranges
 373     dataHist->GetXaxis()->SetRange(0, 0);
 374     dataHist->GetYaxis()->SetRange(0, 0);
 375     dataHist->GetZaxis()->SetRange(0, 0);
 376
 377     // vtx cut
 378     Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5);
 379     Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5);
 380     dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
 381     Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
 382
 383     if (nEvents > 0)
 384     {
 385       // eta cut
 386       dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8));
 387       Float_t etaWidth = 1.6;
 388
 389       TH1D* ptHist = static_cast<TH1D*> (dataHist->Project3D("ze"));
 390
 391       for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
 392       {
 393         Float_t binSize = fPtDist->GetBinWidth(i);
 394         fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
 395         fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
 396       }
 397
 398       delete ptHist;
 399     }
 400     else
 401       printf("ERROR: nEvents is 0!\n");
 402   }
 403
 404   // reset all ranges
 405   dataHist->GetXaxis()->SetRange(0, 0);
 406   dataHist->GetYaxis()->SetRange(0, 0);
 407   dataHist->GetZaxis()->SetRange(0, 0);
 408
 409   // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD)
 410   Int_t ptLowBin = 1;
 411   if (ptCut > 0 && (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS))
 412     ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
 413
 414   //new TCanvas; dataHist->DrawCopy();
 415
 416   //dataHist->Sumw2();
 417   dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
 418   printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
 419   TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
 420
 421   //new TCanvas; vtxVsEta->Draw("COLZ");
 422
 423   dataHist->GetZaxis()->SetRange(0, 0);
 424   vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
 425   vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
 426
 427   if (vtxVsEta == 0)
 428   {
 429     printf("ERROR: pt/multiplicity integration failed\n");
 430     return;
 431   }
 432
 433   //new TCanvas(tag, tag, 800, 600);  vtxVsEta->DrawCopy("COLZ");
 434
 435   // clear result histograms
 436   for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
 437   {
 438     fdNdEta[vertexPos]->Reset();
 439     fdNdEtaPtCutOffCorrected[vertexPos]->Reset();
 440   }
 441
 442   //const Float_t vertexRangeBegin[kVertexBinning] = { -9.99,  -9.99,  0.01 };
 443   //const Float_t vertexRangeEnd[kVertexBinning]   = {  9.99,  -0.01,  9.99 };
 444   const Float_t vertexRangeBegin[kVertexBinning] = { fvtxMin,  fvtxMin,  0.01 };
 445   const Float_t vertexRangeEnd[kVertexBinning]   = { fvtxMax,  -0.01,  fvtxMax };
 446
 447   for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
 448   {
 449     // loop over vertex ranges
 450     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
 451     {
 452       Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]);
 453       Int_t vertexBinEnd   = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
 454
 455       const Int_t *binBegin = 0;
 456       const Int_t maxBins = 20;
 457       const Int_t maxBins1 = 40;
 458
 459       if (dataHist->GetNbinsY() != maxBins)
 460         AliFatal(Form("Binning of acceptance is different from data histogram: data=%d, acceptance=%d", dataHist->GetNbinsY(), maxBins));
 461
 462       // adjust acceptance range
 463       // produce with drawPlots.C: DetermineAcceptance(...)
 464
 465       //const Int_t binBeginSPD[maxBins] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
 466       //const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1};
 467       //const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, -1, 14, 10, 8, 7, 6, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1};
 468
 469       const Int_t binBeginSPD[maxBins1] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4};
 470       const Int_t binBeginTPC[maxBins1] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1};
 471       const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1};
 472       if (fAnalysisMode & AliPWG0Helper::kSPD)
 473       {
 474         binBegin = binBeginSPD;
 475       }
 476       else if (fAnalysisMode & AliPWG0Helper::kTPC)
 477       {
 478         binBegin = binBeginTPC;
 479       }
 480       else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
 481       {
 482         binBegin = binBeginTPCITS;
 483       }
 484
 485       Int_t vtxBegin = 1;
 486       Int_t vtxEnd   = maxBins;
 487
 488       if (binBegin)
 489       {
 490         vtxBegin = binBegin[iEta - 1];
 491         vtxEnd = (Int_t) vtxVsEta->GetNbinsX() + 1 - binBegin[maxBins - iEta];
 492       }
 493       else
 494         Printf("WARNING: No acceptance applied!");
 495
 496       // eta range not accessible
 497       if (vtxBegin == -1)
 498         continue;
 499
 500       //Printf("%d %d | %d %d", vtxBegin, vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)), vtxEnd, vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta)));
 501       //vtxBegin = vertexHist->GetXaxis()->FindBin(GetVtxMin(eta));
 502       //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta));
 503
 504       //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
 505       //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
 506       vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
 507       vertexBinEnd =   TMath::Min(vertexBinEnd, vtxEnd);
 508       //Printf(" after:  %d %d", vertexBinBegin, vertexBinEnd);
 509
 510       // no data for this bin
 511       if (vertexBinBegin > vertexBinEnd)
 512       {
 513         //Printf("Bin empty. Skipped");
 514         continue;
 515       }
 516
 517       Float_t totalEvents = 0;
 518       Float_t sum = 0;
 519       Float_t sumError2 = 0;
 520       Float_t unusedTracks = 0;
 521       Float_t unusedEvents = 0;
 522       for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++)
 523       {
 524         if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd)
 525         {
 526           if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
 527           {
 528             sum += vtxVsEta->GetBinContent(iVtx, iEta);
 529
 530             if (sumError2 > 10e30)
 531               Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
 532
 533             sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
 534           }
 535           totalEvents += vertexHist->GetBinContent(iVtx);
 536         }
 537         else
 538         {
 539           unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta);
 540           unusedEvents += vertexHist->GetBinContent(iVtx);
 541         }
 542       }
 543
 544       if (totalEvents == 0)
 545       {
 546         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
 547         continue;
 548       }
 549
 550       Float_t ptCutOffCorrection = 1;
 551
 552       // find pt cut off correction factor
 553       if ((fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) && (fAnalysisMode & AliPWG0Helper::kFieldOn))
 554       {
 555         if (correction && ptCut > 0)
 556             ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd);
 557
 558         if (ptCutOffCorrection <= 0)
 559         {
 560             printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
 561             continue;
 562         }
 563       }
 564
 565       //printf("Eta: %d (%f) Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f \n", iEta, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection);
 566
 567       Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
 568       if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX())
 569       {
 570         Float_t dndeta = sum / totalEvents;
 571         Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
 572
 573         // correct for additional combinatorics
 574         Float_t combCorr = 0;
 575         if (combinatoricsCorrection)
 576         {
 577           combCorr = combinatoricsCorrection->GetBinContent(combinatoricsCorrection->GetXaxis()->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta)));
 578           dndeta *= combCorr;
 579           error *= combCorr;
 580         }
 581
 582         dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
 583         error  = error / fdNdEta[vertexPos]->GetBinWidth(bin);
 584
 585         fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
 586         fdNdEta[vertexPos]->SetBinError(bin, error);
 587
 588         dndeta /= ptCutOffCorrection;
 589         error  /= ptCutOffCorrection;
 590
 591         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
 592         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
 593
 594         //Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events) (comb. corr: %f)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents, combCorr);
 595       }
 596     }
 597   }
 598 }
 599
 600 //____________________________________________________________________
 601 Float_t dNdEtaAnalysis::GetVtxMin(Float_t eta)
 602 {
 603   // helper function for the SPD acceptance
 604   // the function returns the beginning of the acceptance window in z vertex position as function of eta
 605   // to get the maximum: -GetVtxMin(-eta)
 606
 607   Float_t a[2] = { -15, 0 };
 608   Float_t b[2] = { 0, -1.4 };
 609   Float_t c[2] = { 15, -2.2 };
 610
 611   Float_t meanAB[2];
 612   meanAB[0] = (b[0] + a[0]) / 2;
 613   meanAB[1] = (b[1] + a[1]) / 2;
 614
 615   Float_t meanBC[2];
 616   meanBC[0] = (c[0] + b[0]) / 2;
 617   meanBC[1] = (c[1] + b[1]) / 2;
 618
 619   Float_t mAB = (b[1] - a[1]) / (b[0] - a[0]);
 620   Float_t mBC = (c[1] - b[1]) / (c[0] - b[0]);
 621
 622   Float_t bAB = meanAB[1] - mAB * meanAB[0];
 623   Float_t bBC = meanBC[1] - mBC * meanBC[0];
 624
 625   if (eta > b[1])
 626     return 1.0 / mAB * eta - bAB / mAB;
 627
 628   return 1.0 / mBC * eta - bBC / mBC;
 629 }
 630
 631 //____________________________________________________________________
 632 void dNdEtaAnalysis::SaveHistograms()
 633 {
 634   // save the histograms to a directory with the name of this class (retrieved from TNamed)
 635
 636   gDirectory->mkdir(GetName());
 637   gDirectory->cd(GetName());
 638
 639   if (fData)
 640   {
 641     fData->SaveHistograms();
 642   }
 643   else
 644     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
 645
 646   if (fMult)
 647   {
 648     fMult->Write();
 649   }
 650   else
 651     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n");
 652
 653   if (fPtDist)
 654     fPtDist       ->Write();
 655   else
 656     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
 657
 658   for (Int_t i=0; i<kVertexBinning; ++i)
 659   {
 660     if (fdNdEta[i])
 661       fdNdEta[i]->Write();
 662     else
 663       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
 664
 665     if (fdNdEtaPtCutOffCorrected[i])
 666       fdNdEtaPtCutOffCorrected[i]->Write();
 667     else
 668       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
 669   }
 670
 671   TNamed named("fTag", fTag.Data());
 672   named.Write();
 673
 674   TParameter<Int_t> param("fAnalysisMode", fAnalysisMode);
 675   param.Write();
 676
 677   gDirectory->cd("../");
 678 }
 679
 680 void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
 681 {
 682   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
 683
 684   if (!dir)
 685     dir = GetName();
 686
 687   gDirectory->cd(dir);
 688
 689   fData->LoadHistograms();
 690   fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName()));
 691
 692   for (Int_t i=0; i<kVertexBinning; ++i)
 693   {
 694     fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
 695     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
 696   }
 697
 698   fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
 699
 700   if (static_cast<TNamed*> (gDirectory->Get("fTag")))
 701     fTag = (static_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle();
 702
 703   if (static_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))
 704     fAnalysisMode = (AliPWG0Helper::AnalysisMode) (static_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal();
 705
 706   gDirectory->cd("../");
 707 }
 708
 709 //____________________________________________________________________
 710 void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
 711 {
 712   // draws the histograms
 713
 714   if (!simple)
 715   {
 716     if (fData)
 717       fData->DrawHistograms(GetName());
 718
 719     TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
 720     canvas->Divide(2, 1);
 721
 722     canvas->cd(1);
 723     if (fdNdEtaPtCutOffCorrected[0])
 724       fdNdEtaPtCutOffCorrected[0]->DrawCopy();
 725
 726     if (fdNdEta[0])
 727     {
 728       fdNdEta[0]->SetLineColor(kRed);
 729       fdNdEta[0]->DrawCopy("SAME");
 730     }
 731
 732     canvas->cd(2);
 733     if (fPtDist)
 734       fPtDist->DrawCopy();
 735   }
 736
 737     // histograms for different vertices?
 738   if (kVertexBinning > 0)
 739   {
 740       // doesnt work, but i dont get it, giving up...
 741     TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450);
 742     TCanvas* canvas3 = 0;
 743     if (!simple)
 744       canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450);
 745
 746     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
 747     //printf("%d\n", yPads);
 748     //canvas2->Divide(2, yPads);
 749
 750     TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
 751
 752     for (Int_t i=0; i<kVertexBinning; ++i)
 753     {
 754       if (fdNdEtaPtCutOffCorrected[i])
 755       {
 756         canvas2->cd();
 757
 758         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
 759         fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
 760         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
 761       }
 762       if (canvas3 && fdNdEta[i])
 763       {
 764         canvas3->cd();
 765
 766         fdNdEta[i]->SetLineColor(i+1);
 767         fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
 768       }
 769     }
 770
 771     canvas2->cd();
 772     legend->Draw();
 773     canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
 774
 775     if (canvas3)
 776     {
 777       canvas3->cd();
 778       legend->Draw();
 779     }
 780   }
 781
 782   if (kVertexBinning == 3)
 783   {
 784      TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone"));
 785      TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2"));
 786
 787      if (clone && clone2)
 788      {
 789         TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
 790
 791         clone->Divide(fdNdEtaPtCutOffCorrected[0]);
 792         clone->GetYaxis()->SetRangeUser(0.95, 1.05);
 793         clone->DrawCopy();
 794
 795         clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
 796         clone2->DrawCopy("SAME");
 797
 798         TLine* line = new TLine(-1, 1, 1, 1);
 799         line->Draw();
 800
 801         canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName()));
 802      }
 803    }
 804 }
 805
 806 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
 807 {
 808   // Merges a list of dNdEtaAnalysis objects with this one.
 809   // This is needed for PROOF.
 810   // Returns the number of merged objects (including this)
 811
 812   if (!list)
 813     return 0;
 814
 815   if (list->IsEmpty())
 816     return 1;
 817
 818   TIterator* iter = list->MakeIterator();
 819   TObject* obj;
 820
 821   // sub collections
 822   const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning
 823   TList* collections[nCollections];
 824   for (Int_t i=0; i<nCollections; ++i)
 825     collections[i] = new TList;
 826
 827   Int_t count = 0;
 828   while ((obj = iter->Next()))
 829   {
 830     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
 831     if (entry == 0)
 832       continue;
 833
 834     collections[0]->Add(entry->fData);
 835     collections[1]->Add(entry->fMult);
 836     collections[2]->Add(entry->fPtDist);
 837
 838     for (Int_t i=0; i<kVertexBinning; ++i)
 839     {
 840       collections[3+i]->Add(entry->fdNdEta[i]);
 841       collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
 842     }
 843
 844     ++count;
 845   }
 846
 847   fData->Merge(collections[0]);
 848   fMult->Merge(collections[1]);
 849   fPtDist->Merge(collections[2]);
 850
 851   for (Int_t i=0; i<kVertexBinning; ++i)
 852   {
 853     fdNdEta[i]->Merge(collections[3+i]);
 854     fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]);
 855   }
 856
 857   for (Int_t i=0; i<nCollections; ++i)
 858     delete collections[i];
 859
 860   return count+1;
 861 }
 862
 863