adding flag for non-field data
[u/mrichter/AliRoot.git] / 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