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