Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEtofPIDqa.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15 //
16 // Class AliHFEtofPIDqa
17 //
18 // Monitoring TOF PID in the HFE PID montioring framework. The following
19 // quantities are monitored:
20 //   TOF time distance to electron hypothesis (Number of sigmas)
21 //   TOF time distance to the pion hypothesis (Absolute value)
22 //   TPC dE/dx (Number of sigmas, as control histogram)
23 // (Always as function of momentum, particle species and centrality 
24 // before and after cut)
25 // More information about the PID monitoring framework can be found in
26 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
27 //
28 // Author:
29 //
30 //   Markus Fasel <M.Fasel@gsi.de>
31 #include <TClass.h>
32 #include <TH2.h>
33 #include <THnSparse.h>
34 #include <TString.h>
35
36 #include "AliLog.h"
37 #include "AliPID.h"
38 #include "AliPIDResponse.h"
39 #include "AliVParticle.h"
40
41 #include "AliHFEcollection.h"
42 #include "AliHFEpid.h"
43 #include "AliHFEpidBase.h"
44 #include "AliHFEpidQAmanager.h"
45 #include "AliHFEpidTPC.h"
46 #include "AliHFEpidTRD.h"
47 #include "AliHFEpidTOF.h"
48 #include "AliHFEtools.h"
49 #include "AliHFEtofPIDqa.h"
50
51 ClassImp(AliHFEpidTOF)
52
53 //_________________________________________________________
54 AliHFEtofPIDqa::AliHFEtofPIDqa():
55     AliHFEdetPIDqa()
56   , fHistos(NULL)
57 {
58   //
59   // Dummy constructor
60   //
61 }
62
63 //_________________________________________________________
64 AliHFEtofPIDqa::AliHFEtofPIDqa(const char* name):
65     AliHFEdetPIDqa(name, "QA for TOF")
66   , fHistos(NULL)
67 {
68   //
69   // Default constructor
70   //
71 }
72
73 //_________________________________________________________
74 AliHFEtofPIDqa::AliHFEtofPIDqa(const AliHFEtofPIDqa &o):
75     AliHFEdetPIDqa(o)
76   , fHistos()
77 {
78   //
79   // Copy constructor
80   //
81   o.Copy(*this);
82 }
83
84 //_________________________________________________________
85 AliHFEtofPIDqa &AliHFEtofPIDqa::operator=(const AliHFEtofPIDqa &o){
86   //
87   // Do assignment
88   //
89   AliHFEdetPIDqa::operator=(o);
90   if(&o != this) o.Copy(*this);
91   return *this;
92 }
93
94 //_________________________________________________________
95 AliHFEtofPIDqa::~AliHFEtofPIDqa(){
96   //
97   // Destructor
98   //
99   if(fHistos) delete fHistos;
100 }
101
102 //_________________________________________________________
103 void AliHFEtofPIDqa::Copy(TObject &o) const {
104   //
105   // Make copy
106   //
107   AliHFEtofPIDqa &target = dynamic_cast<AliHFEtofPIDqa &>(o);
108   if(target.fHistos){
109     delete target.fHistos;
110     target.fHistos = NULL;
111   }
112   if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
113 }
114
115 //_________________________________________________________
116 Long64_t AliHFEtofPIDqa::Merge(TCollection *coll){
117   //
118   // Merge with other objects
119   //
120   if(!coll) return 0;
121   if(coll->IsEmpty()) return 1;
122
123   TIter it(coll);
124   AliHFEtofPIDqa *refQA = NULL;
125   TObject *o = NULL;
126   Long64_t count = 0;
127   TList listHistos;
128   while((o = it())){
129     refQA = dynamic_cast<AliHFEtofPIDqa *>(o);
130     if(!refQA) continue;
131
132     listHistos.Add(refQA->fHistos);
133     count++; 
134   }
135   fHistos->Merge(&listHistos);
136   return count + 1;
137 }
138
139 //_________________________________________________________
140 void AliHFEtofPIDqa::Initialize(){
141   //
142   // Define Histograms
143   //
144
145   fHistos = new AliHFEcollection("tofqahistos", "Collection of TOF QA histograms");
146
147   // Make common binning
148   const Int_t kPIDbins = AliPID::kSPECIES + 1;
149   const Int_t kSteps = 2;
150   const Int_t kCentralityBins = 11;
151   const Double_t kMinPID = -1;
152   const Double_t kMinP = 0.;
153   const Double_t kMaxPID = (Double_t)AliPID::kSPECIES;
154   const Double_t kMaxP = 20.;
155   // Quantities where one can switch between low and high resolution
156   Int_t kPbins = fQAmanager->HasHighResolutionHistos() ?  1000 : 100;
157   Int_t kSigmaBins = fQAmanager->HasHighResolutionHistos() ? 1400 : 140;
158
159   // 1st histogram: TOF sigmas: (species, p nsigma, step)
160   Int_t nBinsSigma[5] = {kPIDbins, kPbins, kSigmaBins, kSteps, kCentralityBins};
161   Double_t minSigma[5] = {kMinPID, kMinP, -12., 0, 0};
162   Double_t maxSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
163   fHistos->CreateTHnSparse("tofnSigma", "TOF signal; species; p [GeV/c]; TOF signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
164
165   // 3rd histogram: TPC sigmas to the electron line: (species, p nsigma, step - only filled if apriori PID information available)
166   fHistos->CreateTHnSparse("tofMonitorTPC", "TPC signal; species; p [GeV/c]; TPC signal [a.u.]; Selection Step; Centrality", 5, nBinsSigma, minSigma, maxSigma);
167 }
168
169 //_________________________________________________________
170 void AliHFEtofPIDqa::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step){
171   //
172   // Fill TPC histograms
173   //
174   //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
175   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
176
177   Float_t centrality = track->GetCentrality();
178   Int_t species = track->GetAbInitioPID();
179   if(species >= AliPID::kSPECIES) species = -1;
180
181   AliDebug(1, Form("Monitoring particle of type %d for step %d", species, step));
182   AliHFEpidTOF *tofpid= dynamic_cast<AliHFEpidTOF *>(fQAmanager->GetDetectorPID(AliHFEpid::kTOFpid));
183   AliHFEpidTRD *trdpid= dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
184   
185   const AliPIDResponse *pidResponse = tofpid ? tofpid->GetPIDResponse() : NULL;
186   Double_t contentSignal[5];
187   contentSignal[0] = species;
188   contentSignal[1] = track->GetRecTrack()->P();
189   contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTOF(track->GetRecTrack(), AliPID::kElectron): 0.;
190   contentSignal[3] = step;
191   contentSignal[4] = centrality;
192   fHistos->Fill("tofnSigma", contentSignal);
193   //if(species > -1){
194   contentSignal[1] =trdpid ? trdpid->GetP(track->GetRecTrack(), anatype) : 0.;
195   contentSignal[2] = pidResponse ? pidResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron) : 0.;
196     fHistos->Fill("tofMonitorTPC", contentSignal);
197   //}
198 }
199
200 //_________________________________________________________
201 TH2 *AliHFEtofPIDqa::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
202   //
203   // Get the TPC control histogram for the TRD selection step (either before or after PID)
204   //
205   THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("tofMonitorTPC"));
206   if(!histo){
207     AliError("QA histogram monitoring TPC nSigma not available");
208     return NULL;
209   }
210   if(species > -1 && species < AliPID::kSPECIES){
211     // cut on species (if available)
212     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
213   }
214   TString centname, centtitle;
215   Bool_t hasCentralityInfo = kTRUE;
216   if(centralityClass > -1){
217     if(histo->GetNdimensions() < 5){
218       AliError("Centrality Information not available");
219       centname = centtitle = "MinBias";
220       hasCentralityInfo = kFALSE;
221     } else {
222       // Project centrality classes
223       // -1 is Min. Bias
224       histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
225       centname = Form("Cent%d", centralityClass);
226       centtitle = Form("Centrality %d", centralityClass);
227     }
228   } else {
229     histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
230     centname = centtitle = "MinBias";
231     hasCentralityInfo = kTRUE;
232   }
233   histo->GetAxis(3)->SetRange(step + 1, step + 1); 
234
235   TH2 *hSpec = histo->Projection(2, 1);
236   // construct title and name
237   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
238   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
239   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
240   TString histname = Form("hSigmaTPC%s%s%s", specID.Data(), stepname.Data(), centname.Data());
241   TString histtitle = Form("TPC Sigma for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
242   hSpec->SetName(histname.Data());
243   hSpec->SetTitle(histtitle.Data());
244
245   // Unset range on the original histogram
246   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
247   histo->GetAxis(3)->SetRange(0, histo->GetAxis(3)->GetNbins());
248   if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
249   return hSpec; 
250 }
251
252 //_________________________________________________________
253 TH2 *AliHFEtofPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
254   //
255   // Plot the Spectrum
256   //
257   THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tofnSigma"));
258   if(!hSignal) return NULL;
259   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
260   if(species >= 0 && species < AliPID::kSPECIES)
261     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
262   TString centname, centtitle;
263   Bool_t hasCentralityInfo = kTRUE;
264   if(centralityClass > -1){
265     if(hSignal->GetNdimensions() < 5){
266       AliError("Centrality Information not available");
267       centname = centtitle = "MinBias";
268       hasCentralityInfo = kFALSE;
269     } else {
270       // Project centrality classes
271       // -1 is Min. Bias
272       AliInfo(Form("Centrality Class: %d", centralityClass));
273       hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
274       centname = Form("Cent%d", centralityClass);
275       centtitle = Form("Centrality %d", centralityClass);
276     }
277   } else {
278     hSignal->GetAxis(4)->SetRange(1, hSignal->GetAxis(4)->GetNbins()-1);
279     centname = centtitle = "MinBias";
280     hasCentralityInfo = kTRUE;
281   }
282   TH2 *hTmp = hSignal->Projection(2,1);
283   TString hname = Form("hTOFsigma%s%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after", centname.Data()), 
284           htitle = Form("Time-of-flight spectrum[#sigma] %s selection %s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after", centtitle.Data());
285   if(species > -1){
286     hname += AliPID::ParticleName(species);
287     htitle += Form(" for %ss", AliPID::ParticleName(species));
288   }
289   hTmp->SetName(hname.Data());
290   hTmp->SetTitle(htitle.Data());
291   hTmp->SetStats(kFALSE);
292   hTmp->GetXaxis()->SetTitle("p [GeV/c]");
293   hTmp->GetYaxis()->SetTitle("TOF time|_{el} - expected time|_{el} [#sigma]");
294   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
295   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
296   if(hasCentralityInfo)hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
297   return hTmp;
298 }
299
300 //_________________________________________________________
301 TH1 *AliHFEtofPIDqa::GetHistogram(const char *name){
302   // 
303   // Get the histogram
304   //
305   if(!fHistos) return NULL;
306   TObject *histo = fHistos->Get(name);
307   if(!histo->InheritsFrom("TH1")) return NULL;
308   return dynamic_cast<TH1 *>(histo);
309 }
310