adding flag for non-field data
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
CommitLineData
dc740de4 1/* $Id$ */
2
75ec0f41 3#include "dNdEtaAnalysis.h"
4
ceb5d1b5 5#include <TFile.h>
45e97e28 6#include <TH3F.h>
74fd10b3 7#include <TH2F.h>
8#include <TH1F.h>
ceb5d1b5 9#include <TMath.h>
10#include <TCanvas.h>
7029240a 11#include <TCollection.h>
12#include <TIterator.h>
13#include <TList.h>
fcf2fb36 14#include <TLegend.h>
74fd10b3 15#include <TLine.h>
745d6088 16#include <TParameter.h>
fcf2fb36 17
45e97e28 18#include "AlidNdEtaCorrection.h"
74fd10b3 19#include <AliCorrection.h>
20#include <AliPWG0Helper.h>
21#include <AliCorrectionMatrix2D.h>
22#include <AliCorrectionMatrix3D.h>
ceb5d1b5 23
75ec0f41 24//____________________________________________________________________
b7f4a1fd 25ClassImp(dNdEtaAnalysis)
75ec0f41 26
27//____________________________________________________________________
1afae8ff 28dNdEtaAnalysis::dNdEtaAnalysis() :
29 TNamed(),
30 fData(0),
7307d52c 31 fMult(0),
745d6088 32 fPtDist(0),
33 fAnalysisMode(AliPWG0Helper::kInvalid),
34 fTag()
1afae8ff 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//____________________________________________________________________
a6e0ebfe 46dNdEtaAnalysis::dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
16e24ca3 47 TNamed(name, title),
45e97e28 48 fData(0),
7307d52c 49 fMult(0),
745d6088 50 fPtDist(0),
51 fAnalysisMode(analysisMode),
52 fTag()
7029240a 53{
6bf0714d 54 // constructor
4dd2ad81 55
770a1f1d 56 fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode);
45e97e28 57
74fd10b3 58 // do not add this hists to the directory
59 Bool_t oldStatus = TH1::AddDirectoryStatus();
60 TH1::AddDirectory(kFALSE);
1afae8ff 61
7307d52c 62 fMult = new TH1F("TriggeredMultiplicity", "Triggered Events;raw multiplicity;entries", 1000, -0.5, 999.5);
63
567160d6 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());
74fd10b3 66
67 fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
1afae8ff 68
45e97e28 69 for (Int_t i=1; i<kVertexBinning; ++i)
7029240a 70 {
74fd10b3 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)));
7029240a 73 }
75ec0f41 74
567160d6 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());
74fd10b3 76
77 TH1::AddDirectory(oldStatus);
75ec0f41 78}
79
80//____________________________________________________________________
16e24ca3 81dNdEtaAnalysis::~dNdEtaAnalysis()
82{
83 // destructor
84
0ab29cfa 85 if (fData)
86 {
87 delete fData;
88 fData = 0;
89 }
16e24ca3 90
7307d52c 91 if (fMult)
92 {
93 delete fMult;
94 fMult = 0;
95 }
96
16e24ca3 97 for (Int_t i=0; i<kVertexBinning; ++i)
98 {
0ab29cfa 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 }
16e24ca3 109 }
1afae8ff 110
0ab29cfa 111 if (fPtDist)
112 {
113 delete fPtDist;
114 fPtDist = 0;
115 }
16e24ca3 116}
117
118//_____________________________________________________________________________
119dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
120 TNamed(c),
45e97e28 121 fData(0),
c61a7285 122 fMult(0),
745d6088 123 fPtDist(0),
124 fAnalysisMode(AliPWG0Helper::kInvalid),
125 fTag()
16e24ca3 126{
127 //
128 // dNdEtaAnalysis copy constructor
129 //
130
131 ((dNdEtaAnalysis &) c).Copy(*this);
132}
133
134//_____________________________________________________________________________
135dNdEtaAnalysis &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//_____________________________________________________________________________
146void dNdEtaAnalysis::Copy(TObject &c) const
147{
148 //
149 // Copy function
150 //
151
152 dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
153
74fd10b3 154 target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
7307d52c 155 target.fMult = dynamic_cast<TH1F*> (fMult->Clone());
16e24ca3 156
157 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 158 {
74fd10b3 159 target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
160 target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
1afae8ff 161 }
16e24ca3 162
74fd10b3 163 target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
45e97e28 164
745d6088 165 target.fAnalysisMode = fAnalysisMode;
166 target.fTag = fTag;
167
16e24ca3 168 TNamed::Copy((TNamed &) c);
169}
170
171//____________________________________________________________________
74fd10b3 172void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
6bf0714d 173{
174 // fills a track into the histograms
175
74fd10b3 176 fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
75ec0f41 177}
178
179//____________________________________________________________________
74fd10b3 180void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
6bf0714d 181{
182 // fills an event into the histograms
183
74fd10b3 184 fData->GetEventCorrection()->FillMeas(vtx, n);
75ec0f41 185}
186
187//____________________________________________________________________
7307d52c 188void dNdEtaAnalysis::FillTriggeredEvent(Float_t n)
189{
190 // fills a triggered event into the histograms
191
192 fMult->Fill(n);
193}
194
195//____________________________________________________________________
0fc41645 196void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag)
fcf2fb36 197{
74fd10b3 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
3dfa46a4 201 // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
74fd10b3 202 //
203
51f6de65 204 fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut);
745d6088 205 Printf("\n\n%s", fTag.Data());
3dfa46a4 206
74fd10b3 207 // set corrections to 1
208 fData->SetCorrectionToUnity();
209
210 if (correction && correctionType != AlidNdEtaCorrection::kNone)
211 {
5a6310fe 212 TH3* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
213 TH2* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
74fd10b3 214
215 if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
216 trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
6bf0714d 217
dd367a14 218 if (correctionType >= AlidNdEtaCorrection::kVertexReco)
74fd10b3 219 {
220 trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
221 eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
7307d52c 222
0fc41645 223 // set bin with multiplicity 0 to unity (correction has no meaning in this bin)
7307d52c 224 for (Int_t i=0; i<=eventCorr->GetNbinsX()+1; i++)
225 eventCorr->SetBinContent(i, 1, 1);
dd367a14 226 }
74fd10b3 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
fcf2fb36 252 printf("INFO: No correction applied\n");
253
567160d6 254 TH2F* rawMeasured = (TH2F*) fData->GetEventCorrection()->GetMeasuredHistogram()->Clone("rawMeasured");
255
74fd10b3 256 fData->Multiply();
7307d52c 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
a7f69e56 265 // at the same time calculate expectation from MC (not used, just to check the value)
7307d52c 266
567160d6 267 //TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
7307d52c 268 TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
269
567160d6 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
7307d52c 274 //new TCanvas; correctedEvents->DrawCopy("TEXT");
275
276 // start above 0 mult. bin with integration
567160d6 277 TH1* vertexDist = rawMeasured->ProjectionX("vertexdist_measured", 2, rawMeasured->GetNbinsY()+1);
7307d52c 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
51f6de65 285 TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
286 kineBias->Reset();
287
a7f69e56 288 // loop over vertex bins
567160d6 289 for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
7307d52c 290 {
745d6088 291 Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
7307d52c 292 Double_t events = alpha * triggeredEventsWith0Mult;
293
567160d6 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);
51f6de65 299 kineBias->SetBinContent(i, fZ);
567160d6 300
301 events *= fZ;
302
7307d52c 303 // multiply with trigger correction if set above
304 events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
305
a7f69e56 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);
7307d52c 311
312 correctedEvents->SetBinContent(i, 1, events);
313 }
a7f69e56 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));
7307d52c 316
317 //new TCanvas; correctedEvents->DrawCopy("TEXT");
69b09e3b 318 //new TCanvas; kineBias->DrawCopy();
7307d52c 319 }
320
3dfa46a4 321 fData->PrintInfo(ptCut);
74fd10b3 322
5a6310fe 323 TH3* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
74fd10b3 324
325 // integrate multiplicity axis out (include under/overflow bins!!!)
5a6310fe 326 TH2* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
3dfa46a4 327
7307d52c 328 TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
fcf2fb36 329
1afae8ff 330 // create pt hist
a7f69e56 331 if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
1afae8ff 332 {
333 // reset all ranges
74fd10b3 334 dataHist->GetXaxis()->SetRange(0, 0);
335 dataHist->GetYaxis()->SetRange(0, 0);
336 dataHist->GetZaxis()->SetRange(0, 0);
1afae8ff 337
338 // vtx cut
74fd10b3 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);
1afae8ff 343
74fd10b3 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;
1afae8ff 349
74fd10b3 350 TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
1afae8ff 351
74fd10b3 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 }
1afae8ff 358
74fd10b3 359 delete ptHist;
360 }
361 else
362 printf("ERROR: nEvents is 0!\n");
1afae8ff 363 }
364
365 // reset all ranges
74fd10b3 366 dataHist->GetXaxis()->SetRange(0, 0);
367 dataHist->GetYaxis()->SetRange(0, 0);
368 dataHist->GetZaxis()->SetRange(0, 0);
1afae8ff 369
51f6de65 370 // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD)
1afae8ff 371 Int_t ptLowBin = 1;
a7f69e56 372 if (ptCut > 0 && (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS))
74fd10b3 373 ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
69b09e3b 374
375 //new TCanvas; dataHist->DrawCopy();
1afae8ff 376
69b09e3b 377 //dataHist->Sumw2();
7307d52c 378 dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
51f6de65 379 printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
74fd10b3 380 TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
69b09e3b 381
382 //new TCanvas; vtxVsEta->Draw("COLZ");
74fd10b3 383
384 dataHist->GetZaxis()->SetRange(0, 0);
385 vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
386 vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
1afae8ff 387
847489f7 388 if (vtxVsEta == 0)
389 {
51f6de65 390 printf("ERROR: pt/multiplicity integration failed\n");
847489f7 391 return;
fcf2fb36 392 }
393
51f6de65 394 //new TCanvas(tag, tag, 800, 600); vtxVsEta->DrawCopy("COLZ");
395
567160d6 396 // clear result histograms
397 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
398 {
399 fdNdEta[vertexPos]->Reset();
400 fdNdEtaPtCutOffCorrected[vertexPos]->Reset();
401 }
144ff489 402
567160d6 403 const Float_t vertexRangeBegin[kVertexBinning] = { -9.99, -9.99, 0.01 };
404 const Float_t vertexRangeEnd[kVertexBinning] = { 9.99, -0.01, 9.99 };
b4b9cacc 405
74fd10b3 406 for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
5af55649 407 {
745d6088 408 // loop over vertex ranges
7029240a 409 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
410 {
745d6088 411 Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]);
412 Int_t vertexBinEnd = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
fcf2fb36 413
ea441adf 414 const Int_t *binBegin = 0;
69b09e3b 415 const Int_t maxBins = 60;
1c15d51a 416
417 // adjust acceptance range
69b09e3b 418 // produce with drawPlots.C: DetermineAcceptance(...)
a7f69e56 419 if (fAnalysisMode & AliPWG0Helper::kSPD)
fcf2fb36 420 {
1c15d51a 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
69b09e3b 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
1c15d51a 423
1cbdb1a6 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
567160d6 425
1c15d51a 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
69b09e3b 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
1c15d51a 428
429 binBegin = binBeginSPD;
fcf2fb36 430 }
a7f69e56 431 else if (fAnalysisMode & AliPWG0Helper::kTPC)
ea441adf 432 {
69b09e3b 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
c17301f3 435
ea441adf 436 binBegin = binBeginTPC;
437 }
a7f69e56 438 else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
ea441adf 439 {
440 // TODO create map
441 }
442
51f6de65 443 Int_t vtxBegin = 1;
69b09e3b 444 Int_t vtxEnd = maxBins;
51f6de65 445
446 if (binBegin)
447 {
448 vtxBegin = binBegin[iEta - 1];
69b09e3b 449 vtxEnd = 18 + 1 - binBegin[maxBins - iEta];
51f6de65 450 }
451 else
452 Printf("WARNING: No acceptance applied!");
ea441adf 453
454 // eta range not accessible
455 if (vtxBegin == -1)
456 continue;
457
ea441adf 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));
1c15d51a 461
b9dbe907 462 //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
51f6de65 463 //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
ea441adf 464 vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
465 vertexBinEnd = TMath::Min(vertexBinEnd, vtxEnd);
51f6de65 466 //Printf(" after: %d %d", vertexBinBegin, vertexBinEnd);
7029240a 467
745d6088 468 // no data for this bin
469 if (vertexBinBegin > vertexBinEnd)
567160d6 470 {
471 //Printf("Bin empty. Skipped");
745d6088 472 continue;
567160d6 473 }
b4b9cacc 474
1c15d51a 475 Float_t totalEvents = 0;
5af55649 476 Float_t sum = 0;
477 Float_t sumError2 = 0;
1c15d51a 478 Float_t unusedTracks = 0;
479 Float_t unusedEvents = 0;
480 for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++)
567160d6 481 {
1c15d51a 482 if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd)
5af55649 483 {
1c15d51a 484 if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
485 {
486 sum += vtxVsEta->GetBinContent(iVtx, iEta);
2e88424e 487
1c15d51a 488 if (sumError2 > 10e30)
489 Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
2e88424e 490
1c15d51a 491 sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
492 }
493 totalEvents += vertexHist->GetBinContent(iVtx);
5af55649 494 }
1c15d51a 495 else
496 {
497 unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta);
69b09e3b 498 unusedEvents += vertexHist->GetBinContent(iVtx);
1c15d51a 499 }
500 }
501
502 if (totalEvents == 0)
503 {
504 printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
505 continue;
5af55649 506 }
7029240a 507
1afae8ff 508 Float_t ptCutOffCorrection = 1;
1afae8ff 509
51f6de65 510 // find pt cut off correction factor
a7f69e56 511 if ((fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) && (fAnalysisMode & AliPWG0Helper::kFieldOn))
847489f7 512 {
51f6de65 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 }
847489f7 521 }
522
69b09e3b 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);
9e952c39 524
144ff489 525 Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
745d6088 526 if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX())
144ff489 527 {
7307d52c 528 Float_t dndeta = sum / totalEvents;
529 Float_t error = TMath::Sqrt(sumError2) / totalEvents;
7029240a 530
144ff489 531 dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
532 error = error / fdNdEta[vertexPos]->GetBinWidth(bin);
7029240a 533
144ff489 534 fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
535 fdNdEta[vertexPos]->SetBinError(bin, error);
1afae8ff 536
144ff489 537 dndeta /= ptCutOffCorrection;
538 error /= ptCutOffCorrection;
1afae8ff 539
144ff489 540 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
541 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
2e88424e 542
51f6de65 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);
144ff489 544 }
75ec0f41 545 }
847489f7 546 }
75ec0f41 547}
548
75ec0f41 549//____________________________________________________________________
567160d6 550Float_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//____________________________________________________________________
6bf0714d 581void dNdEtaAnalysis::SaveHistograms()
582{
583 // save the histograms to a directory with the name of this class (retrieved from TNamed)
75ec0f41 584
7029240a 585 gDirectory->mkdir(GetName());
586 gDirectory->cd(GetName());
5fbd0b17 587
1afae8ff 588 if (fData)
589 {
74fd10b3 590 fData->SaveHistograms();
1afae8ff 591 }
592 else
593 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
594
7307d52c 595 if (fMult)
596 {
597 fMult->Write();
598 }
599 else
600 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n");
601
d09fb536 602 if (fPtDist)
603 fPtDist ->Write();
604 else
605 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
606
7029240a 607 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 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 }
75ec0f41 619
745d6088 620 TNamed named("fTag", fTag.Data());
621 named.Write();
622
623 TParameter<Int_t> param("fAnalysisMode", fAnalysisMode);
624 param.Write();
625
75ec0f41 626 gDirectory->cd("../");
627}
628
74fd10b3 629void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
5fbd0b17 630{
6bf0714d 631 // loads the histograms from a directory with the name of this class (retrieved from TNamed)
632
74fd10b3 633 if (!dir)
634 dir = GetName();
5fbd0b17 635
74fd10b3 636 gDirectory->cd(dir);
5fbd0b17 637
74fd10b3 638 fData->LoadHistograms();
7307d52c 639 fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName()));
5fbd0b17 640
641 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 642 {
74fd10b3 643 fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
644 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
1afae8ff 645 }
5fbd0b17 646
74fd10b3 647 fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
0ab29cfa 648
1c15d51a 649 if (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))
650 fTag = (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle();
745d6088 651
1c15d51a 652 if (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))
653 fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal();
745d6088 654
5fbd0b17 655 gDirectory->cd("../");
656}
657
ceb5d1b5 658//____________________________________________________________________
74fd10b3 659void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
ceb5d1b5 660{
6bf0714d 661 // draws the histograms
7307d52c 662
74fd10b3 663 if (!simple)
664 {
665 if (fData)
666 fData->DrawHistograms(GetName());
6bf0714d 667
74fd10b3 668 TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
669 canvas->Divide(2, 1);
ceb5d1b5 670
74fd10b3 671 canvas->cd(1);
672 if (fdNdEtaPtCutOffCorrected[0])
144ff489 673 fdNdEtaPtCutOffCorrected[0]->DrawCopy();
ceb5d1b5 674
74fd10b3 675 if (fdNdEta[0])
676 {
677 fdNdEta[0]->SetLineColor(kRed);
144ff489 678 fdNdEta[0]->DrawCopy("SAME");
74fd10b3 679 }
1afae8ff 680
74fd10b3 681 canvas->cd(2);
682 if (fPtDist)
144ff489 683 fPtDist->DrawCopy();
1afae8ff 684 }
685
fcf2fb36 686 // histograms for different vertices?
687 if (kVertexBinning > 0)
688 {
689 // doesnt work, but i dont get it, giving up...
74fd10b3 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
fcf2fb36 695 //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
696 //printf("%d\n", yPads);
697 //canvas2->Divide(2, yPads);
698
0448e811 699 TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
fcf2fb36 700
5af55649 701 for (Int_t i=0; i<kVertexBinning; ++i)
fcf2fb36 702 {
1afae8ff 703 if (fdNdEtaPtCutOffCorrected[i])
fcf2fb36 704 {
74fd10b3 705 canvas2->cd();
706
1afae8ff 707 fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
144ff489 708 fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
1afae8ff 709 legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
fcf2fb36 710 }
74fd10b3 711 if (canvas3 && fdNdEta[i])
712 {
713 canvas3->cd();
714
715 fdNdEta[i]->SetLineColor(i+1);
144ff489 716 fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
74fd10b3 717 }
fcf2fb36 718 }
719
74fd10b3 720 canvas2->cd();
fcf2fb36 721 legend->Draw();
74fd10b3 722 canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
723
724 if (canvas3)
725 {
726 canvas3->cd();
727 legend->Draw();
728 }
fcf2fb36 729 }
7307d52c 730
74fd10b3 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"));
7307d52c 735
74fd10b3 736 if (clone && clone2)
737 {
738 TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
7307d52c 739
74fd10b3 740 clone->Divide(fdNdEtaPtCutOffCorrected[0]);
741 clone->GetYaxis()->SetRangeUser(0.95, 1.05);
144ff489 742 clone->DrawCopy();
7307d52c 743
74fd10b3 744 clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
144ff489 745 clone2->DrawCopy("SAME");
74fd10b3 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 }
7029240a 753}
754
755Long64_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
7307d52c 771 const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning
7029240a 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
45e97e28 783 collections[0]->Add(entry->fData);
0fc41645 784 collections[1]->Add(entry->fMult);
785 collections[2]->Add(entry->fPtDist);
7029240a 786
787 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 788 {
7307d52c 789 collections[3+i]->Add(entry->fdNdEta[i]);
790 collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
1afae8ff 791 }
7029240a 792
793 ++count;
794 }
795
45e97e28 796 fData->Merge(collections[0]);
7307d52c 797 fMult->Merge(collections[1]);
798 fPtDist->Merge(collections[2]);
799
7029240a 800 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 801 {
7307d52c 802 fdNdEta[i]->Merge(collections[3+i]);
803 fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]);
1afae8ff 804 }
7029240a 805
806 for (Int_t i=0; i<nCollections; ++i)
807 delete collections[i];
808
809 return count+1;
ceb5d1b5 810}
567160d6 811
812