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