]> 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
Fixed Coding Conventions for class AliACORDEv1 (Mario RC)
[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 //
7f0c5c1e 220
563113d0 221 for(Int_t i=0; i< kNdets ; i++ ) {
222 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
223 if(fhResp[i][iP])delete fhResp[i][iP];
224 if(fhProb[i][iP])delete fhProb[i][iP];
225 }
226 }
227
228 for(Int_t j=0; j< AliPID::kSPECIES; j++){
229 if(fhCombResp[j])delete fhCombResp[j];
230 if(fhCombProb[j])delete fhCombProb[j];
231
232 }
563113d0 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;}
7f0c5c1e 246
247 if(dets.Contains("ALL")) for(Int_t i=0; i< kNdets ; i++) fDets[i]=kTRUE;
563113d0 248}
249//__________________________________
250void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES])
251{
252 //
253 // Sets the a priori concentrations
254 //
255 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriors[i]=r[i];
256}
257//__________________________________
258void AliCFTrackCutPid::SetPriorFunctions(TF1 *func[AliPID::kSPECIES])
259{
260 //
261 // Sets the momentu dependent a priori concentrations
262 //
263
264 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriorsFunc[i]=func[i];
265 fIsPpriors = kTRUE;
266}
267//____________________________________________
268void AliCFTrackCutPid::SetANDstatus(TString dets)
269{
270 //
271 //Sets the detectors to be in AND for the combined PID
272 //
273 if(dets.Contains("ITS") && fDets[kITS]) {fDetsInAnd[kITS]=kTRUE;}
274 if(dets.Contains("TPC") && fDets[kTPC]) {fDetsInAnd[kTPC]=kTRUE;}
275 if(dets.Contains("TRD") && fDets[kTRD]) {fDetsInAnd[kTRD]=kTRUE;}
276 if(dets.Contains("TOF") && fDets[kTOF]) {fDetsInAnd[kTOF]=kTRUE;}
277 if(dets.Contains("HMPID") && fDets[kHMPID]) {fDetsInAnd[kHMPID]=kTRUE;}
278
279 fIsDetAND = kTRUE;
280}
281//
282void AliCFTrackCutPid::SetDetectorProbabilityRestriction(TString det, Int_t iPart, Double_t upperprob)
283{
284 //
285 // Sets the detector, the particle and the probability
286 // limit.
287 //
288
289 if(det.Contains("ITS")) fDetRestr = kITS;
290 if(det.Contains("TPC")) fDetRestr = kTPC;
291 if(det.Contains("TRD")) fDetRestr = kTRD;
292 if(det.Contains("TOF")) fDetRestr = kTOF;
293 if(det.Contains("HMPID")) fDetRestr = kHMPID;
294 fiPartRestr = iPart;
295 fDetProbRestr = upperprob;
296}
297//__________________________________
298void AliCFTrackCutPid::TrackInfo(const AliESDtrack *pTrk, ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
299{
300 //
301 // Loads the responses of the XXX chosen detectors for the pTrk
302 // and the corresponding trk status. The final trk status is also loaded.
303 //
304 if(fDets[kITS]) {
305 pTrk->GetITSpid(pid[kITS]);
306 status[kITS]=AliESDtrack::kITSpid;
307 }
308 if(fDets[kTPC]) {
309 pTrk->GetTPCpid(pid[kTPC]);
310 status[kTPC]=AliESDtrack::kTPCpid;
311 }
312 if(fDets[kTRD]) {
313 pTrk->GetTRDpid(pid[kTRD]);
314 status[kTRD]=AliESDtrack::kTRDpid;
315 }
316 if(fDets[kTOF]) {
317 pTrk->GetTOFpid(pid[kTOF]);
318 status[kTOF]=AliESDtrack::kTOFpid;
319 }
320 if(fDets[kHMPID]) {
321 pTrk->GetHMPIDpid(pid[kHMPID]);
322 status[kHMPID]=AliESDtrack::kHMPIDpid;
323 }
324 if(fDetRestr>=0){
325 if(fDetRestr == kITS) pTrk->GetITSpid(pid[kITS]);
326 if(fDetRestr == kTPC) pTrk->GetITSpid(pid[kTPC]);
327 if(fDetRestr == kTRD) pTrk->GetITSpid(pid[kTRD]);
328 if(fDetRestr == kTOF) pTrk->GetITSpid(pid[kTOF]);
329 if(fDetRestr == kHMPID) pTrk->GetITSpid(pid[kHMPID]);
330 }
331
332 status[kNdets]=pTrk->GetStatus();
333 pTrk->GetESDpid(pid[kNdets]);
334}
335//__________________________________
336void AliCFTrackCutPid::SetPPriors(AliESDtrack *pTrk)
337{
338 //
339 //sets the mommentum dependent a priori concentrations
340 //
341
342 for(Int_t i=0; i< AliPID::kSPECIES; i++) {
343 if(pTrk->P()>fPriorsFunc[i]->GetXmin() && pTrk->P() < fPriorsFunc[i]->GetXmax()) fPriors[i]=fPriorsFunc[i]->Eval(pTrk->P());
344 else {AliInfo("the track momentum is not in the function range. Priors are equal") fPriors[i] = 0.2;}
345 }
346}
347//______________________________________
348ULong_t AliCFTrackCutPid::StatusForAND(ULong_t status[kNdets+1]) const
349{
350 //
351 //In case of AND of more detectors the AND-detector status combination.
352 //is calculated and also returned
353 //
354
355 ULong_t andstatus=0;
356 for(Int_t i=0; i< kNdets; i++) {
357 if(!fDetsInAnd[i]) continue;
358 andstatus = andstatus | status[i];
359 AliDebug(1,Form("trk status %lu %i AND-status combination: %lu",status[i],i,andstatus));
360 }
361 return andstatus;
362}
363//_______________________________________
364Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
365{
366 // Identifies the track if its probability is higher than the cut
367 // value. The default value is fCut=0.2, therefore the most probable
368 // one is identified by default. Here all the checks on how to identify
369 // the track are performed (single detector or combined PID,..., detector
370 // restriction probability
371 // Returns: integer corresponding to the identified particle
372
373 Int_t iPart=-1;
374
375 if(!fgIsComb){
376 Bool_t isDet=kFALSE;
377 for(Int_t i=0; i<kNdets; i++){
378 if(!fDets[i]) continue;
379 isDet=kTRUE;
380 AliDebug(1,Form("trk status %lu %i-det-pid status %lu -> combination: %lu",status[kNdets],i,status[i],status[kNdets]&status[i]));
381 if(!(status[kNdets]&status[i])){
382 iPart=-10;
383 AliDebug(1,Form("detector %i -> pid trk status not ok",i));
384 }
385 else {
386 AliDebug(1,Form("resp : %f %f %f %f %f",pid[i][0],pid[i][1],pid[i][2],pid[i][3],pid[i][4]));
387 if(fIsQAOn) iPart = IdentifyQA(pid[i],i);
388 else iPart = Identify(pid[i]);
389 }
390 }
391 if(!isDet){
392 AliDebug(1,Form(" !!! No detector selected, the ESD-pid response is considered"));
393 iPart = Identify(pid[kNdets]);
394 }
395 }else{
396 Double_t calcprob[5];
397 CombPID(status,pid,calcprob);
398 iPart = Identify(calcprob);
399 }
400
401
402 AliDebug(1,Form("selected particle: %i",iPart));
403
404 if(iPart >=0 && fiPartRestr>=0) {
405 AliPID restr(pid[fDetRestr]);
406 restr.SetPriors(fPriors);
407 AliDebug(1,Form("setted upper limit: %f det %i : probability %f ",fDetProbRestr,fDetRestr,restr.GetProbability((AliPID::EParticleType)fiPartRestr)));
408 if(restr.GetProbability((AliPID::EParticleType)fiPartRestr) > fDetProbRestr) {
409 iPart = kDetRestr;
410 AliDebug(1,"\n\n the detector restrictions refused the ID \n\n");
411 }
412 }//cross checks with one detector probability
413
414 AliDebug(1,Form("after the check the selected particle is %i",iPart));
415
416 return iPart;
417}
418//__________________________________
419Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const
420{
421 // Checks if there are no equal values and if the valus
422 // difference between the selected particle and the others i
423 // is higher than a lower limit.
424 // Returns: kTRUE= is acceptable
425
7f0c5c1e 426 AliDebug(2,Form("input particle: %i",iPsel));
563113d0 427 Bool_t ck=kTRUE;
428
429 if(iPsel<0) ck=kFALSE;
430
431 else {
432 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
433 if(j!=iPsel && TMath::Abs(p[j]-p[iPsel])<minDiff) ck=kFALSE;
434 }
435 if(!ck) AliDebug(1,"the values are too close ");
436 }
437 return ck;
438}
439//___________________________________________
440Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const
441{
442 //
443 // The identification is actually performed here with possible
444 // checks on the det responses and/or probabilities
445 //
446 Int_t iPart = -1;
447
448 AliPID getpid(pid,kTRUE);
449
450 if(fgIsComb) {
451 Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
452 getpid.SetPriors(priors);
453 }
454 else getpid.SetPriors(fPriors);
455
456
457 AliPID::EParticleType sel = getpid.GetMostProbable();
458 Double_t probability[AliPID::kSPECIES];
459 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
460
461
462 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
7f0c5c1e 463 AliDebug(2,Form("calc response : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
464 AliDebug(2,Form("probabilities : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
465
563113d0 466 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
467 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
7f0c5c1e 468
563113d0 469 return iPart;
470}
471//___________________________________________
472Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t idets) const
473{
474 //
475 // The same as Identify, but with the QA histo filling
476 //
477 //
478
479 Int_t iPart = -1;
480 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
481
482 AliPID getpid(pid,kTRUE);
483
484 if(fgIsComb) {
485 Double_t priors[5]={0.2,0.2,0.2,0.2,0.2};
486 getpid.SetPriors(priors);
487 }
488 else getpid.SetPriors(fPriors);
489
490 AliPID::EParticleType sel = getpid.GetMostProbable();
491 Double_t probability[AliPID::kSPECIES];
492 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) {
493 probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
494 fhProb[idets][iP]->Fill(probability[iP]);
495 }
496
497 AliPID toresp(pid,kTRUE); Double_t qapriors[5]={0.2,0.2,0.2,0.2,0.2};
498 toresp.SetPriors(qapriors);
499 for(Int_t iPr=0; iPr<AliPID::kSPECIES; iPr++) fhResp[idets][iPr]->Fill(toresp.GetProbability((AliPID::EParticleType)iPr));
500
501 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
502 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
503 AliDebug(1,Form("probab : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
504 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
505 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
506 return iPart;
507}
508//___________________________________________
509Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
510 //
511 // method for the pid-cut selction
512 //
513
514 if (!track) return kFALSE ;
515 TString className(track->ClassName());
516 if (className.CompareTo("AliESDtrack") != 0) {
517 AliError("obj must point to a AliESDtrack ");
518 return kFALSE ;
519 }
520
521 AliESDtrack *esdTrack = (AliESDtrack *)track;
522 ULong_t status[kNdets+1]={0,0,0,0,0,0};
523 Double_t pid[kNdets+1][AliPID::kSPECIES];
524 TrackInfo(esdTrack,status,pid);
525 if(fIsPpriors) SetPPriors(esdTrack);
526 if(GetID(status,pid)==fgParticleType) return kTRUE;
527 else return kFALSE;
528}
529//__________________________________
530void AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const
531{
532 //
533 // Calculates the combined PID according to the chosen detectors.
534 // and provides the array of probabilities
535 //
536
537 Bool_t isdet=kFALSE;
538 Double_t sum=0.;
539 Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.};
540 Double_t comb[AliPID::kSPECIES]={0.,0.,0.,0.,0.};
541 Double_t priors[AliPID::kSPECIES]={0.2,0.2,0.2,0.2,0.2};
7f0c5c1e 542
563113d0 543 ULong_t andstatus =0;
544 if(fIsDetAND) {
545 andstatus = StatusForAND(status);
546 AliDebug(1,Form("AND combination %lu",andstatus));
547 }
548 for(Int_t j=0; j<AliPID::kSPECIES; j++){
549 for(Int_t i=0; i< kNdets; i++){
550 if(!fDets[i]) continue;
551 if(status[kNdets]&status[i]) {
552 if(fIsDetAND) {
553 ULong_t checkstatus = status[kNdets]&andstatus;
554 if(checkstatus != andstatus) continue;
555 else {
556 prod[j]*=pid[i][j];
557 isdet = kTRUE;
558 AliDebug(1,Form("-----> trk status %lu and status %lu -> trk-ANDdetector status combination %lu",status[kNdets],andstatus,status[kNdets]&andstatus));
7f0c5c1e 559 AliDebug(1,Form("In det %i -> particle %i response is %f",i,j,pid[i][j]));
563113d0 560 }
561 }
562 else {
563 prod[j]*=pid[i][j];
564 isdet=kTRUE;
7f0c5c1e 565 AliDebug(2,Form("In det %i -> particle %i response is %f",i,j,pid[i][j]));
563113d0 566 }
567 }//combined mode
568 }//loop on dets
569 }//loop on species
570
571 if(fIsQAOn) {
572 for(Int_t iqa =0; iqa < kNdets; iqa++){
573 if(!fDets[iqa]) continue;
574 AliPID normresp(pid[iqa]);
575 normresp.SetPriors(priors);
576 for(Int_t ip =0; ip< AliPID::kSPECIES; ip++){
577 if(!fhResp[iqa][ip]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",iqa,ip));}
578 else fhResp[iqa][ip]->Fill(normresp.GetProbability((AliPID::EParticleType)ip));
579 }//loop on part
580 }//loop on dets
581 }//if qa
582
7f0c5c1e 583
584 Double_t add = 0; for(Int_t isumm=0; isumm<5; isumm++) add+=prod[isumm];
585 if(add>0) AliDebug(1,Form("calculated comb response: %f %f %f %f %f",prod[0]/add,prod[1]/add,prod[2]/add,prod[3]/add,prod[4]/add));
586 AliDebug(1,Form("the ESDpid response: %f %f %f %f %f",pid[kNdets][0],pid[kNdets][1],pid[kNdets][2],pid[kNdets][3],pid[kNdets][4]));
587
563113d0 588 if(!isdet) {
7f0c5c1e 589 AliWarning("\n !! No detector found for the combined pid --> responses are from the ESDpid !! \n");
590 Double_t sumesdpid=0;
591 for(Int_t nn=0; nn<AliPID::kSPECIES; nn++) {
592 combpid[nn]=pid[kNdets][nn];
593 sumesdpid+=fPriors[nn]*combpid[nn];
594 if(fIsQAOn) {
595 if(!fhCombResp[nn]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",nn));
596 else fhCombResp[nn]->Fill(combpid[nn]);
597 }//fIsQAOn
598 }//loop on species
599 if(sumesdpid > 0) for(Int_t ih = 0; ih < AliPID::kSPECIES; ih++) fhCombProb[ih]->Fill(fPriors[ih]*combpid[ih]/sumesdpid);
600 else AliDebug(1,"priors or ESDpid are zero, please check them");
601 }// end no det
563113d0 602
603 else{
604 for(Int_t k=0; k<AliPID::kSPECIES; k++){
605 if(fIsQAOn) {
7f0c5c1e 606 if(!fhCombResp[k]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called \n",k));
563113d0 607 else fhCombResp[k]->Fill(prod[k]);
608 }
7f0c5c1e 609 AliDebug(1,Form("species: %i priors %f and combined prod %f",k,fPriors[k],prod[k]));
563113d0 610 comb[k]=fPriors[k]*prod[k];
7f0c5c1e 611 AliDebug(1,Form("combined probability %i %f",k,comb[k]));
563113d0 612 sum+=comb[k];
613 }
614
615 if(sum == 0) {
616 AliDebug(1,"Check the detector responses or the priors, their combined products are zero");
617 return;
618 }
619 for(Int_t n=0; n<AliPID::kSPECIES; n++) {
620 combpid[n]=fPriors[n]*prod[n]/sum;
621 if(fIsQAOn) {
7f0c5c1e 622 if(!fhCombProb[n]) AliDebug(1,Form("no fhCombRespi[%i] defined, check if pidcut->Init() was called",n));
623 else fhCombProb[n]->Fill(combpid[n]);
563113d0 624 }
625 }
7f0c5c1e 626 }//end else
563113d0 627}
628//__________________________________________
629//
630//QA part
631//_________________________________________
632void AliCFTrackCutPid::InitialiseHisto()
633{
634 //
635 //QA histo initialization
636 //
637 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++){
638 fhCombResp[iP]=0x0;
639 fhCombProb[iP]=0x0;
640 for(Int_t iDet =0; iDet<kNdets; iDet++){
641 fhResp[iDet][iP]=0x0;
642 fhProb[iDet][iP]=0x0;
643 }
644 }
645}
646//______________________________________________
647void AliCFTrackCutPid::Init()
648{
649 //
650 // initialises QA histograms
651 //
652
653 if(fIsQAOn) DefineHistograms();
654}
655//_________________________________________________
656void AliCFTrackCutPid::DefineHistograms()
657{
658 //
659 //QA histo booking
660 //
7f0c5c1e 661 char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"};
662 char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"};
563113d0 663
7f0c5c1e 664 for(Int_t iDet =0; iDet< kNdets; iDet++)
563113d0 665 {
666 if(!fDets[iDet]) continue;
667 for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){
668 fhResp[iDet][iP] = new TH1F(Form("rDet%iPart%i",iDet,iP),Form("%s %s response ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
669 fhProb[iDet][iP] = new TH1F(Form("pDet%iPart%i",iDet,iP),Form("%s %s probability ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
670 }
671 }
672
7f0c5c1e 673
674if(fgIsComb)
675 {
563113d0 676 for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
677 {
678 fhCombResp[iPart] = new TH1F(Form("rCombPart%i",iPart),Form(" %s combined response ",partic[iPart]),fNbins,fXmin,fXmax);
679 fhCombProb[iPart] = new TH1F(Form("pCombPart%i",iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax);
680 }
681 }
682}
683//___________________________________________________
684
685void AliCFTrackCutPid::AddQAHistograms(TList *qalist) const
686{
687 //
688 // adds QA histograms in a TList
689 //
690 if(!fIsQAOn) return;
691
692 if(fgIsComb){
693 for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){
694 qalist->Add(fhCombResp[iPart]);
695 qalist->Add(fhCombProb[iPart]);
696 }
697 }
698
699 for(Int_t iDet=0; iDet<kNdets; iDet++){
7f0c5c1e 700 if(!fDets[iDet]) continue;
563113d0 701 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
7f0c5c1e 702 qalist->Add(fhResp[iDet][iP]);
703 qalist->Add(fhProb[iDet][iP]);
563113d0 704 }
705 }
7f0c5c1e 706
563113d0 707}