try to fix coverity 24179, 24178, 24125, 24124, check cast to pointers not being...
[u/mrichter/AliRoot.git] / CORRFW / AliCFTrackCutPid.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //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
60 #include "AliCFTrackCutPid.h"
61 #include "AliLog.h"
62 #include <TMath.h>
63 #include <TList.h>
64
65 ClassImp(AliCFTrackCutPid) 
66
67 //__________________________________________________________________________________
68 // CUT ON TRACK PID
69 //__________________________________________________________________________________
70 AliCFTrackCutPid::AliCFTrackCutPid() :
71   AliCFCutBase(),
72   fCut(0.2),
73   fMinDiffResponse(0.0001),
74   fMinDiffProbability(0.001),
75   fgParticleType(10),
76   fgIsComb(kTRUE),
77   fgIsAOD(kFALSE),
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   fProbThreshold(0.)
89
90   //
91   //Default constructor 
92   //
93   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
94     fPriors[j]=0.2;
95   }
96   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
97     fPriorsFunc[j]=0x0;
98   }
99   for(Int_t jDet=0; jDet< kNdets; jDet++)  {
100     fDets[jDet]=kFALSE;
101     fDetsInAnd[jDet]=kFALSE;
102   }
103   
104   InitialiseHisto();
105 }
106 //______________________________
107 AliCFTrackCutPid::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),
114   fgIsAOD(kFALSE),
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),
124   fDetProbRestr(1),
125   fProbThreshold(0.)
126 {
127   //
128   //Constructor
129   // 
130   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
131     fPriors[j]=0.2;
132   }
133   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
134     fPriorsFunc[j]=0x0;
135   }
136   for(Int_t jDet=0; jDet< kNdets; jDet++)  {
137     fDets[jDet]=kFALSE;
138     fDetsInAnd[jDet]=kFALSE;
139   }
140   
141   InitialiseHisto();
142 }
143 //______________________________
144 AliCFTrackCutPid::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),
151   fgIsAOD(c.fgIsAOD),
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),
161   fDetProbRestr(c.fDetProbRestr),
162   fProbThreshold(c.fProbThreshold)
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   }
175   for(Int_t j=0; j< AliPID::kSPECIES; j++){
176     fPriors[j]=c.fPriors[j];
177   }
178   for(Int_t j=0; j< AliPID::kSPECIES; j++){
179     fPriorsFunc[j]=c.fPriorsFunc[j];
180     fhCombResp[j]=c.fhCombResp[j];
181     fhCombProb[j]=c.fhCombProb[j];
182   }
183 }
184 //______________________________
185 AliCFTrackCutPid& 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;
194     this->fMinDiffProbability=c.fMinDiffProbability;
195     this->fgParticleType=c.fgParticleType;  
196     this->fgIsComb=c.fgIsComb;
197     this->fgIsAOD=c.fgIsAOD;
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;
208     this->fProbThreshold=c.fProbThreshold;
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
218     for(Int_t j=0; j< AliPID::kSPECIES; j++){
219       this->fPriors[j]=c.fPriors[j];
220     }
221     for(Int_t j=0; j< AliPID::kSPECIES; j++){
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 //__________________________________________________________________________________
231 AliCFTrackCutPid::~AliCFTrackCutPid() {
232   //
233   //dtor
234   //
235
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   }
248 }
249 //__________________________________
250 void 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;}
261
262   if(dets.Contains("ALL")) for(Int_t i=0; i< kNdets ; i++) fDets[i]=kTRUE;
263 }
264 //__________________________________
265 void 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 //__________________________________
273 void 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 //____________________________________________
283 void 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 //
297 void 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 //__________________________________
313 void 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 //__________________________________
351 void 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());
359     else {AliInfo("the track momentum is not in the function range. Priors are equal"); fPriors[i] = 0.2;}   
360   }
361 }   
362 //______________________________________
363 ULong_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 //_______________________________________
379 Int_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];
412     CombPID(status,pid,calcprob); 
413     iPart = Identify(calcprob);
414     
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 }
434 //_________________________________________________________________________________
435 Int_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];
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   }
447   return Identify(combpid);
448 }
449 //__________________________________
450 Bool_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   
457   AliDebug(2,Form("input particle: %i",iPsel));
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 //___________________________________________
471 Int_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   //
477
478   Int_t iPart = -1;
479  
480   AliDebug(2,Form("calc response bef: %f  %f  %f  %f  %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
481   AliDebug(2,Form("priors           : %f  %f  %f  %f  %f",fPriors[0],fPriors[1],fPriors[2],fPriors[3],fPriors[4]));
482
483   AliPID getpid(pid,kTRUE);
484   getpid.SetPriors(fPriors); 
485   
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]));
490     if(fIsQAOn) fhCombProb[iP]->Fill(probability[iP]);
491   }
492   
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   }
502
503   if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp;
504   if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb;
505
506   return iPart;
507 }
508 //___________________________________________
509 Int_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);
520   getpid.SetPriors(fPriors);
521       
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   
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};
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 //___________________________________________
542 Bool_t AliCFTrackCutPid::IsSelected(TObject *track){
543   //
544   //  method for the pid-cut selction
545   //
546   Bool_t sel = kFALSE;
547   
548   if (!track) return kFALSE ;
549   TString className(track->ClassName());
550   if (className.CompareTo("AliESDtrack") == 0) {
551   AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); 
552   if (!esdTrack) return kFALSE;
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);
557   if(GetID(status,pid)==fgParticleType) sel = kTRUE;
558   }
559
560   if (className.CompareTo("AliAODTrack") == 0) {
561   AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
562   if (!aodtrack) return kFALSE ;
563   if(GetAODID(aodtrack) == fgParticleType) sel = kTRUE;
564   }
565
566  return sel;
567
568 }
569 //__________________________________
570 void  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;
578   Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.};
579   
580   ULong_t andstatus =0;
581   if(fIsDetAND) {
582     andstatus = StatusForAND(status);
583     AliDebug(1,Form("AND combination %lu",andstatus));
584   } 
585   
586   //Products of single detector responses
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));
598             AliDebug(1,Form("In det %i ->  particle %i response is %f",i,j,pid[i][j]));
599           }
600         } else {
601           prod[j]*=pid[i][j];
602           isdet=kTRUE;
603           AliDebug(2,Form("In det %i ->  particle %i response is %f",i,j,pid[i][j]));
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            }
616         }
617       }//combined mode
618     }//loop on dets
619   }//loop on species
620
621   
622    //no detectors found, then go to ESD pid...
623   if(!isdet) {
624     AliWarning("\n !! No detector found for the combined pid --> responses are from the ESDpid !! \n");
625     Double_t sumesdpid=0;
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     }   
641     
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  
659 }
660 //__________________________________________
661 //
662 //QA part
663 //_________________________________________
664 void 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 //______________________________________________
679 void AliCFTrackCutPid::DefineHistograms()
680 {
681   //
682   //QA histo booking
683   //
684
685   if(fgIsAOD){
686      const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"}; 
687
688      for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
689       {
690        fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response  (AODTrack)  ",partic[iPart]),fNbins,fXmin,fXmax);
691        fhCombResp[iPart]->SetXTitle(Form(" %s combined response ",partic[iPart]));
692        fhCombResp[iPart]->SetYTitle("entries");
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]));
696        fhCombProb[iPart]->SetYTitle("entries");
697        AliDebug(1,Form(  "%s is booked!!",fhCombProb[iPart]->GetName()));
698      }
699    }
700
701
702   else {
703    const char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"};
704    const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"};
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++){
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);
711         fhResp[iDet][iP]->SetXTitle(Form(" %s response ",partic[iP]));
712         fhResp[iDet][iP]->SetYTitle("entries");
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]));
715         fhProb[iDet][iP]->SetYTitle("entries");
716        }
717      }
718     
719
720   if(fgIsComb)
721      { 
722       for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++)
723         {
724           fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response    ",partic[iPart]),fNbins,fXmin,fXmax);
725           fhCombResp[iPart]->SetXTitle(Form(" %s response ",partic[iPart]));
726           fhCombResp[iPart]->SetYTitle("entries");
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);
729           fhCombProb[iPart]->SetXTitle(Form(" %s response ",partic[iPart]));
730           fhCombProb[iPart]->SetYTitle("entries");
731           AliDebug(1,Form(  "%s is booked!!",fhCombProb[iPart]->GetName()));
732         }
733      }
734    }
735
736 }
737 //___________________________________________________
738
739 void AliCFTrackCutPid::AddQAHistograms(TList *qalist) 
740 {
741   //
742   // adds QA histograms in a TList
743   //
744   if(!fIsQAOn) return;
745   DefineHistograms();
746   
747   if(fgIsComb || fgIsAOD){
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++){
755     if(!fDets[iDet]) continue;
756     for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
757       qalist->Add(fhResp[iDet][iP]);
758       qalist->Add(fhProb[iDet][iP]);
759     }
760   }  
761
762 }
763