3 #include "dNdEtaAnalysis.h"
11 #include <TCollection.h>
12 #include <TIterator.h>
17 #include "AlidNdEtaCorrection.h"
18 #include <AliCorrection.h>
19 #include <AliPWG0Helper.h>
20 #include <AliCorrectionMatrix2D.h>
21 #include <AliCorrectionMatrix3D.h>
23 //____________________________________________________________________
24 ClassImp(dNdEtaAnalysis)
26 //____________________________________________________________________
27 dNdEtaAnalysis::dNdEtaAnalysis() :
33 // default constructor
35 for (Int_t i=0; i<kVertexBinning; ++i)
38 fdNdEtaPtCutOffCorrected[i] = 0;
42 //____________________________________________________________________
43 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
51 // TODO this binning has to be the same than in AliCorrection, somehow passed?!
52 Float_t binLimitsPt[] = {0.0, 0.05, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
54 fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode);
56 // do not add this hists to the directory
57 Bool_t oldStatus = TH1::AddDirectoryStatus();
58 TH1::AddDirectory(kFALSE);
60 fMult = new TH1F("TriggeredMultiplicity", "Triggered Events;raw multiplicity;entries", 1000, -0.5, 999.5);
62 fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", 40, -2, 2);
64 fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
66 for (Int_t i=1; i<kVertexBinning; ++i)
68 fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
69 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
72 fPtDist = new TH1F("Pt", "p_{T} distribution;p_{T} [GeV/c];#frac{dN}{d#eta dp_{T}} [c/GeV]", 28, binLimitsPt);
74 TH1::AddDirectory(oldStatus);
77 //____________________________________________________________________
78 dNdEtaAnalysis::~dNdEtaAnalysis()
94 for (Int_t i=0; i<kVertexBinning; ++i)
101 if (fdNdEtaPtCutOffCorrected[i])
103 delete fdNdEtaPtCutOffCorrected[i];
104 fdNdEtaPtCutOffCorrected[i] = 0;
115 //_____________________________________________________________________________
116 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
123 // dNdEtaAnalysis copy constructor
126 ((dNdEtaAnalysis &) c).Copy(*this);
129 //_____________________________________________________________________________
130 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
133 // Assignment operator
136 if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
140 //_____________________________________________________________________________
141 void dNdEtaAnalysis::Copy(TObject &c) const
147 dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
149 target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
150 target.fMult = dynamic_cast<TH1F*> (fMult->Clone());
152 for (Int_t i=0; i<kVertexBinning; ++i)
154 target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
155 target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
158 target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
160 TNamed::Copy((TNamed &) c);
163 //____________________________________________________________________
164 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
166 // fills a track into the histograms
168 fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
171 //____________________________________________________________________
172 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
174 // fills an event into the histograms
176 fData->GetEventCorrection()->FillMeas(vtx, n);
179 //____________________________________________________________________
180 void dNdEtaAnalysis::FillTriggeredEvent(Float_t n)
182 // fills a triggered event into the histograms
187 //____________________________________________________________________
188 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, Int_t multCut)
191 // correct with the given correction values and calculate dNdEta and pT distribution
192 // the corrections that are applied can be steered by the flag correctionType
193 // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
196 // TODO put tag somewhere which corrections have been applied
200 Printf("ERROR: A bigger multiplicity cut than 1 is not possible in the current implementation");
204 // set corrections to 1
205 fData->SetCorrectionToUnity();
207 if (correction && correctionType != AlidNdEtaCorrection::kNone)
209 TH3F* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
210 TH2F* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
212 if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
213 trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
215 if (correctionType >= AlidNdEtaCorrection::kVertexReco)
217 trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
218 eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
220 // set bin with multiplicity 0 to 1 (correction has no meaning in this bin)
221 for (Int_t i=0; i<=eventCorr->GetNbinsX()+1; i++)
222 eventCorr->SetBinContent(i, 1, 1);
225 switch (correctionType)
227 case AlidNdEtaCorrection::kINEL :
229 trackCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetCorrectionHistogram());
230 eventCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram());
233 case AlidNdEtaCorrection::kNSD :
235 trackCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetTrackCorrection()->GetCorrectionHistogram());
236 eventCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetEventCorrection()->GetCorrectionHistogram());
239 case AlidNdEtaCorrection::kND :
241 trackCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetTrackCorrection()->GetCorrectionHistogram());
242 eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram());
249 printf("INFO: No correction applied\n");
253 if (correctionType >= AlidNdEtaCorrection::kVertexReco)
255 // There are no events with vertex that have 0 multiplicity, therefore
256 // populate bin with 0 multiplicity with the following idea:
257 // alpha = triggered events with vertex at a given vertex position / all triggered events with vertex
258 // triggered events without vertex and 0 multiplicity at a given vertex position = alpha * all triggered events with 0 multiplicity
259 // afterwards we still correct for the trigger efficiency
261 TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
262 TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
264 //new TCanvas; correctedEvents->DrawCopy("TEXT");
266 // start above 0 mult. bin with integration
267 TH1* vertexDist = measuredEvents->ProjectionX("vertexdist_measured", 2);
268 //new TCanvas; vertexDist->DrawCopy();
270 Int_t allEventsWithVertex = (Int_t) vertexDist->Integral(0, vertexDist->GetNbinsX()+1); // include under/overflow!
271 Int_t triggeredEventsWith0Mult = (Int_t) fMult->GetBinContent(1);
273 Printf("%d triggered events with 0 mult. -- %d triggered events with vertex", triggeredEventsWith0Mult, allEventsWithVertex);
275 for (Int_t i = 1; i <= measuredEvents->GetNbinsX(); i++)
277 Double_t alpha = vertexDist->GetBinContent(i) / allEventsWithVertex;
278 Double_t events = alpha * triggeredEventsWith0Mult;
280 // multiply with trigger correction if set above
281 events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
283 Printf("Bin %d, alpha is %.2f, number of events with 0 mult. are %.2f", i, alpha, events);
285 correctedEvents->SetBinContent(i, 1, events);
288 //new TCanvas; correctedEvents->DrawCopy("TEXT");
291 fData->PrintInfo(ptCut);
293 TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
295 // integrate multiplicity axis out (include under/overflow bins!!!)
296 TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
298 TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
303 dataHist->GetXaxis()->SetRange(0, 0);
304 dataHist->GetYaxis()->SetRange(0, 0);
305 dataHist->GetZaxis()->SetRange(0, 0);
308 Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5);
309 Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5);
310 dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
311 Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
316 dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8));
317 Float_t etaWidth = 1.6;
319 TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
321 for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
323 Float_t binSize = fPtDist->GetBinWidth(i);
324 fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
325 fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
331 printf("ERROR: nEvents is 0!\n");
335 dataHist->GetXaxis()->SetRange(0, 0);
336 dataHist->GetYaxis()->SetRange(0, 0);
337 dataHist->GetZaxis()->SetRange(0, 0);
339 // integrate over pt (with pt cut)
342 ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
344 dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
345 printf("pt range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
346 TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
348 dataHist->GetZaxis()->SetRange(0, 0);
349 vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
350 vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
354 printf("ERROR: pt integration failed\n");
358 //new TCanvas; vtxVsEta->DrawCopy("COLZ");
359 //vtxVsEta->Rebin2D(1, 4);
361 const Float_t vertexRange = 9.99;
363 for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
365 // do we have several histograms for different vertex positions?
366 Int_t vertexBinGlobalBegin = vertexHist->GetXaxis()->FindBin(-vertexRange);
367 Int_t vertexBinWidth = (vertexHist->GetXaxis()->FindBin(vertexRange) - vertexBinGlobalBegin + 1) / (kVertexBinning-1);
368 //printf("vertexBinGlobalBegin = %d, vertexBinWidth = %d\n", vertexBinGlobalBegin, vertexBinWidth);
369 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
371 Int_t vertexBinBegin = vertexBinGlobalBegin;
372 Int_t vertexBinEnd = vertexBinGlobalBegin + vertexBinWidth * (kVertexBinning-1);
374 // the first histogram is always for the whole vertex range
377 vertexBinBegin = vertexBinGlobalBegin + vertexBinWidth * (vertexPos-1);
378 vertexBinEnd = vertexBinBegin + vertexBinWidth;
381 //printf("vertexBinBegin = %d, vertexBinEnd = %d\n", vertexBinBegin, vertexBinEnd);
383 Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd - 1);
384 if (totalEvents == 0)
386 printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
391 Float_t sumError2 = 0;
392 for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
394 if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
396 sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
398 if (sumError2 > 10e30)
399 printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
401 sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
405 Float_t ptCutOffCorrection = 1;
406 if (correction && ptCut > 0)
407 ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
409 if (ptCutOffCorrection <= 0)
411 printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
415 //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);
417 Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
418 if (bin > 0 && bin < fdNdEta[vertexPos]->GetNbinsX())
420 Float_t dndeta = sum / totalEvents;
421 Float_t error = TMath::Sqrt(sumError2) / totalEvents;
423 dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
424 error = error / fdNdEta[vertexPos]->GetBinWidth(bin);
426 fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
427 fdNdEta[vertexPos]->SetBinError(bin, error);
429 dndeta /= ptCutOffCorrection;
430 error /= ptCutOffCorrection;
432 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
433 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
435 //Printf("Bin %d has dN/deta = %f", bin, dndeta);
440 //new TCanvas; fdNdEta[0]->DrawCopy();
443 //____________________________________________________________________
444 void dNdEtaAnalysis::SaveHistograms()
446 // save the histograms to a directory with the name of this class (retrieved from TNamed)
448 gDirectory->mkdir(GetName());
449 gDirectory->cd(GetName());
453 fData->SaveHistograms();
456 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
463 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n");
468 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
470 for (Int_t i=0; i<kVertexBinning; ++i)
475 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
477 if (fdNdEtaPtCutOffCorrected[i])
478 fdNdEtaPtCutOffCorrected[i]->Write();
480 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
483 gDirectory->cd("../");
486 void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
488 // loads the histograms from a directory with the name of this class (retrieved from TNamed)
495 fData->LoadHistograms();
496 fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName()));
498 for (Int_t i=0; i<kVertexBinning; ++i)
500 fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
501 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
504 fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
506 gDirectory->cd("../");
509 //____________________________________________________________________
510 void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
512 // draws the histograms
517 fData->DrawHistograms(GetName());
519 TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
520 canvas->Divide(2, 1);
523 if (fdNdEtaPtCutOffCorrected[0])
524 fdNdEtaPtCutOffCorrected[0]->DrawCopy();
528 fdNdEta[0]->SetLineColor(kRed);
529 fdNdEta[0]->DrawCopy("SAME");
537 // histograms for different vertices?
538 if (kVertexBinning > 0)
540 // doesnt work, but i dont get it, giving up...
541 TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450);
542 TCanvas* canvas3 = 0;
544 canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450);
546 //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
547 //printf("%d\n", yPads);
548 //canvas2->Divide(2, yPads);
550 TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
552 for (Int_t i=0; i<kVertexBinning; ++i)
554 if (fdNdEtaPtCutOffCorrected[i])
558 fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
559 fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
560 legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
562 if (canvas3 && fdNdEta[i])
566 fdNdEta[i]->SetLineColor(i+1);
567 fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
573 canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
582 if (kVertexBinning == 3)
584 TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone"));
585 TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2"));
589 TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
591 clone->Divide(fdNdEtaPtCutOffCorrected[0]);
592 clone->GetYaxis()->SetRangeUser(0.95, 1.05);
595 clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
596 clone2->DrawCopy("SAME");
598 TLine* line = new TLine(-1, 1, 1, 1);
601 canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName()));
606 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
608 // Merges a list of dNdEtaAnalysis objects with this one.
609 // This is needed for PROOF.
610 // Returns the number of merged objects (including this)
618 TIterator* iter = list->MakeIterator();
622 const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning
623 TList* collections[nCollections];
624 for (Int_t i=0; i<nCollections; ++i)
625 collections[i] = new TList;
628 while ((obj = iter->Next()))
630 dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
634 collections[0]->Add(entry->fData);
635 collections[2]->Add(entry->fMult);
636 collections[3]->Add(entry->fPtDist);
638 for (Int_t i=0; i<kVertexBinning; ++i)
640 collections[3+i]->Add(entry->fdNdEta[i]);
641 collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
647 fData->Merge(collections[0]);
648 fMult->Merge(collections[1]);
649 fPtDist->Merge(collections[2]);
651 for (Int_t i=0; i<kVertexBinning; ++i)
653 fdNdEta[i]->Merge(collections[3+i]);
654 fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]);
657 for (Int_t i=0; i<nCollections; ++i)
658 delete collections[i];