]> git.uio.no Git - u/mrichter/AliRoot.git/blame - CORRFW/AliCFTrackCutPid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Moved calibration here and require coincidence of HG and LG channels
[u/mrichter/AliRoot.git] / CORRFW / AliCFTrackCutPid.cxx
CommitLineData
563113d0 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//This class is intended to provide a selection on
17//the PID for single charged particles as electrons, muons
18//pions, kaons and protons. The user is supposed to set only one particle
19//to look at. The class at the moment uses one single ESD track.
20//The identification is done in Identify(), whereas the GetID() checks
21//the track status or if the combined PID should be applied.
22//The Correction Framework follows the Annalysis framework.
23//The main method of this class is the IsSelected function which returns
24//true (false) if the ESD track is (not) identified as the previously
25//setted particle.
26//
27//
28//usage:
29//
30// -----ID by one detector------
31//AliCFTrackCutPid *pCut = new AliCFTrackCutPid("pion_sel","pion_sel");
32//Double_t priors[5]={...};
33//pCut->SetPartycleType(AliPID::kPion,kFALSE)
34//pCut->SetDetectors("DET"); ^^^^^^
35//pCut->SetPriors(priors);
36//
37// -----ID combining more detectors----
38//AliCFTrackCutPid *pCut = new AliCFTrackCutPid("pion_sel","pion_sel");
39//Double_t priors[5]={...};
40//pCut->SetPartycleType(AliPID::kPion,kTRUE)
41//pCut->SetDetectors("DET1 DET2 .."); ^^^^^
42//pCut->SetPriors(priors)
43//
44//Comments:
45//The user can choose not to identify a track if one residual beteween
46//the identified particle probability and one among all the other
47//probabilties is smaller than a predefined value.
48//The same can be done for the detector responses.
49//This happens by setting:
50//
51//pCut->SetMinDiffProb(0.005,kTRUE) //for probabilities
52//
53//pCut->SetMinDiffResp(0.005,kTRUE) //for det responses
54//
55//Annalisa.Mastroserio@ba.infn.it
56//
57
58
59#include "AliCFTrackCutPid.h"
60#include "AliLog.h"
61#include <TMath.h>
62#include <TList.h>
63
64ClassImp(AliCFTrackCutPid)
65
66//__________________________________________________________________________________
67// CUT ON TRACK PID
68//__________________________________________________________________________________
69AliCFTrackCutPid::AliCFTrackCutPid() :
70 AliCFCutBase(),
71 fCut(0.2),
72 fMinDiffResponse(0.0001),
73 fMinDiffProbability(0.001),
74 fgParticleType(10),
75 fgIsComb(kTRUE),
76 fCheckResponse(kFALSE),
77 fCheckSelection(kTRUE),
78 fIsPpriors(kFALSE),
79 fIsDetAND(kFALSE),
80 fXmin(-0.005),
81 fXmax(1.005),
82 fNbins(101),
83 fDetRestr(-1),
84 fiPartRestr(-1),
85 fDetProbRestr(1)
86{
87 //
88 //Default constructor
89 //
90 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
91 fPriors[j]=0.2;
92 fPriorsFunc[j]=0x0;
93 }
94
95 for(Int_t jDet=0; jDet< kNdets; jDet++) {
96 fDets[jDet]=kFALSE;
97 fDetsInAnd[jDet]=kFALSE;
98 }
99
100 InitialiseHisto();
101}
102//______________________________
103AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) :
104 AliCFCutBase(name,title),
105 fCut(0.2),
106 fMinDiffResponse(0.0001),
107 fMinDiffProbability(0.001),
108 fgParticleType(10),
109 fgIsComb(kTRUE),
110 fCheckResponse(kFALSE),
111 fCheckSelection(kTRUE),
112 fIsPpriors(kFALSE),
113 fIsDetAND(kFALSE),
114 fXmin(-0.005),
115 fXmax(1.005),
116 fNbins(101),
117 fDetRestr(-1),
118 fiPartRestr(-1),
119 fDetProbRestr(1)
120{
121 //
122 //Constructor
123 //
124 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
125 fPriors[j]=0.2;
126 fPriorsFunc[j]=0x0;
127 }
128
129 for(Int_t jDet=0; jDet< kNdets; jDet++) {
130 fDets[jDet]=kFALSE;
131 fDetsInAnd[jDet]=kFALSE;
132 }
133
134 InitialiseHisto();
135}
136//______________________________
137AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) :
138 AliCFCutBase(c),
139 fCut(c.fCut),
140 fMinDiffResponse(c.fMinDiffResponse),
141 fMinDiffProbability(c.fMinDiffProbability),
142 fgParticleType(c.fgParticleType),
143 fgIsComb(c.fgIsComb),
144 fCheckResponse(c.fCheckResponse),
145 fCheckSelection(c.fCheckSelection),
146 fIsPpriors(c.fIsPpriors),
147 fIsDetAND(c.fIsDetAND),
148 fXmin(c.fXmin),
149 fXmax(c.fXmax),
150 fNbins(c.fNbins),
151 fDetRestr(c.fDetRestr),
152 fiPartRestr(c.fiPartRestr),
153 fDetProbRestr(c.fDetProbRestr)
154{
155 //
156 //Copy constructor
157 //
158 for(Int_t i=0; i< kNdets ; i++ ) {
159 fDets[i]=c.fDets[i];
160 fDetsInAnd[i]=c.fDetsInAnd[i];
161 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
162 fhResp[i][iP]=c.fhResp[i][iP];
163 fhProb[i][iP]=c.fhProb[i][iP];
164 }
165 }
166 for(Int_t j=0; j< AliPID::kSPECIES; j++){
167 fPriors[j]=c.fPriors[j];
168 fPriorsFunc[j]=c.fPriorsFunc[j];
169 fhCombResp[j]=c.fhCombResp[j];
170 fhCombProb[j]=c.fhCombProb[j];
171 }
172}
173//______________________________
174AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c)
175{
176 //
177 // Assignment operator
178 //
179 if (this != &c) {
180 AliCFCutBase::operator=(c) ;
181 this->fCut=c.fCut;
182 this->fMinDiffResponse=c.fMinDiffResponse;
183 this->fMinDiffProbability=c.fMinDiffProbability;
184 this->fgParticleType=c.fgParticleType;
185 this->fgIsComb=c.fgIsComb;
186 this->fCheckResponse=c.fCheckResponse;
187 this->fCheckSelection=c.fCheckSelection;
188 this->fIsPpriors=c.fIsPpriors;
189 this->fIsDetAND=c.fIsDetAND;
190 this->fXmin=c.fXmin;
191 this->fXmax=c.fXmax;
192 this->fNbins=c.fNbins;
193 this->fDetRestr=c.fDetRestr;
194 this->fiPartRestr=c.fiPartRestr;
195 this->fDetProbRestr=c.fDetProbRestr;
196
197 for(Int_t i=0; i< kNdets ; i++ ) {
198 this->fDets[i]=c.fDets[i];
199 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
200 this->fhResp[i][iP]=c.fhResp[i][iP];
201 this->fhProb[i][iP]=c.fhProb[i][iP];
202 }
203 }
204
205 for(Int_t j=0; j< AliPID::kSPECIES; j++){
206 this->fPriors[j]=c.fPriors[j];
207 this->fhCombResp[j]=c.fhCombResp[j];
208 this->fhCombProb[j]=c.fhCombProb[j];
209 this-> fPriorsFunc[j]=c.fPriorsFunc[j];
210
211 }
212 }
213 return *this ;
214}
215//__________________________________________________________________________________
216AliCFTrackCutPid::~AliCFTrackCutPid() {
217 //
218 //dtor
219 //
220 for(Int_t i=0; i< kNdets ; i++ ) {
221 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
222 if(fhResp[i][iP])delete fhResp[i][iP];
223 if(fhProb[i][iP])delete fhProb[i][iP];
224 }
225 }
226
227 for(Int_t j=0; j< AliPID::kSPECIES; j++){
228 if(fhCombResp[j])delete fhCombResp[j];
229 if(fhCombProb[j])delete fhCombProb[j];
230
231 }
232
233}
234//__________________________________
235void AliCFTrackCutPid::SetDetectors(TString dets)
236{
237 //
238 // The string of detectors is translated into
239 // the respective booelan data members
240
241 if(dets.Contains("ITS")) {fDets[kITS]=kTRUE;}
242 if(dets.Contains("TPC")) {fDets[kTPC]=kTRUE;}
243 if(dets.Contains("TRD")) {fDets[kTRD]=kTRUE;}
244 if(dets.Contains("TOF")) {fDets[kTOF]=kTRUE;}
245 if(dets.Contains("HMPID")) {fDets[kHMPID]=kTRUE;}
246}
247//__________________________________
248void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES])
249{
250 //
251 // Sets the a priori concentrations
252 //
253 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriors[i]=r[i];
254}
255//__________________________________
256void AliCFTrackCutPid::SetPriorFunctions(TF1 *func[AliPID::kSPECIES])
257{
258 //
259 // Sets the momentu dependent a priori concentrations
260 //
261
262 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriorsFunc[i]=func[i];
263 fIsPpriors = kTRUE;
264}
265//____________________________________________
266void AliCFTrackCutPid::SetANDstatus(TString dets)
267{
268 //
269 //Sets the detectors to be in AND for the combined PID
270 //
271 if(dets.Contains("ITS") && fDets[kITS]) {fDetsInAnd[kITS]=kTRUE;}
272 if(dets.Contains("TPC") && fDets[kTPC]) {fDetsInAnd[kTPC]=kTRUE;}
273 if(dets.Contains("TRD") && fDets[kTRD]) {fDetsInAnd[kTRD]=kTRUE;}
274 if(dets.Contains("TOF") && fDets[kTOF]) {fDetsInAnd[kTOF]=kTRUE;}
275 if(dets.Contains("HMPID") && fDets[kHMPID]) {fDetsInAnd[kHMPID]=kTRUE;}
276
277 fIsDetAND = kTRUE;
278}
279//
280void AliCFTrackCutPid::SetDetectorProbabilityRestriction(TString det, Int_t iPart, Double_t upperprob)
281{
282 //
283 // Sets the detector, the particle and the probability
284 // limit.
285 //
286
287 if(det.Contains("ITS")) fDetRestr = kITS;
288 if(det.Contains("TPC")) fDetRestr = kTPC;
289 if(det.Contains("TRD")) fDetRestr = kTRD;
290 if(det.Contains("TOF")) fDetRestr = kTOF;
291 if(det.Contains("HMPID")) fDetRestr = kHMPID;
292 fiPartRestr = iPart;
293 fDetProbRestr = upperprob;
294}
295//__________________________________
296void AliCFTrackCutPid::TrackInfo(const AliESDtrack *pTrk, ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
297{
298 //
299 // Loads the responses of the XXX chosen detectors for the pTrk
300 // and the corresponding trk status. The final trk status is also loaded.
301 //
302 if(fDets[kITS]) {
303 pTrk->GetITSpid(pid[kITS]);
304 status[kITS]=AliESDtrack::kITSpid;
305 }
306 if(fDets[kTPC]) {
307 pTrk->GetTPCpid(pid[kTPC]);
308 status[kTPC]=AliESDtrack::kTPCpid;
309 }
310 if(fDets[kTRD]) {
311 pTrk->GetTRDpid(pid[kTRD]);
312 status[kTRD]=AliESDtrack::kTRDpid;
313 }
314 if(fDets[kTOF]) {
315 pTrk->GetTOFpid(pid[kTOF]);
316 status[kTOF]=AliESDtrack::kTOFpid;
317 }
318 if(fDets[kHMPID]) {
319 pTrk->GetHMPIDpid(pid[kHMPID]);
320 status[kHMPID]=AliESDtrack::kHMPIDpid;
321 }
322 if(fDetRestr>=0){
323 if(fDetRestr == kITS) pTrk->GetITSpid(pid[kITS]);
324 if(fDetRestr == kTPC) pTrk->GetITSpid(pid[kTPC]);
325 if(fDetRestr == kTRD) pTrk->GetITSpid(pid[kTRD]);
326 if(fDetRestr == kTOF) pTrk->GetITSpid(pid[kTOF]);
327 if(fDetRestr == kHMPID) pTrk->GetITSpid(pid[kHMPID]);
328 }
329
330 status[kNdets]=pTrk->GetStatus();
331 pTrk->GetESDpid(pid[kNdets]);
332}
333//__________________________________
334void AliCFTrackCutPid::SetPPriors(AliESDtrack *pTrk)
335{
336 //
337 //sets the mommentum dependent a priori concentrations
338 //
339
340 for(Int_t i=0; i< AliPID::kSPECIES; i++) {
341 if(pTrk->P()>fPriorsFunc[i]->GetXmin() && pTrk->P() < fPriorsFunc[i]->GetXmax()) fPriors[i]=fPriorsFunc[i]->Eval(pTrk->P());
342 else {AliInfo("the track momentum is not in the function range. Priors are equal") fPriors[i] = 0.2;}
343 }
344}
345//______________________________________
346ULong_t AliCFTrackCutPid::StatusForAND(ULong_t status[kNdets+1]) const
347{
348 //
349 //In case of AND of more detectors the AND-detector status combination.
350 //is calculated and also returned
351 //
352
353 ULong_t andstatus=0;
354 for(Int_t i=0; i< kNdets; i++) {
355 if(!fDetsInAnd[i]) continue;
356 andstatus = andstatus | status[i];
357 AliDebug(1,Form("trk status %lu %i AND-status combination: %lu",status[i],i,andstatus));
358 }
359 return andstatus;
360}
361//_______________________________________
362Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
363{
364 // Identifies the track if its probability is higher than the cut
365 // value. The default value is fCut=0.2, therefore the most probable
366 // one is identified by default. Here all the checks on how to identify
367 // the track are performed (single detector or combined PID,..., detector
368 // restriction probability
369 // Returns: integer corresponding to the identified particle
370
371 Int_t iPart=-1;
372
373 if(!fgIsComb){
374 Bool_t isDet=kFALSE;
375 for(Int_t i=0; i<kNdets; i++){
376 if(!fDets[i]) continue;
377 isDet=kTRUE;
378 AliDebug(1,Form("trk status %lu %i-det-pid status %lu -> combination: %lu",status[kNdets],i,status[i],status[kNdets]&status[i]));
379 if(!(status[kNdets]&status[i])){
380 iPart=-10;
381 AliDebug(1,Form("detector %i -> pid trk status not ok",i));
382 }
383 else {
384 AliDebug(1,Form("resp : %f %f %f %f %f",pid[i][0],pid[i][1],pid[i][2],pid[i][3],pid[i][4]));
385 if(fIsQAOn) iPart = IdentifyQA(pid[i],i);
386 else iPart = Identify(pid[i]);
387 }
388 }
389 if(!isDet){
390 AliDebug(1,Form(" !!! No detector selected, the ESD-pid response is considered"));
391 iPart = Identify(pid[kNdets]);
392 }
393 }else{
394 Double_t calcprob[5];
395 CombPID(status,pid,calcprob);
396 iPart = Identify(calcprob);
397 }
398
399
400 AliDebug(1,Form("selected particle: %i",iPart));
401
402 if(iPart >=0 && fiPartRestr>=0) {
403 AliPID restr(pid[fDetRestr]);
404 restr.SetPriors(fPriors);
405 AliDebug(1,Form("setted upper limit: %f det %i : probability %f ",fDetProbRestr,fDetRestr,restr.GetProbability((AliPID::EParticleType)fiPartRestr)));
406 if(restr.GetProbability((AliPID::EParticleType)fiPartRestr) > fDetProbRestr) {
407 iPart = kDetRestr;
408 AliDebug(1,"\n\n the detector restrictions refused the ID \n\n");
409 }
410 }//cross checks with one detector probability
411
412 AliDebug(1,Form("after the check the selected particle is %i",iPart));
413
414 return iPart;
415}
416//__________________________________
417Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const
418{
419 // Checks if there are no equal values and if the valus
420 // difference between the selected particle and the others i
421 // is higher than a lower limit.
422 // Returns: kTRUE= is acceptable
423
424 AliDebug(1,Form("input particle: %i",iPsel));
425 Bool_t ck=kTRUE;
426
427 if(iPsel<0) ck=kFALSE;
428
429 else {
430 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
431 if(j!=iPsel && TMath::Abs(p[j]-p[iPsel])<minDiff) ck=kFALSE;
432 }
433 if(!ck) AliDebug(1,"the values are too close ");
434 }
435 return ck;
436}
437//___________________________________________
438Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const
439{
440 //
441 // The identification is actually performed here with possible
442 // checks on the det responses and/or probabilities
443 //
444 Int_t iPart = -1;
445
446 AliPID getpid(pid,kTRUE);
447
448 if(fgIsComb) {
449 Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
450 getpid.SetPriors(priors);
451 }
452 else getpid.SetPriors(fPriors);
453
454
455 AliPID::EParticleType sel = getpid.GetMostProbable();
456 Double_t probability[AliPID::kSPECIES];
457 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
458
459
460 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
461 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
462 AliDebug(1,Form("probab : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
463 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
464 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
465 return iPart;
466}
467//___________________________________________
468Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t idets) const
469{
470 //
471 // The same as Identify, but with the QA histo filling
472 //
473 //
474
475 Int_t iPart = -1;
476 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
477
478 AliPID getpid(pid,kTRUE);
479
480 if(fgIsComb) {
481 Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
482 getpid.SetPriors(priors);
483 }
484 else getpid.SetPriors(fPriors);
485
486 AliPID::EParticleType sel = getpid.GetMostProbable();
487 Double_t probability[AliPID::kSPECIES];
488 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) {
489 probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
490 fhProb[idets][iP]->Fill(probability[iP]);
491 }
492
493 AliPID toresp(pid,kTRUE); Double_t qapriors[5]={0.2,0.2,0.2,0.2,0.2};
494 toresp.SetPriors(qapriors);
495 for(Int_t iPr=0; iPr<AliPID::kSPECIES; iPr++) fhResp[idets][iPr]->Fill(toresp.GetProbability((AliPID::EParticleType)iPr));
496
497 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
498 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
499 AliDebug(1,Form("probab : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
500 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
501 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
502 return iPart;
503}
504//___________________________________________
505Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
506 //
507 // method for the pid-cut selction
508 //
509
510 if (!track) return kFALSE ;
511 TString className(track->ClassName());
512 if (className.CompareTo("AliESDtrack") != 0) {
513 AliError("obj must point to a AliESDtrack ");
514 return kFALSE ;
515 }
516
517 AliESDtrack *esdTrack = (AliESDtrack *)track;
518 ULong_t status[kNdets+1]={0,0,0,0,0,0};
519 Double_t pid[kNdets+1][AliPID::kSPECIES];
520 TrackInfo(esdTrack,status,pid);
521 if(fIsPpriors) SetPPriors(esdTrack);
522 if(GetID(status,pid)==fgParticleType) return kTRUE;
523 else return kFALSE;
524}
525//__________________________________
526void AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const
527{
528 //
529 // Calculates the combined PID according to the chosen detectors.
530 // and provides the array of probabilities
531 //
532
533 Bool_t isdet=kFALSE;
534 Double_t sum=0.;
535 Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.};
536 Double_t comb[AliPID::kSPECIES]={0.,0.,0.,0.,0.};
537 Double_t priors[AliPID::kSPECIES]={0.2,0.2,0.2,0.2,0.2};
538 ULong_t andstatus =0;
539 if(fIsDetAND) {
540 andstatus = StatusForAND(status);
541 AliDebug(1,Form("AND combination %lu",andstatus));
542 }
543 for(Int_t j=0; j<AliPID::kSPECIES; j++){
544 for(Int_t i=0; i< kNdets; i++){
545 if(!fDets[i]) continue;
546 if(status[kNdets]&status[i]) {
547 if(fIsDetAND) {
548 ULong_t checkstatus = status[kNdets]&andstatus;
549 if(checkstatus != andstatus) continue;
550 else {
551 prod[j]*=pid[i][j];
552 isdet = kTRUE;
553 AliDebug(1,Form("-----> trk status %lu and status %lu -> trk-ANDdetector status combination %lu",status[kNdets],andstatus,status[kNdets]&andstatus));
554 AliDebug(1,Form("In det loop %i -> particle %i response is %f",i,j,pid[i][j]));
555 }
556 }
557 else {
558 prod[j]*=pid[i][j];
559 isdet=kTRUE;
560 AliDebug(1,Form("In det loop %i -> particle %i response is %f",i,j,pid[i][j]));
561 }
562 }//combined mode
563 }//loop on dets
564 }//loop on species
565
566 if(fIsQAOn) {
567 for(Int_t iqa =0; iqa < kNdets; iqa++){
568 if(!fDets[iqa]) continue;
569 AliPID normresp(pid[iqa]);
570 normresp.SetPriors(priors);
571 for(Int_t ip =0; ip< AliPID::kSPECIES; ip++){
572 if(!fhResp[iqa][ip]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",iqa,ip));}
573 else fhResp[iqa][ip]->Fill(normresp.GetProbability((AliPID::EParticleType)ip));
574 }//loop on part
575 }//loop on dets
576 }//if qa
577
578 if(!isdet) {
579 AliDebug(1,"No proper status for the combined pid -> probabilities are set to zero");
580 for(Int_t nn=0; nn<AliPID::kSPECIES; nn++) {combpid[nn]=0;}
581 }
582
583 else{
584 for(Int_t k=0; k<AliPID::kSPECIES; k++){
585 if(fIsQAOn) {
586 if(!fhCombResp[k]) AliDebug(1,Form("no fhCombResp[%i], check if pidcut->Init() was called",k));
587 else fhCombResp[k]->Fill(prod[k]);
588 }
589 AliDebug(1,Form("species %i priors %f and prod %f",k,fPriors[k],prod[k]));
590 comb[k]=fPriors[k]*prod[k];
591 AliDebug(1,Form("comb %i %f",k,comb[k]));
592 sum+=comb[k];
593 }
594
595 if(sum == 0) {
596 AliDebug(1,"Check the detector responses or the priors, their combined products are zero");
597 return;
598 }
599 for(Int_t n=0; n<AliPID::kSPECIES; n++) {
600 combpid[n]=fPriors[n]*prod[n]/sum;
601 if(fIsQAOn) {
602 if(!fhCombProb[n]) Printf("no fhCombResp[%i], check if pidcut->Init() was called",n);
603 fhCombProb[n]->Fill(combpid[n]);
604 }
605 }
606 }
607}
608//__________________________________________
609//
610//QA part
611//_________________________________________
612void AliCFTrackCutPid::InitialiseHisto()
613{
614 //
615 //QA histo initialization
616 //
617 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++){
618 fhCombResp[iP]=0x0;
619 fhCombProb[iP]=0x0;
620 for(Int_t iDet =0; iDet<kNdets; iDet++){
621 fhResp[iDet][iP]=0x0;
622 fhProb[iDet][iP]=0x0;
623 }
624 }
625}
626//______________________________________________
627void AliCFTrackCutPid::Init()
628{
629 //
630 // initialises QA histograms
631 //
632
633 if(fIsQAOn) DefineHistograms();
634}
635//_________________________________________________
636void AliCFTrackCutPid::DefineHistograms()
637{
638 //
639 //QA histo booking
640 //
641 char *detect[5]={"ITS","TPC","TRD","TOF","HMPID"};
642 char *partic[5]={"electron","muon","pion","kaon","proton"};
643
644 for(Int_t iDet =0; iDet< kNdets; iDet++)
645 {
646 if(!fDets[iDet]) continue;
647 for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){
648 fhResp[iDet][iP] = new TH1F(Form("rDet%iPart%i",iDet,iP),Form("%s %s response ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
649 fhProb[iDet][iP] = new TH1F(Form("pDet%iPart%i",iDet,iP),Form("%s %s probability ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
650 }
651 }
652
653 if(fgIsComb){
654 for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
655 {
656 fhCombResp[iPart] = new TH1F(Form("rCombPart%i",iPart),Form(" %s combined response ",partic[iPart]),fNbins,fXmin,fXmax);
657 fhCombProb[iPart] = new TH1F(Form("pCombPart%i",iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax);
658 }
659 }
660}
661//___________________________________________________
662
663void AliCFTrackCutPid::AddQAHistograms(TList *qalist) const
664{
665 //
666 // adds QA histograms in a TList
667 //
668 if(!fIsQAOn) return;
669
670 if(fgIsComb){
671 for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){
672 qalist->Add(fhCombResp[iPart]);
673 qalist->Add(fhCombProb[iPart]);
674 }
675 }
676
677 for(Int_t iDet=0; iDet<kNdets; iDet++){
678 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
679 if(!fgIsComb)qalist->Add(fhResp[iDet][iP]);
680 if(!fgIsComb)qalist->Add(fhProb[iDet][iP]);
681 }
682 }
683}