2 /**************************************************************************
3 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
17 //-----------------------------------------------------------------
18 // AliSpectraAODHistoManager class
19 //-----------------------------------------------------------------
28 #include "AliAnalysisTask.h"
29 #include "AliAnalysisManager.h"
30 #include "AliAODTrack.h"
31 #include "AliAODMCParticle.h"
32 #include "AliAODEvent.h"
33 #include "AliAODInputHandler.h"
34 #include "AliAnalysisTaskESDfilter.h"
35 #include "AliAnalysisDataContainer.h"
36 #include "AliSpectraAODHistoManager.h"
42 ClassImp(AliSpectraAODHistoManager)
45 using namespace AliSpectraNameSpace;
46 #include "HistogramNames.cxx" // generate this automatically running createNames.py
48 const char * kParticleSpecies[] =
59 AliSpectraAODHistoManager::AliSpectraAODHistoManager(const char *name,Int_t nrebin): TNamed(name, "AOD Spectra Histo Manager"), fOutputList(0), fNRebin(0)
63 fOutputList = new TList;
64 fOutputList->SetOwner();
65 Bool_t oldStatus = TH1::AddDirectoryStatus();
66 TH1::AddDirectory(kFALSE);
67 for (Int_t ihist = 0; ihist < kNHist ; ihist++)
69 if (ihist <= kNPtGenHist) BookPtGenHistogram(kHistName[ihist]); // PT histos
70 if (ihist > kNPtGenHist && ihist <= kNPtGenAllChHist) BookPtGenAllChHistogram(kHistName[ihist]); // PT histos
71 if (ihist > kNPtGenAllChHist && ihist <= kNPtRecHist) BookPtRecHistogram(kHistName[ihist]); // PT histos
72 if (ihist > kNPtRecHist && ihist <= kNPtRecAllChHist) BookPtRecAllChHistogram(kHistName[ihist]); // PT histos
73 if (ihist > kNPtRecAllChHist && ihist <= kNHistPID) BookPIDHistogram(kHistName[ihist]); // PID histos
74 if (ihist > kNHistPID && ihist <= kNHistNSig) BookNSigHistogram(kHistName[ihist]); // NSigmaSep histos
77 TH1::AddDirectory(oldStatus);
80 //_______________________________________________________
82 TH2F* AliSpectraAODHistoManager::BookPtGenHistogram(const char * name)
84 // Return a pt histogram with predefined binning, set the ID and add it to the output list
85 AliInfo(Form("Booking pt gen histogram %s", name));
88 const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
91 TH2F * hist = new TH2F(name,Form("P_{T} distribution (%s)", name),nbinsTempl,templBins,2,-0.5,1.5);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
92 hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
93 hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
94 hist->SetMarkerStyle(kFullCircle);
96 fOutputList->Add(hist);
101 //_______________________________________________________
103 TH2F* AliSpectraAODHistoManager::BookPtGenAllChHistogram(const char * name)
105 // Return a pt histogram with predefined binning, set the ID and add it to the output list
106 AliInfo(Form("Booking pt gen histogram - no PID %s", name));
109 const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7,7.5,8,8.5,9,9.5,10};
112 TH2F * hist = new TH2F(name,Form("P_{T} distribution (All Ch) (%s)", name),nbinsTempl,templBins,2,-0.5,1.5);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
113 hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
114 hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
115 hist->SetMarkerStyle(kFullCircle);
117 fOutputList->Add(hist);
123 //_______________________________________________________
124 TH2F* AliSpectraAODHistoManager::BookPtRecHistogram(const char * name)
126 // Return a pt histogram with predefined binning, set the ID and add it to the output list
127 AliInfo(Form("Booking pt rec histogram %s, rebin:%d", name, fNRebin));
130 const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
133 TH2F * hist = new TH2F(name,Form("reconstructed P_{T} distribution (%s)", name),nbinsTempl,templBins,3000,-3,3);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
134 if(fNRebin!=0)hist->RebinY(fNRebin);
135 hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
136 hist->GetYaxis()->SetTitle("DCA xy");
137 hist->SetMarkerStyle(kFullCircle);
139 fOutputList->Add(hist);
144 //_______________________________________________________
145 TH2F* AliSpectraAODHistoManager::BookPtRecAllChHistogram(const char * name)
147 // Return a pt histogram with predefined binning, set the ID and add it to the output list
148 AliInfo(Form("Booking pt rec histogram %s, rebin:%d", name, fNRebin));
151 const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7,7.5,8,8.5,9,9.5,10};
154 TH2F * hist = new TH2F(name,Form("reconstructed P_{T} distribution (All Ch) (%s)", name),nbinsTempl,templBins,3000,-3,3);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
155 if(fNRebin!=0)hist->RebinY(fNRebin);
156 hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
157 hist->GetYaxis()->SetTitle("DCA xy");
158 hist->SetMarkerStyle(kFullCircle);
160 fOutputList->Add(hist);
165 //_____________________________________________________________________________
167 TH2F* AliSpectraAODHistoManager::BookPIDHistogram(const char * name)
169 // Return a pt histogram with predefined binning, set the ID and add it to the output list
170 AliInfo(Form("Booking PID histogram %s, rebin:%d", name, fNRebin));
172 TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 1000, 0, 2.5, 2000, -1000, 1000);
174 hist->RebinY(fNRebin);
175 hist->RebinX(fNRebin);
177 hist->GetXaxis()->SetTitle("(GeV / c)");
178 hist->GetYaxis()->SetTitle("PID signal");
180 fOutputList->Add(hist);
185 //_____________________________________________________________________________
187 TH2F* AliSpectraAODHistoManager::BookNSigHistogram(const char * name)
189 // Return a pt histogram with predefined binning, set the ID and add it to the output list
190 AliInfo(Form("Booking NSigma histogram %s, rebin:%d", name, fNRebin));
192 TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 1000, 0, 2.5, 2000,-40, 40);
194 hist->RebinY(fNRebin);
195 hist->RebinX(fNRebin);
197 hist->GetXaxis()->SetTitle("P (GeV / c)");
198 hist->GetYaxis()->SetTitle("TPC");
200 fOutputList->Add(hist);
205 //_____________________________________________________________________________
207 TH1F* AliSpectraAODHistoManager::GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA)
209 // //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
210 // //if minDCA=-1 && maxDCA=-1 the projection is done using the full DCA range
211 TH2F *hist=(TH2F*)fOutputList->FindObject(name);
213 Printf("--- Projecting %s on Xaxis[%f,%f]:",name,minDCA,maxDCA);
214 if(minDCA==-1 && maxDCA==-1){
215 outhist=(TH1F*)hist->ProjectionX("_px",0,-1,"e");
216 Printf("Full Range");
218 Int_t firstbin=hist->GetYaxis()->FindBin(minDCA);
219 Int_t lastbin=hist->GetYaxis()->FindBin(maxDCA);
220 Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
221 outhist=(TH1F*)hist->ProjectionX("_px",firstbin,lastbin,"e");
223 Printf("Entries outhist: %.0f Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
227 //_____________________________________________________________________________
229 TH1F* AliSpectraAODHistoManager::GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt)
231 // //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
232 // //if minPt=-1 && maxPt=-1 the projection is done using the full DCA range
233 TH2F *hist=(TH2F*)fOutputList->FindObject(name);
235 Printf("--- Projecting %s on Yaxis[%f,%f]:",name,minPt,maxPt);
236 if(minPt==-1 && maxPt==-1){
237 outhist=(TH1F*)hist->ProjectionY("_py",0,-1,"e");
238 Printf("Full Range");
240 Int_t firstbin=hist->GetXaxis()->FindBin(minPt);
241 Int_t lastbin=hist->GetXaxis()->FindBin(maxPt);
242 Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
243 outhist=(TH1F*)hist->ProjectionY("_py",firstbin,lastbin,"e");
244 Printf("GetDCAHistogram1D(%s) BinRange:%d %d Pt Range: %f %f",hist->GetName(),firstbin,lastbin,hist->GetXaxis()->GetBinLowEdge(firstbin),hist->GetXaxis()->GetBinLowEdge(firstbin)+hist->GetXaxis()->GetBinWidth(lastbin));
246 Printf("Entries outhist: %.0f Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
250 //_____________________________________________________________________________
252 Long64_t AliSpectraAODHistoManager::Merge(TCollection* list)
254 // Merge a list of AliSpectraAODHistoManager objects with this.
255 // Returns the number of merged objects (including this).
257 // AliInfo("Merging");
265 TIterator* iter = list->MakeIterator();
268 // collections of all histograms
273 while ((obj = iter->Next())) {
274 AliSpectraAODHistoManager* entry = dynamic_cast<AliSpectraAODHistoManager*> (obj);
278 TList * hlist = entry->GetOutputList();
279 collections.Add(hlist);
283 fOutputList->Merge(&collections);
291 TH1* AliSpectraAODHistoManager::GetHistogram1D(UInt_t histoType, UInt_t particleType, UInt_t charge) {
292 // GetHistogram using particle ID and histogram type
295 if (particleType == kSpUndefined) {
296 AliError ("Trying to get histo for undefined particle");
301 case kHistPtGenTruePrimary:
302 baseId = kHistPtGenTruePrimaryPionPlus;
304 case kHistPtRecSigma:
305 baseId = kHistPtRecSigmaPionPlus;
308 baseId = kHistPtRecTruePionPlus;
310 case kHistPtRecTruePrimary:
311 baseId = kHistPtRecTruePrimaryPionPlus;
313 case kHistPtRecPrimary:
314 baseId = kHistPtRecPrimaryPionPlus;
316 case kHistPtRecSigmaPrimary:
317 baseId = kHistPtRecSigmaPrimaryPionPlus;
319 case kHistPtRecSigmaSecondaryMaterial:
320 baseId = kHistPtRecSigmaSecondaryMaterialPionPlus;
322 case kHistPtRecSigmaSecondaryWeakDecay:
323 baseId = kHistPtRecSigmaSecondaryWeakDecayPionPlus;
326 baseId = kHistNSigPionTPC;
329 baseId = kHistNSigPionTOF;
331 case kHistNSigTPCTOF:
332 baseId = kHistNSigPionTPCTOF;
339 AliFatal(Form("Wrong histogram type %d", histoType));
341 //cout << "T[" << histoType << "] ID["<< baseId <<"] P["<<particleType<<"] C[" << charge
342 // << " --> ["<< baseId + particleType + 3*(charge) <<"] = " ;
344 baseId = baseId + particleType + 3*(charge);
346 //cout << GetHistogram(baseId)->GetName() << endl;
348 return GetHistogram(baseId);
351 TH2* AliSpectraAODHistoManager::GetHistogram2D(UInt_t histoType, UInt_t particleType, UInt_t charge){
352 // returns histo based on ids, casting it to TH2*
353 return (TH2*) GetHistogram1D(histoType,particleType,charge);