]>
Commit | Line | Data |
---|---|---|
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 | 25 | ClassImp(dNdEtaAnalysis) |
75ec0f41 | 26 | |
1afae8ff | 27 | //____________________________________________________________________ |
28 | dNdEtaAnalysis::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 | ||
75ec0f41 | 45 | //____________________________________________________________________ |
770a1f1d | 46 | dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, 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 | ||
16e24ca3 | 80 | //____________________________________________________________________ |
81 | dNdEtaAnalysis::~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 | //_____________________________________________________________________________ | |
119 | dNdEtaAnalysis::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 | //_____________________________________________________________________________ | |
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 | ||
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 | ||
75ec0f41 | 171 | //____________________________________________________________________ |
74fd10b3 | 172 | void 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 | 180 | void 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 | ||
7307d52c | 187 | //____________________________________________________________________ |
188 | void dNdEtaAnalysis::FillTriggeredEvent(Float_t n) | |
189 | { | |
190 | // fills a triggered event into the histograms | |
191 | ||
192 | fMult->Fill(n); | |
193 | } | |
194 | ||
75ec0f41 | 195 | //____________________________________________________________________ |
0fc41645 | 196 | void 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 | |
265 | ||
567160d6 | 266 | //TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram(); |
7307d52c | 267 | TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram(); |
268 | ||
567160d6 | 269 | TH2* eTrig = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram(); |
270 | TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram(); | |
51f6de65 | 271 | //TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, rawMeasured->GetNbinsY()+1); |
567160d6 | 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 | ||
567160d6 | 288 | for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++) |
7307d52c | 289 | { |
745d6088 | 290 | Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex; |
7307d52c | 291 | Double_t events = alpha * triggeredEventsWith0Mult; |
292 | ||
567160d6 | 293 | if (eTrigVtx_projx->GetBinContent(i) == 0) |
294 | continue; | |
295 | ||
296 | Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) * | |
297 | eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1); | |
51f6de65 | 298 | kineBias->SetBinContent(i, fZ); |
567160d6 | 299 | |
300 | events *= fZ; | |
301 | ||
7307d52c | 302 | // multiply with trigger correction if set above |
303 | events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1); | |
304 | ||
567160d6 | 305 | Printf("Bin %d, alpha is %.2f, fZ is %.3f, number of events with 0 mult.: %.2f", i, alpha * 100., fZ, events); |
7307d52c | 306 | |
307 | correctedEvents->SetBinContent(i, 1, events); | |
308 | } | |
309 | ||
310 | //new TCanvas; correctedEvents->DrawCopy("TEXT"); | |
51f6de65 | 311 | //new TCanvas; kineBias->Draw(); |
7307d52c | 312 | } |
313 | ||
3dfa46a4 | 314 | fData->PrintInfo(ptCut); |
74fd10b3 | 315 | |
5a6310fe | 316 | TH3* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram(); |
74fd10b3 | 317 | |
318 | // integrate multiplicity axis out (include under/overflow bins!!!) | |
5a6310fe | 319 | TH2* tmp = fData->GetEventCorrection()->GetGeneratedHistogram(); |
3dfa46a4 | 320 | |
7307d52c | 321 | TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e"); |
fcf2fb36 | 322 | |
1afae8ff | 323 | // create pt hist |
51f6de65 | 324 | if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) |
1afae8ff | 325 | { |
326 | // reset all ranges | |
74fd10b3 | 327 | dataHist->GetXaxis()->SetRange(0, 0); |
328 | dataHist->GetYaxis()->SetRange(0, 0); | |
329 | dataHist->GetZaxis()->SetRange(0, 0); | |
1afae8ff | 330 | |
331 | // vtx cut | |
74fd10b3 | 332 | Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5); |
333 | Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5); | |
334 | dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd); | |
335 | Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd); | |
1afae8ff | 336 | |
74fd10b3 | 337 | if (nEvents > 0) |
338 | { | |
339 | // eta cut | |
340 | dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8)); | |
341 | Float_t etaWidth = 1.6; | |
1afae8ff | 342 | |
74fd10b3 | 343 | TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze")); |
1afae8ff | 344 | |
74fd10b3 | 345 | for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i) |
346 | { | |
347 | Float_t binSize = fPtDist->GetBinWidth(i); | |
348 | fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth); | |
349 | fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth); | |
350 | } | |
1afae8ff | 351 | |
74fd10b3 | 352 | delete ptHist; |
353 | } | |
354 | else | |
355 | printf("ERROR: nEvents is 0!\n"); | |
1afae8ff | 356 | } |
357 | ||
358 | // reset all ranges | |
74fd10b3 | 359 | dataHist->GetXaxis()->SetRange(0, 0); |
360 | dataHist->GetYaxis()->SetRange(0, 0); | |
361 | dataHist->GetZaxis()->SetRange(0, 0); | |
1afae8ff | 362 | |
51f6de65 | 363 | // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD) |
1afae8ff | 364 | Int_t ptLowBin = 1; |
51f6de65 | 365 | if (ptCut > 0 && fAnalysisMode != AliPWG0Helper::kSPD) |
74fd10b3 | 366 | ptLowBin = dataHist->GetZaxis()->FindBin(ptCut); |
1afae8ff | 367 | |
7307d52c | 368 | dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1); |
51f6de65 | 369 | printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1); |
74fd10b3 | 370 | TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e")); |
371 | ||
372 | dataHist->GetZaxis()->SetRange(0, 0); | |
373 | vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle()); | |
374 | vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle()); | |
1afae8ff | 375 | |
847489f7 | 376 | if (vtxVsEta == 0) |
377 | { | |
51f6de65 | 378 | printf("ERROR: pt/multiplicity integration failed\n"); |
847489f7 | 379 | return; |
fcf2fb36 | 380 | } |
381 | ||
51f6de65 | 382 | //new TCanvas(tag, tag, 800, 600); vtxVsEta->DrawCopy("COLZ"); |
383 | ||
567160d6 | 384 | // clear result histograms |
385 | for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos) | |
386 | { | |
387 | fdNdEta[vertexPos]->Reset(); | |
388 | fdNdEtaPtCutOffCorrected[vertexPos]->Reset(); | |
389 | } | |
144ff489 | 390 | |
567160d6 | 391 | const Float_t vertexRangeBegin[kVertexBinning] = { -9.99, -9.99, 0.01 }; |
392 | const Float_t vertexRangeEnd[kVertexBinning] = { 9.99, -0.01, 9.99 }; | |
b4b9cacc | 393 | |
74fd10b3 | 394 | for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++) |
5af55649 | 395 | { |
745d6088 | 396 | // loop over vertex ranges |
7029240a | 397 | for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos) |
398 | { | |
745d6088 | 399 | Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]); |
400 | Int_t vertexBinEnd = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]); | |
fcf2fb36 | 401 | |
ea441adf | 402 | const Int_t *binBegin = 0; |
1c15d51a | 403 | |
404 | // adjust acceptance range | |
745d6088 | 405 | if (fAnalysisMode == AliPWG0Helper::kSPD) |
fcf2fb36 | 406 | { |
1c15d51a | 407 | //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 |
ea441adf | 408 | 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 | 409 | |
1cbdb1a6 | 410 | //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 | 411 | |
1c15d51a | 412 | //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 |
413 | ||
414 | binBegin = binBeginSPD; | |
fcf2fb36 | 415 | } |
ea441adf | 416 | else if (fAnalysisMode == AliPWG0Helper::kTPC) |
417 | { | |
c17301f3 | 418 | 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 |
419 | ||
ea441adf | 420 | binBegin = binBeginTPC; |
421 | } | |
422 | else if (fAnalysisMode == AliPWG0Helper::kTPCITS) | |
423 | { | |
424 | // TODO create map | |
425 | } | |
426 | ||
51f6de65 | 427 | Int_t vtxBegin = 1; |
428 | Int_t vtxEnd = 30; | |
429 | ||
430 | if (binBegin) | |
431 | { | |
432 | vtxBegin = binBegin[iEta - 1]; | |
433 | vtxEnd = 18 + 1 - binBegin[30 - iEta]; | |
434 | } | |
435 | else | |
436 | Printf("WARNING: No acceptance applied!"); | |
ea441adf | 437 | |
438 | // eta range not accessible | |
439 | if (vtxBegin == -1) | |
440 | continue; | |
441 | ||
ea441adf | 442 | //Printf("%d %d | %d %d", vtxBegin, vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)), vtxEnd, vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta))); |
443 | //vtxBegin = vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)); | |
444 | //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta)); | |
1c15d51a | 445 | |
446 | Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta); | |
51f6de65 | 447 | //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd); |
ea441adf | 448 | vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin); |
449 | vertexBinEnd = TMath::Min(vertexBinEnd, vtxEnd); | |
51f6de65 | 450 | //Printf(" after: %d %d", vertexBinBegin, vertexBinEnd); |
7029240a | 451 | |
745d6088 | 452 | // no data for this bin |
453 | if (vertexBinBegin > vertexBinEnd) | |
567160d6 | 454 | { |
455 | //Printf("Bin empty. Skipped"); | |
745d6088 | 456 | continue; |
567160d6 | 457 | } |
b4b9cacc | 458 | |
1c15d51a | 459 | Float_t totalEvents = 0; |
5af55649 | 460 | Float_t sum = 0; |
461 | Float_t sumError2 = 0; | |
1c15d51a | 462 | Float_t unusedTracks = 0; |
463 | Float_t unusedEvents = 0; | |
464 | for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++) | |
567160d6 | 465 | { |
1c15d51a | 466 | if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd) |
5af55649 | 467 | { |
1c15d51a | 468 | if (vtxVsEta->GetBinContent(iVtx, iEta) != 0) |
469 | { | |
470 | sum += vtxVsEta->GetBinContent(iVtx, iEta); | |
2e88424e | 471 | |
1c15d51a | 472 | if (sumError2 > 10e30) |
473 | Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... "); | |
2e88424e | 474 | |
1c15d51a | 475 | sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2); |
476 | } | |
477 | totalEvents += vertexHist->GetBinContent(iVtx); | |
5af55649 | 478 | } |
1c15d51a | 479 | else |
480 | { | |
481 | unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta); | |
482 | unusedEvents +=vertexHist->GetBinContent(iVtx); | |
483 | } | |
484 | } | |
485 | ||
486 | if (totalEvents == 0) | |
487 | { | |
488 | printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd); | |
489 | continue; | |
5af55649 | 490 | } |
7029240a | 491 | |
1afae8ff | 492 | Float_t ptCutOffCorrection = 1; |
1afae8ff | 493 | |
51f6de65 | 494 | // find pt cut off correction factor |
495 | if (fAnalysisMode != AliPWG0Helper::kSPD) | |
847489f7 | 496 | { |
51f6de65 | 497 | if (correction && ptCut > 0) |
498 | ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd); | |
499 | ||
500 | if (ptCutOffCorrection <= 0) | |
501 | { | |
502 | printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd); | |
503 | continue; | |
504 | } | |
847489f7 | 505 | } |
506 | ||
9e952c39 | 507 | //printf("Eta: %d Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f\n", iEta, vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection); |
508 | ||
144ff489 | 509 | Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta)); |
745d6088 | 510 | if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX()) |
144ff489 | 511 | { |
7307d52c | 512 | Float_t dndeta = sum / totalEvents; |
513 | Float_t error = TMath::Sqrt(sumError2) / totalEvents; | |
7029240a | 514 | |
144ff489 | 515 | dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin); |
516 | error = error / fdNdEta[vertexPos]->GetBinWidth(bin); | |
7029240a | 517 | |
144ff489 | 518 | fdNdEta[vertexPos]->SetBinContent(bin, dndeta); |
519 | fdNdEta[vertexPos]->SetBinError(bin, error); | |
1afae8ff | 520 | |
144ff489 | 521 | dndeta /= ptCutOffCorrection; |
522 | error /= ptCutOffCorrection; | |
1afae8ff | 523 | |
144ff489 | 524 | fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta); |
525 | fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error); | |
2e88424e | 526 | |
51f6de65 | 527 | //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 | 528 | } |
75ec0f41 | 529 | } |
847489f7 | 530 | } |
75ec0f41 | 531 | } |
532 | ||
567160d6 | 533 | //____________________________________________________________________ |
534 | Float_t dNdEtaAnalysis::GetVtxMin(Float_t eta) | |
535 | { | |
536 | // helper function for the SPD acceptance | |
537 | // the function returns the beginning of the acceptance window in z vertex position as function of eta | |
538 | // to get the maximum: -GetVtxMin(-eta) | |
539 | ||
540 | Float_t a[2] = { -15, 0 }; | |
541 | Float_t b[2] = { 0, -1.4 }; | |
542 | Float_t c[2] = { 15, -2.2 }; | |
543 | ||
544 | Float_t meanAB[2]; | |
545 | meanAB[0] = (b[0] + a[0]) / 2; | |
546 | meanAB[1] = (b[1] + a[1]) / 2; | |
547 | ||
548 | Float_t meanBC[2]; | |
549 | meanBC[0] = (c[0] + b[0]) / 2; | |
550 | meanBC[1] = (c[1] + b[1]) / 2; | |
551 | ||
552 | Float_t mAB = (b[1] - a[1]) / (b[0] - a[0]); | |
553 | Float_t mBC = (c[1] - b[1]) / (c[0] - b[0]); | |
554 | ||
555 | Float_t bAB = meanAB[1] - mAB * meanAB[0]; | |
556 | Float_t bBC = meanBC[1] - mBC * meanBC[0]; | |
557 | ||
558 | if (eta > b[1]) | |
559 | return 1.0 / mAB * eta - bAB / mAB; | |
560 | ||
561 | return 1.0 / mBC * eta - bBC / mBC; | |
562 | } | |
563 | ||
75ec0f41 | 564 | //____________________________________________________________________ |
6bf0714d | 565 | void dNdEtaAnalysis::SaveHistograms() |
566 | { | |
567 | // save the histograms to a directory with the name of this class (retrieved from TNamed) | |
75ec0f41 | 568 | |
7029240a | 569 | gDirectory->mkdir(GetName()); |
570 | gDirectory->cd(GetName()); | |
5fbd0b17 | 571 | |
1afae8ff | 572 | if (fData) |
573 | { | |
74fd10b3 | 574 | fData->SaveHistograms(); |
1afae8ff | 575 | } |
576 | else | |
577 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n"); | |
578 | ||
7307d52c | 579 | if (fMult) |
580 | { | |
581 | fMult->Write(); | |
582 | } | |
583 | else | |
584 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n"); | |
585 | ||
d09fb536 | 586 | if (fPtDist) |
587 | fPtDist ->Write(); | |
588 | else | |
589 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n"); | |
590 | ||
7029240a | 591 | for (Int_t i=0; i<kVertexBinning; ++i) |
1afae8ff | 592 | { |
593 | if (fdNdEta[i]) | |
594 | fdNdEta[i]->Write(); | |
595 | else | |
596 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i); | |
597 | ||
598 | if (fdNdEtaPtCutOffCorrected[i]) | |
599 | fdNdEtaPtCutOffCorrected[i]->Write(); | |
600 | else | |
601 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i); | |
602 | } | |
75ec0f41 | 603 | |
745d6088 | 604 | TNamed named("fTag", fTag.Data()); |
605 | named.Write(); | |
606 | ||
607 | TParameter<Int_t> param("fAnalysisMode", fAnalysisMode); | |
608 | param.Write(); | |
609 | ||
75ec0f41 | 610 | gDirectory->cd("../"); |
611 | } | |
612 | ||
74fd10b3 | 613 | void dNdEtaAnalysis::LoadHistograms(const Char_t* dir) |
5fbd0b17 | 614 | { |
6bf0714d | 615 | // loads the histograms from a directory with the name of this class (retrieved from TNamed) |
616 | ||
74fd10b3 | 617 | if (!dir) |
618 | dir = GetName(); | |
5fbd0b17 | 619 | |
74fd10b3 | 620 | gDirectory->cd(dir); |
5fbd0b17 | 621 | |
74fd10b3 | 622 | fData->LoadHistograms(); |
7307d52c | 623 | fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName())); |
5fbd0b17 | 624 | |
625 | for (Int_t i=0; i<kVertexBinning; ++i) | |
1afae8ff | 626 | { |
74fd10b3 | 627 | fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName())); |
628 | fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName())); | |
1afae8ff | 629 | } |
5fbd0b17 | 630 | |
74fd10b3 | 631 | fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName())); |
0ab29cfa | 632 | |
1c15d51a | 633 | if (dynamic_cast<TNamed*> (gDirectory->Get("fTag"))) |
634 | fTag = (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle(); | |
745d6088 | 635 | |
1c15d51a | 636 | if (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode"))) |
637 | fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal(); | |
745d6088 | 638 | |
5fbd0b17 | 639 | gDirectory->cd("../"); |
640 | } | |
641 | ||
ceb5d1b5 | 642 | //____________________________________________________________________ |
74fd10b3 | 643 | void dNdEtaAnalysis::DrawHistograms(Bool_t simple) |
ceb5d1b5 | 644 | { |
6bf0714d | 645 | // draws the histograms |
7307d52c | 646 | |
74fd10b3 | 647 | if (!simple) |
648 | { | |
649 | if (fData) | |
650 | fData->DrawHistograms(GetName()); | |
6bf0714d | 651 | |
74fd10b3 | 652 | TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400); |
653 | canvas->Divide(2, 1); | |
ceb5d1b5 | 654 | |
74fd10b3 | 655 | canvas->cd(1); |
656 | if (fdNdEtaPtCutOffCorrected[0]) | |
144ff489 | 657 | fdNdEtaPtCutOffCorrected[0]->DrawCopy(); |
ceb5d1b5 | 658 | |
74fd10b3 | 659 | if (fdNdEta[0]) |
660 | { | |
661 | fdNdEta[0]->SetLineColor(kRed); | |
144ff489 | 662 | fdNdEta[0]->DrawCopy("SAME"); |
74fd10b3 | 663 | } |
1afae8ff | 664 | |
74fd10b3 | 665 | canvas->cd(2); |
666 | if (fPtDist) | |
144ff489 | 667 | fPtDist->DrawCopy(); |
1afae8ff | 668 | } |
669 | ||
fcf2fb36 | 670 | // histograms for different vertices? |
671 | if (kVertexBinning > 0) | |
672 | { | |
673 | // doesnt work, but i dont get it, giving up... | |
74fd10b3 | 674 | TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450); |
675 | TCanvas* canvas3 = 0; | |
676 | if (!simple) | |
677 | canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450); | |
678 | ||
fcf2fb36 | 679 | //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2); |
680 | //printf("%d\n", yPads); | |
681 | //canvas2->Divide(2, yPads); | |
682 | ||
0448e811 | 683 | TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4); |
fcf2fb36 | 684 | |
5af55649 | 685 | for (Int_t i=0; i<kVertexBinning; ++i) |
fcf2fb36 | 686 | { |
1afae8ff | 687 | if (fdNdEtaPtCutOffCorrected[i]) |
fcf2fb36 | 688 | { |
74fd10b3 | 689 | canvas2->cd(); |
690 | ||
1afae8ff | 691 | fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1); |
144ff489 | 692 | fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME"); |
1afae8ff | 693 | legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1)); |
fcf2fb36 | 694 | } |
74fd10b3 | 695 | if (canvas3 && fdNdEta[i]) |
696 | { | |
697 | canvas3->cd(); | |
698 | ||
699 | fdNdEta[i]->SetLineColor(i+1); | |
144ff489 | 700 | fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME"); |
74fd10b3 | 701 | } |
fcf2fb36 | 702 | } |
703 | ||
74fd10b3 | 704 | canvas2->cd(); |
fcf2fb36 | 705 | legend->Draw(); |
74fd10b3 | 706 | canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName())); |
707 | ||
708 | if (canvas3) | |
709 | { | |
710 | canvas3->cd(); | |
711 | legend->Draw(); | |
712 | } | |
fcf2fb36 | 713 | } |
7307d52c | 714 | |
74fd10b3 | 715 | if (kVertexBinning == 3) |
716 | { | |
717 | TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone")); | |
718 | TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2")); | |
7307d52c | 719 | |
74fd10b3 | 720 | if (clone && clone2) |
721 | { | |
722 | TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450); | |
7307d52c | 723 | |
74fd10b3 | 724 | clone->Divide(fdNdEtaPtCutOffCorrected[0]); |
725 | clone->GetYaxis()->SetRangeUser(0.95, 1.05); | |
144ff489 | 726 | clone->DrawCopy(); |
7307d52c | 727 | |
74fd10b3 | 728 | clone2->Divide(fdNdEtaPtCutOffCorrected[0]); |
144ff489 | 729 | clone2->DrawCopy("SAME"); |
74fd10b3 | 730 | |
731 | TLine* line = new TLine(-1, 1, 1, 1); | |
732 | line->Draw(); | |
733 | ||
734 | canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName())); | |
735 | } | |
736 | } | |
7029240a | 737 | } |
738 | ||
739 | Long64_t dNdEtaAnalysis::Merge(TCollection* list) | |
740 | { | |
741 | // Merges a list of dNdEtaAnalysis objects with this one. | |
742 | // This is needed for PROOF. | |
743 | // Returns the number of merged objects (including this) | |
744 | ||
745 | if (!list) | |
746 | return 0; | |
747 | ||
748 | if (list->IsEmpty()) | |
749 | return 1; | |
750 | ||
751 | TIterator* iter = list->MakeIterator(); | |
752 | TObject* obj; | |
753 | ||
754 | // sub collections | |
7307d52c | 755 | const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning |
7029240a | 756 | TList* collections[nCollections]; |
757 | for (Int_t i=0; i<nCollections; ++i) | |
758 | collections[i] = new TList; | |
759 | ||
760 | Int_t count = 0; | |
761 | while ((obj = iter->Next())) | |
762 | { | |
763 | dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj); | |
764 | if (entry == 0) | |
765 | continue; | |
766 | ||
45e97e28 | 767 | collections[0]->Add(entry->fData); |
0fc41645 | 768 | collections[1]->Add(entry->fMult); |
769 | collections[2]->Add(entry->fPtDist); | |
7029240a | 770 | |
771 | for (Int_t i=0; i<kVertexBinning; ++i) | |
1afae8ff | 772 | { |
7307d52c | 773 | collections[3+i]->Add(entry->fdNdEta[i]); |
774 | collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]); | |
1afae8ff | 775 | } |
7029240a | 776 | |
777 | ++count; | |
778 | } | |
779 | ||
45e97e28 | 780 | fData->Merge(collections[0]); |
7307d52c | 781 | fMult->Merge(collections[1]); |
782 | fPtDist->Merge(collections[2]); | |
783 | ||
7029240a | 784 | for (Int_t i=0; i<kVertexBinning; ++i) |
1afae8ff | 785 | { |
7307d52c | 786 | fdNdEta[i]->Merge(collections[3+i]); |
787 | fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]); | |
1afae8ff | 788 | } |
7029240a | 789 | |
790 | for (Int_t i=0; i<nCollections; ++i) | |
791 | delete collections[i]; | |
792 | ||
793 | return count+1; | |
ceb5d1b5 | 794 | } |
567160d6 | 795 | |
796 |