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