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() :
32 // default constructor
34 for (Int_t i=0; i<kVertexBinning; ++i)
37 fdNdEtaPtCutOffCorrected[i] = 0;
41 //____________________________________________________________________
42 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
49 // TODO this binning has to be the same than in AliCorrection, somehow passed?!
50 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};
52 fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode);
54 // do not add this hists to the directory
55 Bool_t oldStatus = TH1::AddDirectoryStatus();
56 TH1::AddDirectory(kFALSE);
58 fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", 40, -2, 2);
60 fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
62 for (Int_t i=1; i<kVertexBinning; ++i)
64 fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
65 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
68 fPtDist = new TH1F("Pt", "p_{T} distribution;p_{T} [GeV/c];#frac{dN}{d#eta dp_{T}} [c/GeV]", 28, binLimitsPt);
70 TH1::AddDirectory(oldStatus);
73 //____________________________________________________________________
74 dNdEtaAnalysis::~dNdEtaAnalysis()
84 for (Int_t i=0; i<kVertexBinning; ++i)
91 if (fdNdEtaPtCutOffCorrected[i])
93 delete fdNdEtaPtCutOffCorrected[i];
94 fdNdEtaPtCutOffCorrected[i] = 0;
105 //_____________________________________________________________________________
106 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
112 // dNdEtaAnalysis copy constructor
115 ((dNdEtaAnalysis &) c).Copy(*this);
118 //_____________________________________________________________________________
119 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
122 // Assignment operator
125 if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
129 //_____________________________________________________________________________
130 void dNdEtaAnalysis::Copy(TObject &c) const
136 dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
138 target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
140 for (Int_t i=0; i<kVertexBinning; ++i)
142 target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
143 target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
146 target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
148 TNamed::Copy((TNamed &) c);
151 //____________________________________________________________________
152 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
154 // fills a track into the histograms
156 fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
159 //____________________________________________________________________
160 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
162 // fills an event into the histograms
164 fData->GetEventCorrection()->FillMeas(vtx, n);
167 //____________________________________________________________________
168 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType)
171 // correct with the given correction values and calculate dNdEta and pT distribution
172 // the corrections that are applied can be steered by the flag correctionType
175 // TODO put tag somewhere which corrections have been applied
177 // set corrections to 1
178 fData->SetCorrectionToUnity();
180 if (correction && correctionType != AlidNdEtaCorrection::kNone)
182 TH3F* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
183 TH2F* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
185 if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
186 trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
188 if (correctionType >= AlidNdEtaCorrection::kVertexReco)
190 trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
191 eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
194 switch (correctionType)
196 case AlidNdEtaCorrection::kINEL :
198 trackCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetCorrectionHistogram());
199 eventCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram());
202 case AlidNdEtaCorrection::kNSD :
204 trackCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetTrackCorrection()->GetCorrectionHistogram());
205 eventCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetEventCorrection()->GetCorrectionHistogram());
208 case AlidNdEtaCorrection::kND :
210 trackCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetTrackCorrection()->GetCorrectionHistogram());
211 eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram());
218 printf("INFO: No correction applied\n");
222 TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
224 // integrate multiplicity axis out (include under/overflow bins!!!)
225 TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
226 TH1D* vertexHist = tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
231 dataHist->GetXaxis()->SetRange(0, 0);
232 dataHist->GetYaxis()->SetRange(0, 0);
233 dataHist->GetZaxis()->SetRange(0, 0);
236 Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5);
237 Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5);
238 dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
239 Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
244 dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8));
245 Float_t etaWidth = 1.6;
247 TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
249 for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
251 Float_t binSize = fPtDist->GetBinWidth(i);
252 fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
253 fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
259 printf("ERROR: nEvents is 0!\n");
263 dataHist->GetXaxis()->SetRange(0, 0);
264 dataHist->GetYaxis()->SetRange(0, 0);
265 dataHist->GetZaxis()->SetRange(0, 0);
267 // integrate over pt (with pt cut)
270 ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
272 dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins());
273 printf("range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins());
274 TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
276 dataHist->GetZaxis()->SetRange(0, 0);
277 vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
278 vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
282 printf("ERROR: pt integration failed\n");
286 //new TCanvas; vtxVsEta->DrawCopy();
287 //vtxVsEta->Rebin2D(1, 4);
289 const Float_t vertexRange = 9.99;
291 for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
293 // do we have several histograms for different vertex positions?
294 Int_t vertexBinGlobalBegin = vertexHist->GetXaxis()->FindBin(-vertexRange);
295 Int_t vertexBinWidth = (vertexHist->GetXaxis()->FindBin(vertexRange) - vertexBinGlobalBegin + 1) / (kVertexBinning-1);
296 //printf("vertexBinGlobalBegin = %d, vertexBinWidth = %d\n", vertexBinGlobalBegin, vertexBinWidth);
297 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
300 Int_t vertexBinBegin = vertexBinGlobalBegin;
301 Int_t vertexBinEnd = vertexBinGlobalBegin + vertexBinWidth * (kVertexBinning-1);
303 // the first histogram is always for the whole vertex range
306 vertexBinBegin = vertexBinGlobalBegin + vertexBinWidth * (vertexPos-1);
307 vertexBinEnd = vertexBinBegin + vertexBinWidth;
310 //printf("vertexBinBegin = %d, vertexBinEnd = %d\n", vertexBinBegin, vertexBinEnd);
312 Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd - 1);
313 if (totalEvents == 0)
315 printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
320 Float_t sumError2 = 0;
321 for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
323 if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
325 sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
327 if (sumError2 > 10e30)
328 printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
330 sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
334 Float_t ptCutOffCorrection = 1;
335 if (correction && ptCut > 0)
336 ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
338 if (ptCutOffCorrection <= 0)
340 printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
344 //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);
346 Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
347 if (bin > 0 && bin < fdNdEta[vertexPos]->GetNbinsX())
349 Float_t dndeta = sum / totalEvents;
350 Float_t error = TMath::Sqrt(sumError2) / totalEvents;
352 dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
353 error = error / fdNdEta[vertexPos]->GetBinWidth(bin);
355 fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
356 fdNdEta[vertexPos]->SetBinError(bin, error);
358 dndeta /= ptCutOffCorrection;
359 error /= ptCutOffCorrection;
361 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
362 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
364 //Printf("Bin %d has dN/deta = %f", bin, dndeta);
369 //new TCanvas; fdNdEta[0]->DrawCopy();
372 //____________________________________________________________________
373 void dNdEtaAnalysis::SaveHistograms()
375 // save the histograms to a directory with the name of this class (retrieved from TNamed)
377 gDirectory->mkdir(GetName());
378 gDirectory->cd(GetName());
382 fData->SaveHistograms();
385 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
390 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
392 for (Int_t i=0; i<kVertexBinning; ++i)
397 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
399 if (fdNdEtaPtCutOffCorrected[i])
400 fdNdEtaPtCutOffCorrected[i]->Write();
402 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
405 gDirectory->cd("../");
408 void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
410 // loads the histograms from a directory with the name of this class (retrieved from TNamed)
417 fData->LoadHistograms();
419 for (Int_t i=0; i<kVertexBinning; ++i)
421 fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
422 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
425 fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
427 gDirectory->cd("../");
430 //____________________________________________________________________
431 void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
433 // draws the histograms
438 fData->DrawHistograms(GetName());
440 TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
441 canvas->Divide(2, 1);
444 if (fdNdEtaPtCutOffCorrected[0])
445 fdNdEtaPtCutOffCorrected[0]->DrawCopy();
449 fdNdEta[0]->SetLineColor(kRed);
450 fdNdEta[0]->DrawCopy("SAME");
458 // histograms for different vertices?
459 if (kVertexBinning > 0)
461 // doesnt work, but i dont get it, giving up...
462 TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450);
463 TCanvas* canvas3 = 0;
465 canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450);
467 //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
468 //printf("%d\n", yPads);
469 //canvas2->Divide(2, yPads);
471 TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
473 for (Int_t i=0; i<kVertexBinning; ++i)
475 if (fdNdEtaPtCutOffCorrected[i])
479 fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
480 fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
481 legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
483 if (canvas3 && fdNdEta[i])
487 fdNdEta[i]->SetLineColor(i+1);
488 fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
494 canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
503 if (kVertexBinning == 3)
505 TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone"));
506 TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2"));
510 TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
512 clone->Divide(fdNdEtaPtCutOffCorrected[0]);
513 clone->GetYaxis()->SetRangeUser(0.95, 1.05);
516 clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
517 clone2->DrawCopy("SAME");
519 TLine* line = new TLine(-1, 1, 1, 1);
522 canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName()));
527 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
529 // Merges a list of dNdEtaAnalysis objects with this one.
530 // This is needed for PROOF.
531 // Returns the number of merged objects (including this)
539 TIterator* iter = list->MakeIterator();
543 const Int_t nCollections = 2 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
544 TList* collections[nCollections];
545 for (Int_t i=0; i<nCollections; ++i)
546 collections[i] = new TList;
549 while ((obj = iter->Next()))
551 dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
555 collections[0]->Add(entry->fData);
556 collections[1]->Add(entry->fPtDist);
558 for (Int_t i=0; i<kVertexBinning; ++i)
560 collections[2+i]->Add(entry->fdNdEta[i]);
561 collections[2+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
567 fData->Merge(collections[0]);
568 fPtDist->Merge(collections[1]);
569 for (Int_t i=0; i<kVertexBinning; ++i)
571 fdNdEta[i]->Merge(collections[2+i]);
572 fdNdEta[i]->Merge(collections[2+kVertexBinning+i]);
575 for (Int_t i=0; i<nCollections; ++i)
576 delete collections[i];