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