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