]> 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
added function AliCFContainer::MakeSlice reading user-defined steps.
[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
f9658e11 59
563113d0 60#include "AliCFTrackCutPid.h"
61#include "AliLog.h"
62#include <TMath.h>
63#include <TList.h>
64
65ClassImp(AliCFTrackCutPid)
66
67//__________________________________________________________________________________
68// CUT ON TRACK PID
69//__________________________________________________________________________________
70AliCFTrackCutPid::AliCFTrackCutPid() :
71 AliCFCutBase(),
72 fCut(0.2),
73 fMinDiffResponse(0.0001),
74 fMinDiffProbability(0.001),
75 fgParticleType(10),
76 fgIsComb(kTRUE),
bc6176fb 77 fgIsAOD(kFALSE),
563113d0 78 fCheckResponse(kFALSE),
79 fCheckSelection(kTRUE),
80 fIsPpriors(kFALSE),
81 fIsDetAND(kFALSE),
82 fXmin(-0.005),
83 fXmax(1.005),
84 fNbins(101),
85 fDetRestr(-1),
86 fiPartRestr(-1),
87 fDetProbRestr(1)
88{
89 //
90 //Default constructor
91 //
92 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
93 fPriors[j]=0.2;
94 fPriorsFunc[j]=0x0;
95 }
96
97 for(Int_t jDet=0; jDet< kNdets; jDet++) {
98 fDets[jDet]=kFALSE;
99 fDetsInAnd[jDet]=kFALSE;
100 }
101
102 InitialiseHisto();
103}
104//______________________________
105AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) :
106 AliCFCutBase(name,title),
107 fCut(0.2),
108 fMinDiffResponse(0.0001),
109 fMinDiffProbability(0.001),
110 fgParticleType(10),
111 fgIsComb(kTRUE),
bc6176fb 112 fgIsAOD(kFALSE),
563113d0 113 fCheckResponse(kFALSE),
114 fCheckSelection(kTRUE),
115 fIsPpriors(kFALSE),
116 fIsDetAND(kFALSE),
117 fXmin(-0.005),
118 fXmax(1.005),
119 fNbins(101),
120 fDetRestr(-1),
121 fiPartRestr(-1),
122 fDetProbRestr(1)
123{
124 //
125 //Constructor
126 //
127 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
128 fPriors[j]=0.2;
129 fPriorsFunc[j]=0x0;
130 }
131
132 for(Int_t jDet=0; jDet< kNdets; jDet++) {
133 fDets[jDet]=kFALSE;
134 fDetsInAnd[jDet]=kFALSE;
135 }
136
137 InitialiseHisto();
138}
139//______________________________
140AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) :
141 AliCFCutBase(c),
142 fCut(c.fCut),
143 fMinDiffResponse(c.fMinDiffResponse),
144 fMinDiffProbability(c.fMinDiffProbability),
145 fgParticleType(c.fgParticleType),
146 fgIsComb(c.fgIsComb),
bc6176fb 147 fgIsAOD(c.fgIsAOD),
563113d0 148 fCheckResponse(c.fCheckResponse),
149 fCheckSelection(c.fCheckSelection),
150 fIsPpriors(c.fIsPpriors),
151 fIsDetAND(c.fIsDetAND),
152 fXmin(c.fXmin),
153 fXmax(c.fXmax),
154 fNbins(c.fNbins),
155 fDetRestr(c.fDetRestr),
156 fiPartRestr(c.fiPartRestr),
157 fDetProbRestr(c.fDetProbRestr)
158{
159 //
160 //Copy constructor
161 //
162 for(Int_t i=0; i< kNdets ; i++ ) {
163 fDets[i]=c.fDets[i];
164 fDetsInAnd[i]=c.fDetsInAnd[i];
165 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
166 fhResp[i][iP]=c.fhResp[i][iP];
167 fhProb[i][iP]=c.fhProb[i][iP];
168 }
169 }
170 for(Int_t j=0; j< AliPID::kSPECIES; j++){
171 fPriors[j]=c.fPriors[j];
172 fPriorsFunc[j]=c.fPriorsFunc[j];
173 fhCombResp[j]=c.fhCombResp[j];
174 fhCombProb[j]=c.fhCombProb[j];
175 }
176}
177//______________________________
178AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c)
179{
180 //
181 // Assignment operator
182 //
183 if (this != &c) {
184 AliCFCutBase::operator=(c) ;
185 this->fCut=c.fCut;
186 this->fMinDiffResponse=c.fMinDiffResponse;
bc6176fb 187 this->fMinDiffProbability=c.fMinDiffProbability;
563113d0 188 this->fgParticleType=c.fgParticleType;
189 this->fgIsComb=c.fgIsComb;
bc6176fb 190 this->fgIsAOD=c.fgIsAOD;
563113d0 191 this->fCheckResponse=c.fCheckResponse;
192 this->fCheckSelection=c.fCheckSelection;
193 this->fIsPpriors=c.fIsPpriors;
194 this->fIsDetAND=c.fIsDetAND;
195 this->fXmin=c.fXmin;
196 this->fXmax=c.fXmax;
197 this->fNbins=c.fNbins;
198 this->fDetRestr=c.fDetRestr;
199 this->fiPartRestr=c.fiPartRestr;
200 this->fDetProbRestr=c.fDetProbRestr;
201
202 for(Int_t i=0; i< kNdets ; i++ ) {
203 this->fDets[i]=c.fDets[i];
204 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
205 this->fhResp[i][iP]=c.fhResp[i][iP];
206 this->fhProb[i][iP]=c.fhProb[i][iP];
207 }
208 }
209
210 for(Int_t j=0; j< AliPID::kSPECIES; j++){
211 this->fPriors[j]=c.fPriors[j];
212 this->fhCombResp[j]=c.fhCombResp[j];
213 this->fhCombProb[j]=c.fhCombProb[j];
214 this-> fPriorsFunc[j]=c.fPriorsFunc[j];
215
216 }
217 }
218 return *this ;
219}
220//__________________________________________________________________________________
221AliCFTrackCutPid::~AliCFTrackCutPid() {
222 //
223 //dtor
224 //
7f0c5c1e 225
563113d0 226 for(Int_t i=0; i< kNdets ; i++ ) {
227 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
228 if(fhResp[i][iP])delete fhResp[i][iP];
229 if(fhProb[i][iP])delete fhProb[i][iP];
230 }
231 }
232
233 for(Int_t j=0; j< AliPID::kSPECIES; j++){
234 if(fhCombResp[j])delete fhCombResp[j];
235 if(fhCombProb[j])delete fhCombProb[j];
236
237 }
563113d0 238}
239//__________________________________
240void AliCFTrackCutPid::SetDetectors(TString dets)
241{
242 //
243 // The string of detectors is translated into
244 // the respective booelan data members
245
246 if(dets.Contains("ITS")) {fDets[kITS]=kTRUE;}
247 if(dets.Contains("TPC")) {fDets[kTPC]=kTRUE;}
248 if(dets.Contains("TRD")) {fDets[kTRD]=kTRUE;}
249 if(dets.Contains("TOF")) {fDets[kTOF]=kTRUE;}
250 if(dets.Contains("HMPID")) {fDets[kHMPID]=kTRUE;}
7f0c5c1e 251
252 if(dets.Contains("ALL")) for(Int_t i=0; i< kNdets ; i++) fDets[i]=kTRUE;
563113d0 253}
254//__________________________________
255void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES])
256{
257 //
258 // Sets the a priori concentrations
259 //
260 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriors[i]=r[i];
261}
262//__________________________________
263void AliCFTrackCutPid::SetPriorFunctions(TF1 *func[AliPID::kSPECIES])
264{
265 //
266 // Sets the momentu dependent a priori concentrations
267 //
268
269 for(Int_t i=0; i<AliPID::kSPECIES; i++) fPriorsFunc[i]=func[i];
270 fIsPpriors = kTRUE;
271}
272//____________________________________________
273void AliCFTrackCutPid::SetANDstatus(TString dets)
274{
275 //
276 //Sets the detectors to be in AND for the combined PID
277 //
278 if(dets.Contains("ITS") && fDets[kITS]) {fDetsInAnd[kITS]=kTRUE;}
279 if(dets.Contains("TPC") && fDets[kTPC]) {fDetsInAnd[kTPC]=kTRUE;}
280 if(dets.Contains("TRD") && fDets[kTRD]) {fDetsInAnd[kTRD]=kTRUE;}
281 if(dets.Contains("TOF") && fDets[kTOF]) {fDetsInAnd[kTOF]=kTRUE;}
282 if(dets.Contains("HMPID") && fDets[kHMPID]) {fDetsInAnd[kHMPID]=kTRUE;}
283
284 fIsDetAND = kTRUE;
285}
286//
287void AliCFTrackCutPid::SetDetectorProbabilityRestriction(TString det, Int_t iPart, Double_t upperprob)
288{
289 //
290 // Sets the detector, the particle and the probability
291 // limit.
292 //
293
294 if(det.Contains("ITS")) fDetRestr = kITS;
295 if(det.Contains("TPC")) fDetRestr = kTPC;
296 if(det.Contains("TRD")) fDetRestr = kTRD;
297 if(det.Contains("TOF")) fDetRestr = kTOF;
298 if(det.Contains("HMPID")) fDetRestr = kHMPID;
299 fiPartRestr = iPart;
300 fDetProbRestr = upperprob;
301}
302//__________________________________
303void AliCFTrackCutPid::TrackInfo(const AliESDtrack *pTrk, ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
304{
305 //
306 // Loads the responses of the XXX chosen detectors for the pTrk
307 // and the corresponding trk status. The final trk status is also loaded.
308 //
309 if(fDets[kITS]) {
310 pTrk->GetITSpid(pid[kITS]);
311 status[kITS]=AliESDtrack::kITSpid;
312 }
313 if(fDets[kTPC]) {
314 pTrk->GetTPCpid(pid[kTPC]);
315 status[kTPC]=AliESDtrack::kTPCpid;
316 }
317 if(fDets[kTRD]) {
318 pTrk->GetTRDpid(pid[kTRD]);
319 status[kTRD]=AliESDtrack::kTRDpid;
320 }
321 if(fDets[kTOF]) {
322 pTrk->GetTOFpid(pid[kTOF]);
323 status[kTOF]=AliESDtrack::kTOFpid;
324 }
325 if(fDets[kHMPID]) {
326 pTrk->GetHMPIDpid(pid[kHMPID]);
327 status[kHMPID]=AliESDtrack::kHMPIDpid;
328 }
329 if(fDetRestr>=0){
330 if(fDetRestr == kITS) pTrk->GetITSpid(pid[kITS]);
331 if(fDetRestr == kTPC) pTrk->GetITSpid(pid[kTPC]);
332 if(fDetRestr == kTRD) pTrk->GetITSpid(pid[kTRD]);
333 if(fDetRestr == kTOF) pTrk->GetITSpid(pid[kTOF]);
334 if(fDetRestr == kHMPID) pTrk->GetITSpid(pid[kHMPID]);
335 }
336
337 status[kNdets]=pTrk->GetStatus();
338 pTrk->GetESDpid(pid[kNdets]);
339}
340//__________________________________
341void AliCFTrackCutPid::SetPPriors(AliESDtrack *pTrk)
342{
343 //
344 //sets the mommentum dependent a priori concentrations
345 //
346
347 for(Int_t i=0; i< AliPID::kSPECIES; i++) {
348 if(pTrk->P()>fPriorsFunc[i]->GetXmin() && pTrk->P() < fPriorsFunc[i]->GetXmax()) fPriors[i]=fPriorsFunc[i]->Eval(pTrk->P());
349 else {AliInfo("the track momentum is not in the function range. Priors are equal") fPriors[i] = 0.2;}
350 }
351}
352//______________________________________
353ULong_t AliCFTrackCutPid::StatusForAND(ULong_t status[kNdets+1]) const
354{
355 //
356 //In case of AND of more detectors the AND-detector status combination.
357 //is calculated and also returned
358 //
359
360 ULong_t andstatus=0;
361 for(Int_t i=0; i< kNdets; i++) {
362 if(!fDetsInAnd[i]) continue;
363 andstatus = andstatus | status[i];
364 AliDebug(1,Form("trk status %lu %i AND-status combination: %lu",status[i],i,andstatus));
365 }
366 return andstatus;
367}
368//_______________________________________
369Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES]) const
370{
371 // Identifies the track if its probability is higher than the cut
372 // value. The default value is fCut=0.2, therefore the most probable
373 // one is identified by default. Here all the checks on how to identify
374 // the track are performed (single detector or combined PID,..., detector
375 // restriction probability
376 // Returns: integer corresponding to the identified particle
377
378 Int_t iPart=-1;
379
380 if(!fgIsComb){
381 Bool_t isDet=kFALSE;
382 for(Int_t i=0; i<kNdets; i++){
383 if(!fDets[i]) continue;
384 isDet=kTRUE;
385 AliDebug(1,Form("trk status %lu %i-det-pid status %lu -> combination: %lu",status[kNdets],i,status[i],status[kNdets]&status[i]));
386 if(!(status[kNdets]&status[i])){
387 iPart=-10;
388 AliDebug(1,Form("detector %i -> pid trk status not ok",i));
389 }
390 else {
391 AliDebug(1,Form("resp : %f %f %f %f %f",pid[i][0],pid[i][1],pid[i][2],pid[i][3],pid[i][4]));
392 if(fIsQAOn) iPart = IdentifyQA(pid[i],i);
393 else iPart = Identify(pid[i]);
394 }
395 }
396 if(!isDet){
397 AliDebug(1,Form(" !!! No detector selected, the ESD-pid response is considered"));
398 iPart = Identify(pid[kNdets]);
399 }
400 }else{
401 Double_t calcprob[5];
f9658e11 402 CombPID(status,pid,calcprob);
563113d0 403 iPart = Identify(calcprob);
f9658e11 404
563113d0 405 }
406
407
408 AliDebug(1,Form("selected particle: %i",iPart));
409
410 if(iPart >=0 && fiPartRestr>=0) {
411 AliPID restr(pid[fDetRestr]);
412 restr.SetPriors(fPriors);
413 AliDebug(1,Form("setted upper limit: %f det %i : probability %f ",fDetProbRestr,fDetRestr,restr.GetProbability((AliPID::EParticleType)fiPartRestr)));
414 if(restr.GetProbability((AliPID::EParticleType)fiPartRestr) > fDetProbRestr) {
415 iPart = kDetRestr;
416 AliDebug(1,"\n\n the detector restrictions refused the ID \n\n");
417 }
418 }//cross checks with one detector probability
419
420 AliDebug(1,Form("after the check the selected particle is %i",iPart));
421
422 return iPart;
423}
bc6176fb 424//_________________________________________________________________________________
425Int_t AliCFTrackCutPid::GetAODID(AliAODTrack *aodtrack) const
426{
427//
428// Identifies the AOD Track using the combined pid responses
429//
430
431 Double_t combpid[AliPID::kSPECIES];
f9658e11 432 for(Int_t i=0; i< AliPID::kSPECIES; i++) {
433 combpid[i]= aodtrack->PID()[i];
434 if(!fhCombResp[i]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",i));
435 else fhCombResp[i]->Fill(combpid[i]);
436 }
bc6176fb 437 return Identify(combpid);
438}
563113d0 439//__________________________________
440Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const
441{
442 // Checks if there are no equal values and if the valus
443 // difference between the selected particle and the others i
444 // is higher than a lower limit.
445 // Returns: kTRUE= is acceptable
446
7f0c5c1e 447 AliDebug(2,Form("input particle: %i",iPsel));
563113d0 448 Bool_t ck=kTRUE;
449
450 if(iPsel<0) ck=kFALSE;
451
452 else {
453 for(Int_t j=0; j< AliPID::kSPECIES; j++) {
454 if(j!=iPsel && TMath::Abs(p[j]-p[iPsel])<minDiff) ck=kFALSE;
455 }
456 if(!ck) AliDebug(1,"the values are too close ");
457 }
458 return ck;
459}
460//___________________________________________
461Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const
462{
463 //
464 // The identification is actually performed here with possible
465 // checks on the det responses and/or probabilities
466 //
bc6176fb 467
563113d0 468 Int_t iPart = -1;
bc6176fb 469
470 AliDebug(2,Form("calc response bef: %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
f9658e11 471 AliDebug(2,Form("priors : %f %f %f %f %f",fPriors[0],fPriors[1],fPriors[2],fPriors[3],fPriors[4]));
472
563113d0 473 AliPID getpid(pid,kTRUE);
bc6176fb 474 getpid.SetPriors(fPriors);
563113d0 475
bc6176fb 476 Double_t probability[AliPID::kSPECIES]={0.,0.,0.,0.,0.};
477 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) {
478 probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
479 AliDebug(2,Form("prob %i %f",iP, probability[iP]));
f9658e11 480 if(fIsQAOn) fhCombProb[iP]->Fill(probability[iP]);
bc6176fb 481 }
563113d0 482
483 AliPID::EParticleType sel = getpid.GetMostProbable();
563113d0 484 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
7f0c5c1e 485 AliDebug(2,Form("probabilities : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
486
563113d0 487 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
488 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
7f0c5c1e 489
563113d0 490 return iPart;
491}
492//___________________________________________
493Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t idets) const
494{
495 //
496 // The same as Identify, but with the QA histo filling
497 //
498 //
499
500 Int_t iPart = -1;
501 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
502
503 AliPID getpid(pid,kTRUE);
f9658e11 504 getpid.SetPriors(fPriors);
505
563113d0 506 AliPID::EParticleType sel = getpid.GetMostProbable();
507 Double_t probability[AliPID::kSPECIES];
508 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) {
509 probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
510 fhProb[idets][iP]->Fill(probability[iP]);
511 }
512
513 AliPID toresp(pid,kTRUE); Double_t qapriors[5]={0.2,0.2,0.2,0.2,0.2};
514 toresp.SetPriors(qapriors);
515 for(Int_t iPr=0; iPr<AliPID::kSPECIES; iPr++) fhResp[idets][iPr]->Fill(toresp.GetProbability((AliPID::EParticleType)iPr));
516
517 if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel;
518 AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
519 AliDebug(1,Form("probab : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4]));
520 if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
521 if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
522 return iPart;
523}
524//___________________________________________
525Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
526 //
527 // method for the pid-cut selction
528 //
bc6176fb 529 Bool_t sel = kFALSE;
563113d0 530
531 if (!track) return kFALSE ;
532 TString className(track->ClassName());
bc6176fb 533 if (className.CompareTo("AliESDtrack") == 0) {
f9658e11 534 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
563113d0 535 ULong_t status[kNdets+1]={0,0,0,0,0,0};
536 Double_t pid[kNdets+1][AliPID::kSPECIES];
537 TrackInfo(esdTrack,status,pid);
538 if(fIsPpriors) SetPPriors(esdTrack);
bc6176fb 539 if(GetID(status,pid)==fgParticleType) sel = kTRUE;
540 }
541
542 if (className.CompareTo("AliAODTrack") == 0) {
f9658e11 543 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
bc6176fb 544 if(GetAODID(aodtrack) == fgParticleType) sel = kTRUE;
545 }
546
547 return sel;
548
563113d0 549}
550//__________________________________
551void AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const
552{
553 //
554 // Calculates the combined PID according to the chosen detectors.
555 // and provides the array of probabilities
556 //
557
558 Bool_t isdet=kFALSE;
563113d0 559 Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.};
7f0c5c1e 560
563113d0 561 ULong_t andstatus =0;
562 if(fIsDetAND) {
563 andstatus = StatusForAND(status);
564 AliDebug(1,Form("AND combination %lu",andstatus));
565 }
f9658e11 566
567 //Products of single detector responses
563113d0 568 for(Int_t j=0; j<AliPID::kSPECIES; j++){
569 for(Int_t i=0; i< kNdets; i++){
570 if(!fDets[i]) continue;
571 if(status[kNdets]&status[i]) {
572 if(fIsDetAND) {
573 ULong_t checkstatus = status[kNdets]&andstatus;
574 if(checkstatus != andstatus) continue;
575 else {
576 prod[j]*=pid[i][j];
577 isdet = kTRUE;
578 AliDebug(1,Form("-----> trk status %lu and status %lu -> trk-ANDdetector status combination %lu",status[kNdets],andstatus,status[kNdets]&andstatus));
7f0c5c1e 579 AliDebug(1,Form("In det %i -> particle %i response is %f",i,j,pid[i][j]));
563113d0 580 }
f9658e11 581 } else {
563113d0 582 prod[j]*=pid[i][j];
583 isdet=kTRUE;
7f0c5c1e 584 AliDebug(2,Form("In det %i -> particle %i response is %f",i,j,pid[i][j]));
f9658e11 585
586 if(fIsQAOn){
587 if(!fhResp[i][j]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",i,j));}
588 else fhResp[i][j]->Fill(pid[i][j]);
589
590 if(!fhProb[i][j]) {AliDebug(1,Form("no pointer to the histo fhProb%i%i, check if pidcut->Init() was called",i,j));}
591 else {
592 AliPID detprob(pid[i],kTRUE);
593 detprob.SetPriors(fPriors);
594 fhProb[i][j]->Fill(detprob.GetProbability((AliPID::EParticleType)j));
595 }
596 }
563113d0 597 }
598 }//combined mode
599 }//loop on dets
600 }//loop on species
7f0c5c1e 601
f9658e11 602
603 //no detectors found, then go to ESD pid...
563113d0 604 if(!isdet) {
7f0c5c1e 605 AliWarning("\n !! No detector found for the combined pid --> responses are from the ESDpid !! \n");
606 Double_t sumesdpid=0;
f9658e11 607 for(Int_t nn=0; nn<AliPID::kSPECIES; nn++) sumesdpid+=pid[kNdets][nn];
608 if(sumesdpid<=0) {
609 AliDebug(1,"priors or ESDpid are inconsistent, please check them");
610 return;
611 } else {
612 for(Int_t k=0; k<AliPID::kSPECIES; k++){
613 combpid[k] = pid[kNdets][k]/sumesdpid;
614 if(fIsQAOn) {
615 if(!fhCombResp[k]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",k));
616 else fhCombResp[k]->Fill(combpid[k]);
617 }
618 }//loop on species
619 }
620 return;
621 }
563113d0 622
f9658e11 623 Double_t add = 0; for(Int_t isumm=0; isumm<5; isumm++) add+=prod[isumm];
624 if(add>0) {
625 for(Int_t ip =0; ip < AliPID::kSPECIES; ip++) {
626 combpid[ip] = prod[ip]/add;
627 if(fIsQAOn) {
628 if(!fhCombResp[ip]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",ip));
629 else fhCombResp[ip]->Fill(combpid[ip]);
630
631 }
632 }
633 AliDebug(1,Form("calculated comb response: %f %f %f %f %f",combpid[0],combpid[1],combpid[2],combpid[3],combpid[4]));
634 } else {
635 AliDebug(1,"single detector responses are inconsistent, please check them....");
636 return;
637 }
638 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]));
639
563113d0 640}
641//__________________________________________
642//
643//QA part
644//_________________________________________
645void AliCFTrackCutPid::InitialiseHisto()
646{
647 //
648 //QA histo initialization
649 //
650 for(Int_t iP=0; iP<AliPID::kSPECIES; iP++){
651 fhCombResp[iP]=0x0;
652 fhCombProb[iP]=0x0;
653 for(Int_t iDet =0; iDet<kNdets; iDet++){
654 fhResp[iDet][iP]=0x0;
655 fhProb[iDet][iP]=0x0;
656 }
657 }
658}
659//______________________________________________
563113d0 660void AliCFTrackCutPid::DefineHistograms()
661{
662 //
663 //QA histo booking
664 //
7f0c5c1e 665
bc6176fb 666 if(fgIsAOD){
1f2d22de 667 const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"};
bc6176fb 668
669 for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
563113d0 670 {
ce2a1393 671 fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response (AODTrack) ",partic[iPart]),fNbins,fXmin,fXmax);
f9658e11 672 fhCombResp[iPart]->SetXTitle(Form(" %s combined response ",partic[iPart]));
673 fhCombResp[iPart]->SetYTitle("entries");
ce2a1393 674 AliDebug(1,Form( "%s is booked!!",fhCombResp[iPart]->GetName()));
675 fhCombProb[iPart] = new TH1F(Form("%s_pCombPart%i",GetName(),iPart),Form("%s combined probability (AODTrack) ",partic[iPart]),fNbins,fXmin,fXmax);
676 fhCombProb[iPart]->SetXTitle(Form(" %s combined probability ",partic[iPart]));
f9658e11 677 fhCombProb[iPart]->SetYTitle("entries");
ce2a1393 678 AliDebug(1,Form( "%s is booked!!",fhCombProb[iPart]->GetName()));
bc6176fb 679 }
680 }
681
682
683 else {
1f2d22de 684 const char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"};
685 const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"};
bc6176fb 686
687 for(Int_t iDet =0; iDet< kNdets; iDet++)
688 {
689 if(!fDets[iDet]) continue;
690 for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){
ce2a1393 691 fhResp[iDet][iP] = new TH1F(Form("%s_rDet%iPart%i",GetName(),iDet,iP),Form("%s response for %s ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
f9658e11 692 fhResp[iDet][iP]->SetXTitle(Form(" %s response ",partic[iP]));
693 fhResp[iDet][iP]->SetYTitle("entries");
ce2a1393 694 fhProb[iDet][iP] = new TH1F(Form("%s_pDet%iPart%i",GetName(),iDet,iP),Form("%s calculated probability for %s",detect[iDet],partic[iP]),fNbins,fXmin,fXmax);
695 fhProb[iDet][iP]->SetXTitle(Form(" %s probability ",partic[iP]));
f9658e11 696 fhProb[iDet][iP]->SetYTitle("entries");
bc6176fb 697 }
698 }
699
700
701 if(fgIsComb)
702 {
703 for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
704 {
ce2a1393 705 fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response ",partic[iPart]),fNbins,fXmin,fXmax);
f9658e11 706 fhCombResp[iPart]->SetXTitle(Form(" %s response ",partic[iPart]));
707 fhCombResp[iPart]->SetYTitle("entries");
ce2a1393 708 AliDebug(1,Form( "%s is booked!!",fhCombResp[iPart]->GetName()));
709 fhCombProb[iPart] = new TH1F(Form("%s_pCombPart%i",GetName(),iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax);
f9658e11 710 fhCombProb[iPart]->SetXTitle(Form(" %s response ",partic[iPart]));
711 fhCombProb[iPart]->SetYTitle("entries");
ce2a1393 712 AliDebug(1,Form( "%s is booked!!",fhCombProb[iPart]->GetName()));
bc6176fb 713 }
714 }
715 }
716
563113d0 717}
718//___________________________________________________
719
107a3100 720void AliCFTrackCutPid::AddQAHistograms(TList *qalist)
563113d0 721{
722 //
723 // adds QA histograms in a TList
724 //
725 if(!fIsQAOn) return;
107a3100 726 DefineHistograms();
563113d0 727
f9658e11 728 if(fgIsComb || fgIsAOD){
563113d0 729 for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){
730 qalist->Add(fhCombResp[iPart]);
731 qalist->Add(fhCombProb[iPart]);
732 }
733 }
734
735 for(Int_t iDet=0; iDet<kNdets; iDet++){
7f0c5c1e 736 if(!fDets[iDet]) continue;
563113d0 737 for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
7f0c5c1e 738 qalist->Add(fhResp[iDet][iP]);
739 qalist->Add(fhProb[iDet][iP]);
563113d0 740 }
741 }
7f0c5c1e 742
563113d0 743}
f9658e11 744