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 // AliSpectraBothHistoManager 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 "AliSpectraBothHistoManager.h"
41 ClassImp(AliSpectraBothHistoManager)
44 using namespace AliSpectraNameSpaceBoth;
45 #include "HistogramNamesBoth.cxx" // generate this automatically running createNames.py
48 //const char* kParticleSpecies[] =
59 AliSpectraBothHistoManager::AliSpectraBothHistoManager(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(kHistNameBoth[ihist]); // PT histos
70 if (ihist > kNPtGenHist && ihist <= kNPtGenAllChHist) BookPtGenAllChHistogram(kHistNameBoth[ihist]); // PT histos
71 if (ihist > kNPtGenAllChHist && ihist <= kNPtRecHist) BookPtRecHistogram(kHistNameBoth[ihist]); // PT histos
72 if (ihist > kNPtRecHist && ihist <= kNPtRecAllChHist) BookPtRecAllChHistogram(kHistNameBoth[ihist]); // PT histos
73 if (ihist > kNPtRecAllChHist && ihist <= kNHistPID) BookPIDHistogram(kHistNameBoth[ihist]); // PID histos
74 if (ihist > kNHistPID && ihist <= kNHistNSig) BookNSigHistogram(kHistNameBoth[ihist]); // NSigmaSep histos
75 if(ihist==kHistGenMulvsRawMul) BookGenMulvsRawMulHistogram(kHistNameBoth[ihist]);
78 TH1::AddDirectory(oldStatus);
80 //_______________________________________________________________________________________
81 AliSpectraBothHistoManager::~AliSpectraBothHistoManager()
85 fOutputList->Delete();
89 //_______________________________________________________
91 TH2F* AliSpectraBothHistoManager::BookPtGenHistogram(const char * name)
93 // Return a pt histogram with predefined binning, set the ID and add it to the output list
94 AliInfo(Form("Booking pt gen histogram %s", name));
97 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};
100 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)
101 hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
102 hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
103 hist->SetMarkerStyle(kFullCircle);
105 fOutputList->Add(hist);
110 //_______________________________________________________
112 TH2F* AliSpectraBothHistoManager::BookPtGenAllChHistogram(const char * name)
114 // Return a pt histogram with predefined binning, set the ID and add it to the output list
115 AliInfo(Form("Booking pt gen histogram - no PID %s", name));
118 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};
121 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)
122 hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
123 hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
124 hist->SetMarkerStyle(kFullCircle);
126 fOutputList->Add(hist);
132 //_______________________________________________________
133 TH2F* AliSpectraBothHistoManager::BookPtRecHistogram(const char * name)
135 // Return a pt histogram with predefined binning, set the ID and add it to the output list
136 AliInfo(Form("Booking pt rec histogram %s, rebin:%d", name, fNRebin));
139 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};
142 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)
143 if(fNRebin!=0)hist->RebinY(fNRebin);
144 hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
145 hist->GetYaxis()->SetTitle("DCA xy");
146 hist->SetMarkerStyle(kFullCircle);
148 fOutputList->Add(hist);
153 //_______________________________________________________
154 TH2F* AliSpectraBothHistoManager::BookPtRecAllChHistogram(const char * name)
156 // Return a pt histogram with predefined binning, set the ID and add it to the output list
157 AliInfo(Form("Booking pt rec histogram %s, rebin:%d", name, fNRebin));
160 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};
163 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)
164 if(fNRebin!=0)hist->RebinY(fNRebin);
165 hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
166 hist->GetYaxis()->SetTitle("DCA xy");
167 hist->SetMarkerStyle(kFullCircle);
169 fOutputList->Add(hist);
174 //_____________________________________________________________________________
176 TH2F* AliSpectraBothHistoManager::BookPIDHistogram(const char * name)
178 // Return a pt histogram with predefined binning, set the ID and add it to the output list
179 AliInfo(Form("Booking PID histogram %s, rebin:%d", name, fNRebin));
181 TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 50, 0, 2.5, 200, -1000, 1000);
183 hist->RebinY(fNRebin);
184 // hist->RebinX(fNRebin);
186 if(tmp.Contains("Pt"))
187 hist->GetXaxis()->SetTitle("P_{t} (GeV / c)");
189 hist->GetXaxis()->SetTitle("P (GeV / c)");
191 hist->GetYaxis()->SetTitle("PID signal");
193 fOutputList->Add(hist);
198 //_____________________________________________________________________________
200 TH2F* AliSpectraBothHistoManager::BookNSigHistogram(const char * name)
202 // Return a pt histogram with predefined binning, set the ID and add it to the output list
203 AliInfo(Form("Booking NSigma histogram %s, rebin:%d", name, fNRebin));
207 if(tmp.Contains("TPCTOF"))
212 TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 50, 0, 2.5,nbins,miny, 20);
214 hist->RebinY(fNRebin);
215 //hist->RebinX(fNRebin);
217 if(tmp.Contains("Pt"))
218 hist->GetXaxis()->SetTitle("P_{t} (GeV / c)");
220 hist->GetXaxis()->SetTitle("P (GeV / c)");
221 //hist->GetYaxis()->SetTitle("TPC");
223 fOutputList->Add(hist);
228 //_____________________________________________________________________________
230 TH2F* AliSpectraBothHistoManager::BookGenMulvsRawMulHistogram(const char * name)
232 // Return a pt histogram with predefined binning, set the ID and add it to the output list
233 AliInfo(Form("Booking Gen vs Raw multilicity histogram %s", name));
235 TH2F * hist = new TH2F(name, Form("Gen vs Raw multilicity (%s)", name), 100, -0.5, 99.5, 100, -0.5, 99.5);
236 hist->GetXaxis()->SetTitle("gen mul ");
237 hist->GetYaxis()->SetTitle("raw mul ");
239 fOutputList->Add(hist);
243 //_____________________________________________________________________________
245 TH1F* AliSpectraBothHistoManager::GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA)
247 // //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
248 // //if minDCA=-1 && maxDCA=-1 the projection is done using the full DCA range
249 TH2F *hist=(TH2F*)fOutputList->FindObject(name);
251 AliDebug(2,Form("--- Projecting %s on Xaxis[%f,%f]:",name,minDCA,maxDCA));
252 if(minDCA==-1 && maxDCA==-1){
253 outhist=(TH1F*)hist->ProjectionX("_px",0,-1,"e");
254 AliDebug(2,"Full Range");
256 Int_t firstbin=hist->GetYaxis()->FindBin(minDCA);
257 Int_t lastbin=hist->GetYaxis()->FindBin(maxDCA);
258 AliDebug(2,Form("firstbin: %d lastbin: %d",firstbin,lastbin));
259 outhist=(TH1F*)hist->ProjectionX("_px",firstbin,lastbin,"e");
261 AliDebug(2,Form("Entries outhist: %.0f Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries()));
265 //_____________________________________________________________________________
267 TH1F* AliSpectraBothHistoManager::GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt)
269 // //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
270 // //if minPt=-1 && maxPt=-1 the projection is done using the full DCA range
271 TH2F *hist=(TH2F*)fOutputList->FindObject(name);
273 AliDebug(2,Form("--- Projecting %s on Yaxis[%f,%f]:",name,minPt,maxPt));
274 if(minPt==-1 && maxPt==-1){
275 outhist=(TH1F*)hist->ProjectionY("_py",0,-1,"e");
276 AliDebug(2,"Full Range");
278 Int_t firstbin=hist->GetXaxis()->FindBin(minPt);
279 Int_t lastbin=hist->GetXaxis()->FindBin(maxPt);
280 AliDebug(2,Form("firstbin: %d lastbin: %d",firstbin,lastbin));
281 outhist=(TH1F*)hist->ProjectionY("_py",firstbin,lastbin,"e");
282 AliDebug(2,Form("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)));
284 AliDebug(2,Form("Entries outhist: %.0f Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries()));
288 //_____________________________________________________________________________
290 Long64_t AliSpectraBothHistoManager::Merge(TCollection* list)
292 // Merge a list of AliSpectraBothHistoManager objects with this.
293 // Returns the number of merged objects (including this).
295 // AliInfo("Merging");
303 TIterator* iter = list->MakeIterator();
306 // collections of all histograms
311 while ((obj = iter->Next())) {
312 AliSpectraBothHistoManager* entry = dynamic_cast<AliSpectraBothHistoManager*> (obj);
316 TList * hlist = entry->GetOutputList();
317 collections.Add(hlist);
321 fOutputList->Merge(&collections);
328 //________________________________________________________________________________
329 TH1* AliSpectraBothHistoManager::GetHistogram1D(UInt_t histoType, UInt_t particleType, UInt_t charge) {
330 // GetHistogram using particle ID and histogram type
333 if (particleType == kSpUndefined) {
334 AliError ("Trying to get histo for undefined particle");
339 case kHistPtGenTruePrimary:
340 baseId = kHistPtGenTruePrimaryPionPlus;
342 case kHistPtRecSigma:
343 baseId = kHistPtRecSigmaPionPlus;
346 baseId = kHistPtRecTruePionPlus;
348 case kHistPtRecTruePrimary:
349 baseId = kHistPtRecTruePrimaryPionPlus;
351 case kHistPtRecPrimary:
352 baseId = kHistPtRecPrimaryPionPlus;
354 case kHistPtRecSigmaPrimary:
355 baseId = kHistPtRecSigmaPrimaryPionPlus;
357 case kHistPtRecSigmaSecondaryMaterial:
358 baseId = kHistPtRecSigmaSecondaryMaterialPionPlus;
360 case kHistPtRecSigmaSecondaryWeakDecay:
361 baseId = kHistPtRecSigmaSecondaryWeakDecayPionPlus;
364 baseId = kHistNSigPionTPC;
367 baseId = kHistNSigPionTOF;
369 case kHistNSigTPCTOF:
370 baseId = kHistNSigPionTPCTOF;
377 AliFatal(Form("Wrong histogram type %d", histoType));
379 //cout << "T[" << histoType << "] ID["<< baseId <<"] P["<<particleType<<"] C[" << charge
380 // << " --> ["<< baseId + particleType + 3*(charge) <<"] = " ;
382 baseId = baseId + particleType + 3*(charge);
384 //cout << GetHistogram(baseId)->GetName() << endl;
386 return GetHistogram(baseId);
388 //____________________________________________________________________________________________________
389 TH2* AliSpectraBothHistoManager::GetHistogram2D(UInt_t histoType, UInt_t particleType, UInt_t charge){
390 // returns histo based on ids, casting it to TH2*
391 return (TH2*) GetHistogram1D(histoType,particleType,charge);