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