Updates for the pid caching. By default it is still switched off
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDResponse.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 /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */
17
18 //-----------------------------------------------------------------
19 //        Base class for handling the pid response               //
20 //        functions of all detectors                             //
21 //        and give access to the nsigmas                         //
22 //                                                               //
23 //   Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch //
24 //-----------------------------------------------------------------
25
26 #include <TList.h>
27 #include <TObjArray.h>
28 #include <TPRegexp.h>
29 #include <TF1.h>
30 #include <TSpline.h>
31 #include <TFile.h>
32 #include <TArrayI.h>
33 #include <TArrayF.h>
34
35 #include <AliVEvent.h>
36 #include <AliVTrack.h>
37 #include <AliLog.h>
38 #include <AliPID.h>
39 #include <AliOADBContainer.h>
40 #include <AliTRDPIDResponseObject.h>
41 #include <AliTOFPIDParams.h>
42
43 #include "AliPIDResponse.h"
44 #include "AliDetectorPID.h"
45
46 #include "AliCentrality.h"
47
48 ClassImp(AliPIDResponse);
49
50 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
51 TNamed("PIDResponse","PIDResponse"),
52 fITSResponse(isMC),
53 fTPCResponse(),
54 fTRDResponse(),
55 fTOFResponse(),
56 fEMCALResponse(),
57 fRange(5.),
58 fITSPIDmethod(kITSTruncMean),
59 fIsMC(isMC),
60 fCachePID(kTRUE),
61 fOADBPath(),
62 fCustomTPCpidResponse(),
63 fBeamType("PP"),
64 fLHCperiod(),
65 fMCperiodTPC(),
66 fMCperiodUser(),
67 fCurrentFile(),
68 fRecoPass(0),
69 fRecoPassUser(-1),
70 fRun(0),
71 fOldRun(0),
72 fResT0A(75.),
73 fResT0C(65.),
74 fResT0AC(55.),
75 fArrPidResponseMaster(NULL),
76 fResolutionCorrection(NULL),
77 fOADBvoltageMaps(NULL),
78 fTRDPIDResponseObject(NULL),
79 fTOFtail(1.1),
80 fTOFPIDParams(NULL),
81 fEMCALPIDParams(NULL),
82 fCurrentEvent(NULL),
83 fCurrCentrality(0.0),
84 fTuneMConData(kFALSE)
85 {
86   //
87   // default ctor
88   //
89   AliLog::SetClassDebugLevel("AliPIDResponse",0);
90   AliLog::SetClassDebugLevel("AliESDpid",0);
91   AliLog::SetClassDebugLevel("AliAODpidUtil",0);
92
93 }
94
95 //______________________________________________________________________________
96 AliPIDResponse::~AliPIDResponse()
97 {
98   //
99   // dtor
100   //
101   delete fArrPidResponseMaster;
102   delete fTRDPIDResponseObject;
103   delete fTOFPIDParams;
104 }
105
106 //______________________________________________________________________________
107 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
108 TNamed(other),
109 fITSResponse(other.fITSResponse),
110 fTPCResponse(other.fTPCResponse),
111 fTRDResponse(other.fTRDResponse),
112 fTOFResponse(other.fTOFResponse),
113 fEMCALResponse(other.fEMCALResponse),
114 fRange(other.fRange),
115 fITSPIDmethod(other.fITSPIDmethod),
116 fIsMC(other.fIsMC),
117 fCachePID(other.fCachePID),
118 fOADBPath(other.fOADBPath),
119 fCustomTPCpidResponse(other.fCustomTPCpidResponse),
120 fBeamType("PP"),
121 fLHCperiod(),
122 fMCperiodTPC(),
123 fMCperiodUser(other.fMCperiodUser),
124 fCurrentFile(),
125 fRecoPass(0),
126 fRecoPassUser(other.fRecoPassUser),
127 fRun(0),
128 fOldRun(0),
129 fResT0A(75.),
130 fResT0C(65.),
131 fResT0AC(55.),
132 fArrPidResponseMaster(NULL),
133 fResolutionCorrection(NULL),
134 fOADBvoltageMaps(NULL),
135 fTRDPIDResponseObject(NULL),
136 fTOFtail(1.1),
137 fTOFPIDParams(NULL),
138 fEMCALPIDParams(NULL),
139 fCurrentEvent(NULL),
140 fCurrCentrality(0.0),
141 fTuneMConData(kFALSE)
142 {
143   //
144   // copy ctor
145   //
146 }
147
148 //______________________________________________________________________________
149 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
150 {
151   //
152   // copy ctor
153   //
154   if(this!=&other) {
155     delete fArrPidResponseMaster;
156     TNamed::operator=(other);
157     fITSResponse=other.fITSResponse;
158     fTPCResponse=other.fTPCResponse;
159     fTRDResponse=other.fTRDResponse;
160     fTOFResponse=other.fTOFResponse;
161     fEMCALResponse=other.fEMCALResponse;
162     fRange=other.fRange;
163     fITSPIDmethod=other.fITSPIDmethod;
164     fOADBPath=other.fOADBPath;
165     fCustomTPCpidResponse=other.fCustomTPCpidResponse;
166     fIsMC=other.fIsMC;
167     fCachePID=other.fCachePID;
168     fBeamType="PP";
169     fLHCperiod="";
170     fMCperiodTPC="";
171     fMCperiodUser=other.fMCperiodUser;
172     fCurrentFile="";
173     fRecoPass=0;
174     fRecoPassUser=other.fRecoPassUser;
175     fRun=0;
176     fOldRun=0;
177     fResT0A=75.;
178     fResT0C=65.;
179     fResT0AC=55.;
180     fArrPidResponseMaster=NULL;
181     fResolutionCorrection=NULL;
182     fOADBvoltageMaps=NULL;
183     fTRDPIDResponseObject=NULL;
184     fEMCALPIDParams=NULL;
185     fTOFtail=1.1;
186     fTOFPIDParams=NULL;
187     fCurrentEvent=other.fCurrentEvent;
188
189   }
190   return *this;
191 }
192
193 //______________________________________________________________________________
194 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
195 {
196   //
197   // NumberOfSigmas for 'detCode'
198   //
199
200   switch (detCode){
201     case kDetITS: return NumberOfSigmasITS(track, type); break;
202     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
203     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
204     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
205     default: return -999.;
206   }
207
208 }
209
210 //______________________________________________________________________________
211 Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
212 {
213   //
214   // NumberOfSigmas for 'detCode'
215   //
216   return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
217 }
218
219 //______________________________________________________________________________
220 // public buffered versions of the PID calculation
221 //
222
223 //______________________________________________________________________________
224 Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
225 {
226   //
227   // Calculate the number of sigmas in the ITS
228   //
229   
230   AliVTrack *track=(AliVTrack*)vtrack;
231   
232   // look for cached value first
233   const AliDetectorPID *detPID=track->GetDetectorPID();
234   const EDetector detector=kITS;
235
236   if ( detPID && detPID->HasNumberOfSigmas(detector)){
237     return detPID->GetNumberOfSigmas(detector, type);
238   } else if (fCachePID) {
239     FillTrackDetectorPID(track, detector);
240     detPID=track->GetDetectorPID();
241     return detPID->GetNumberOfSigmas(detector, type);
242   }
243
244   return GetNumberOfSigmasITS(track, type);
245 }
246
247 //______________________________________________________________________________
248 Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
249 {
250   //
251   // Calculate the number of sigmas in the TPC
252   //
253   
254   AliVTrack *track=(AliVTrack*)vtrack;
255   
256   // look for cached value first
257   const AliDetectorPID *detPID=track->GetDetectorPID();
258   const EDetector detector=kTPC;
259   
260   if ( detPID && detPID->HasNumberOfSigmas(detector)){
261     return detPID->GetNumberOfSigmas(detector, type);
262   } else if (fCachePID) {
263     FillTrackDetectorPID(track, detector);
264     detPID=track->GetDetectorPID();
265     return detPID->GetNumberOfSigmas(detector, type);
266   }
267   
268   return GetNumberOfSigmasTPC(track, type);
269 }
270
271 //______________________________________________________________________________
272 Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack, 
273                                            AliPID::EParticleType type,
274                                            AliTPCPIDResponse::ETPCdEdxSource dedxSource) 
275 {
276   //get number of sigmas according the selected TPC gain configuration scenario
277   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
278
279   Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource);
280
281   return nSigma;
282 }
283
284 //______________________________________________________________________________
285 Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
286 {
287   //
288   // Calculate the number of sigmas in the TOF
289   //
290   
291   AliVTrack *track=(AliVTrack*)vtrack;
292   
293   // look for cached value first
294   const AliDetectorPID *detPID=track->GetDetectorPID();
295   const EDetector detector=kTOF;
296   
297   if ( detPID && detPID->HasNumberOfSigmas(detector)){
298     return detPID->GetNumberOfSigmas(detector, type);
299   } else if (fCachePID) {
300     FillTrackDetectorPID(track, detector);
301     detPID=track->GetDetectorPID();
302     return detPID->GetNumberOfSigmas(detector, type);
303   }
304   
305   return GetNumberOfSigmasTOF(track, type);
306 }
307
308 //______________________________________________________________________________
309 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
310 {
311   //
312   // Calculate the number of sigmas in the EMCAL
313   //
314   
315   AliVTrack *track=(AliVTrack*)vtrack;
316
317   // look for cached value first
318   const AliDetectorPID *detPID=track->GetDetectorPID();
319   const EDetector detector=kEMCAL;
320   
321   if ( detPID && detPID->HasNumberOfSigmas(detector)){
322     return detPID->GetNumberOfSigmas(detector, type);
323   } else if (fCachePID) {
324     FillTrackDetectorPID(track, detector);
325     detPID=track->GetDetectorPID();
326     return detPID->GetNumberOfSigmas(detector, type);
327   }
328   
329   return GetNumberOfSigmasEMCAL(track, type);
330 }
331
332 //______________________________________________________________________________
333 Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4])  const
334 {
335   //
336   // emcal nsigma with eop and showershape
337   //
338   AliVTrack *track=(AliVTrack*)vtrack;
339   
340   AliVCluster *matchedClus = NULL;
341
342   Double_t mom     = -1.; 
343   Double_t pt      = -1.; 
344   Double_t EovP    = -1.;
345   Double_t fClsE   = -1.;
346
347   // initialize eop and shower shape parameters
348   eop = -1.;
349   for(Int_t i = 0; i < 4; i++){
350     showershape[i] = -1.;
351   }
352   
353   Int_t nMatchClus = -1;
354   Int_t charge     = 0;
355   
356   // Track matching
357   nMatchClus = track->GetEMCALcluster();
358   if(nMatchClus > -1){
359
360     mom    = track->P();
361     pt     = track->Pt();
362     charge = track->Charge();
363     
364     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
365     
366     if(matchedClus){
367       
368       // matched cluster is EMCAL
369       if(matchedClus->IsEMCAL()){
370         
371         fClsE       = matchedClus->E();
372         EovP        = fClsE/mom;
373         
374         // fill used EMCAL variables here
375         eop            = EovP; // E/p
376         showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
377         showershape[1] = matchedClus->GetM02(); // long axis
378         showershape[2] = matchedClus->GetM20(); // short axis
379         showershape[3] = matchedClus->GetDispersion(); // dispersion
380
381         // look for cached value first
382         const AliDetectorPID *detPID=track->GetDetectorPID();
383         const EDetector detector=kEMCAL;
384         
385         if ( detPID && detPID->HasNumberOfSigmas(detector)){
386           return detPID->GetNumberOfSigmas(detector, type);
387         } else if (fCachePID) {
388           FillTrackDetectorPID(track, detector);
389           detPID=track->GetDetectorPID();
390           return detPID->GetNumberOfSigmas(detector, type);
391         }
392         
393         // NSigma value really meaningful only for electrons!
394         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
395       }
396     }
397   }
398   return -999;
399 }
400
401
402 //______________________________________________________________________________
403 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[])  const
404 {
405   //
406   // Compute PID response of 'detCode'
407   //
408
409   switch (detCode){
410     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
411     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
412     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
413     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
414     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
415     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
416     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
417     default: return kDetNoSignal;
418   }
419 }
420
421 //______________________________________________________________________________
422 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
423 {
424   //
425   // Compute PID response of 'detCode'
426   //
427
428   return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
429 }
430
431 //______________________________________________________________________________
432 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
433 {
434   //
435   // Compute PID response for the ITS
436   //
437
438   // look for cached value first
439   const AliDetectorPID *detPID=track->GetDetectorPID();
440   const EDetector detector=kITS;
441   
442   if ( detPID && detPID->HasRawProbabilitiy(detector)){
443     return detPID->GetRawProbability(detector, p, nSpecies);
444   } else if (fCachePID) {
445     FillTrackDetectorPID(track, detector);
446     detPID=track->GetDetectorPID();
447     return detPID->GetRawProbability(detector, p, nSpecies);
448   }
449   
450   return GetComputeITSProbability(track, nSpecies, p);
451 }
452 //______________________________________________________________________________
453 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
454 {
455   //
456   // Compute PID response for the TPC
457   //
458   
459   // look for cached value first
460   const AliDetectorPID *detPID=track->GetDetectorPID();
461   const EDetector detector=kTPC;
462   
463   if ( detPID && detPID->HasRawProbabilitiy(detector)){
464     return detPID->GetRawProbability(detector, p, nSpecies);
465   } else if (fCachePID) {
466     FillTrackDetectorPID(track, detector);
467     detPID=track->GetDetectorPID();
468     return detPID->GetRawProbability(detector, p, nSpecies);
469   }
470   
471   return GetComputeTPCProbability(track, nSpecies, p);
472 }
473 //______________________________________________________________________________
474 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
475 {
476   //
477   // Compute PID response for the
478   //
479   
480   const AliDetectorPID *detPID=track->GetDetectorPID();
481   const EDetector detector=kTOF;
482   
483   if ( detPID && detPID->HasRawProbabilitiy(detector)){
484     return detPID->GetRawProbability(detector, p, nSpecies);
485   } else if (fCachePID) {
486     FillTrackDetectorPID(track, detector);
487     detPID=track->GetDetectorPID();
488     return detPID->GetRawProbability(detector, p, nSpecies);
489   }
490   
491   return GetComputeTOFProbability(track, nSpecies, p);
492 }
493 //______________________________________________________________________________
494 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
495 {
496   //
497   // Compute PID response for the
498   //
499   
500   const AliDetectorPID *detPID=track->GetDetectorPID();
501   const EDetector detector=kTRD;
502
503   // chacke only for the default method (1d at the moment)
504   if (PIDmethod!=AliTRDPIDResponse::kLQ1D) return GetComputeTRDProbability(track, nSpecies, p, PIDmethod);
505   
506   if ( detPID && detPID->HasRawProbabilitiy(detector)){
507     return detPID->GetRawProbability(detector, p, nSpecies);
508   } else if (fCachePID) {
509     FillTrackDetectorPID(track, detector);
510     detPID=track->GetDetectorPID();
511     return detPID->GetRawProbability(detector, p, nSpecies);
512   }
513   
514   return GetComputeTRDProbability(track, nSpecies, p);
515 }
516 //______________________________________________________________________________
517 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
518 {
519   //
520   // Compute PID response for the EMCAL
521   //
522   
523   const AliDetectorPID *detPID=track->GetDetectorPID();
524   const EDetector detector=kEMCAL;
525   
526   if ( detPID && detPID->HasRawProbabilitiy(detector)){
527     return detPID->GetRawProbability(detector, p, nSpecies);
528   } else if (fCachePID) {
529     FillTrackDetectorPID(track, detector);
530     detPID=track->GetDetectorPID();
531     return detPID->GetRawProbability(detector, p, nSpecies);
532   }
533   
534   return GetComputeEMCALProbability(track, nSpecies, p);
535 }
536 //______________________________________________________________________________
537 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
538 {
539   //
540   // Compute PID response for the PHOS
541   //
542   
543   // look for cached value first
544 //   if (track->GetDetectorPID()){
545 //     return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies);
546 //   }
547   
548   // set flat distribution (no decision)
549   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
550   return kDetNoSignal;
551 }
552 //______________________________________________________________________________
553 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
554 {
555   //
556   // Compute PID response for the HMPID
557   //
558
559   const AliDetectorPID *detPID=track->GetDetectorPID();
560   const EDetector detector=kHMPID;
561
562   if ( detPID && detPID->HasRawProbabilitiy(detector)){
563     return detPID->GetRawProbability(detector, p, nSpecies);
564   } else if (fCachePID) {
565     FillTrackDetectorPID(track, detector);
566     detPID=track->GetDetectorPID();
567     return detPID->GetRawProbability(detector, p, nSpecies);
568   }
569
570   return GetComputeHMPIDProbability(track, nSpecies, p);
571 }
572
573 //______________________________________________________________________________
574 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
575 {
576   //
577   // Apply settings for the current event
578   //
579   fRecoPass=pass;
580   
581
582   fCurrentEvent=NULL;
583   if (!event) return;
584   fCurrentEvent=event;
585   if (run>0) fRun=run;
586   else fRun=event->GetRunNumber();
587   
588   if (fRun!=fOldRun){
589     ExecNewRun();
590     fOldRun=fRun;
591   }
592   
593   //TPC resolution parametrisation PbPb
594   if ( fResolutionCorrection ){
595     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
596     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
597   }
598   
599   //TOF resolution
600   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
601
602
603   // Get and set centrality
604   AliCentrality *centrality = event->GetCentrality();
605   if(centrality){
606     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
607   }
608   else{
609     fCurrCentrality = -1;
610   }
611 }
612
613 //______________________________________________________________________________
614 void AliPIDResponse::ExecNewRun()
615 {
616   //
617   // Things to Execute upon a new run
618   //
619   SetRecoInfo();
620   
621   SetITSParametrisation();
622   
623   SetTPCPidResponseMaster();
624   SetTPCParametrisation();
625
626   SetTRDPidResponseMaster(); 
627   InitializeTRDResponse();
628
629   SetEMCALPidResponseMaster(); 
630   InitializeEMCALResponse();
631   
632   SetTOFPidResponseMaster();
633   InitializeTOFResponse();
634
635   if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
636 }
637
638 //______________________________________________________________________________
639 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
640 {
641   //
642   // Get TPC multiplicity in bins of 150
643   //
644   
645   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
646   Double_t tpcMulti=0.;
647   if(vertexTPC){
648     Double_t vertexContribTPC=vertexTPC->GetNContributors();
649     tpcMulti=vertexContribTPC/150.;
650     if (tpcMulti>20.) tpcMulti=20.;
651   }
652   
653   return tpcMulti;
654 }
655
656 //______________________________________________________________________________
657 void AliPIDResponse::SetRecoInfo()
658 {
659   //
660   // Set reconstruction information
661   //
662   
663   //reset information
664   fLHCperiod="";
665   fMCperiodTPC="";
666   
667   fBeamType="";
668     
669   fBeamType="PP";
670   
671   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
672   TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
673
674   //find the period by run number (UGLY, but not stored in ESD and AOD... )
675   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
676   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
677   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
678   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
679   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
680   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
681   else if (fRun>=136851&&fRun<=139517) {
682     fLHCperiod="LHC10H";
683     fMCperiodTPC="LHC10H8";
684     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
685     fBeamType="PBPB";
686   }
687   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
688   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
689   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
690   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
691   // also for 11e,f use 11d
692   else if (fRun>=160676&&fRun<=162740) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
693   else if (fRun>=162933&&fRun<=165746) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
694   
695   else if (fRun>=166529 && fRun<=170718) {
696     fLHCperiod="LHC11H";
697     fMCperiodTPC="LHC11A10";
698     fBeamType="PBPB";
699     if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
700   }
701   if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
702   // for the moment use LHC12b parameters up to LHC12e
703   if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
704 //   if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
705 //   if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
706 //   if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
707
708 //   if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
709 //   if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
710 //   if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
711 // for the moment use 12g parametrisation for all full gain runs (LHC12f+)
712   if (fRun >= 186636  ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
713
714   //exception new pp MC productions from 2011
715   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
716   // exception for 11f1
717   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
718 }
719
720 //______________________________________________________________________________
721 void AliPIDResponse::SetITSParametrisation()
722 {
723   //
724   // Set the ITS parametrisation
725   //
726 }
727
728 //______________________________________________________________________________
729 void AliPIDResponse::SetTPCPidResponseMaster()
730 {
731   //
732   // Load the TPC pid response functions from the OADB
733   // Load the TPC voltage maps from OADB
734   //
735   //don't load twice for the moment
736    if (fArrPidResponseMaster) return;
737  
738
739   //reset the PID response functions
740   delete fArrPidResponseMaster;
741   fArrPidResponseMaster=NULL;
742   
743   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
744   TFile *f=NULL;
745   if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
746   
747   TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
748   f=TFile::Open(fileNamePIDresponse.Data());
749   if (f && f->IsOpen() && !f->IsZombie()){
750     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
751   }
752   delete f;
753
754   TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
755   f=TFile::Open(fileNameVoltageMaps.Data());
756   if (f && f->IsOpen() && !f->IsZombie()){
757     fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
758   }
759   delete f;
760   
761   if (!fArrPidResponseMaster){
762     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
763     return;
764   }
765   fArrPidResponseMaster->SetOwner();
766
767   if (!fOADBvoltageMaps)
768   {
769     AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
770   }
771   fArrPidResponseMaster->SetOwner();
772 }
773
774 //______________________________________________________________________________
775 void AliPIDResponse::SetTPCParametrisation()
776 {
777   //
778   // Change BB parametrisation for current run
779   //
780   
781   if (fLHCperiod.IsNull()) {
782     AliFatal("No period set, not changing parametrisation");
783     return;
784   }
785   
786   //
787   // Set default parametrisations for data and MC
788   //
789   
790   //data type
791   TString datatype="DATA";
792   //in case of mc fRecoPass is per default 1
793   if (fIsMC) {
794       if(!fTuneMConData) datatype="MC";
795       fRecoPass=1;
796   }
797   
798   //
799   //reset old splines
800   //
801   fTPCResponse.ResetSplines();
802
803   // period
804   TString period=fLHCperiod;
805   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
806
807   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
808   Bool_t found=kFALSE;
809   //
810   //set the new PID splines
811   //
812   if (fArrPidResponseMaster){
813     Int_t recopass = fRecoPass;
814     if(fTuneMConData) recopass = fRecoPassUser;
815     //for MC don't use period information
816     //if (fIsMC) period="[A-Z0-9]*";
817     //for MC use MC period information
818     //pattern for the default entry (valid for all particles)
819     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
820
821     //find particle id ang gain scenario
822     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
823     {
824       TObject *grAll=NULL;
825       TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
826       gainScenario.ToUpper();
827       //loop over entries and filter them
828       for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
829       {
830         TObject *responseFunction=fArrPidResponseMaster->At(iresp);
831         if (responseFunction==NULL) continue;
832         TString responseName=responseFunction->GetName();
833          
834         if (!reg.MatchB(responseName)) continue;
835
836         TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
837         TObject* tmp=NULL;
838         tmp=arr->At(1); if (!tmp) continue;
839         TString particleName=tmp->GetName();
840         tmp=arr->At(3); if (!tmp) continue;
841         TString gainScenarioName=tmp->GetName();
842         delete arr;
843         if (particleName.IsNull()) continue;
844         if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
845         else 
846         {
847           for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
848           {
849             TString particle=AliPID::ParticleName(ispec);
850             particle.ToUpper();
851             //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
852             if ( particle == particleName && gainScenario == gainScenarioName )
853             {
854               fTPCResponse.SetResponseFunction( responseFunction,
855                                                 (AliPID::EParticleType)ispec,
856                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
857               fTPCResponse.SetUseDatabase(kTRUE);
858               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
859               found=kTRUE;
860               // overwrite default with proton spline (for light nuclei)
861               if (ispec==AliPID::kProton) grAll=responseFunction;
862               break;
863             }
864           }
865         }
866       }
867       if (grAll)
868       {
869         for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
870         {
871           if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
872                                                  (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
873           {
874               fTPCResponse.SetResponseFunction( grAll,
875                                                 (AliPID::EParticleType)ispec,
876                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
877               fTPCResponse.SetUseDatabase(kTRUE);
878               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
879               found=kTRUE;
880           }
881         }
882       }
883     }
884   }
885   else AliInfo("no fArrPidResponseMaster");
886
887   if (!found){
888     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
889   }
890
891   //
892   // Setup resolution parametrisation
893   //
894   
895   //default
896   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
897   
898   if (fRun>=122195){
899     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
900   }
901
902   if (fRun>=186636){
903 //   if (fRun>=188356){
904     fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
905   }
906   
907   if (fArrPidResponseMaster)
908   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
909   
910   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
911
912   //read in the voltage map
913   TVectorF* gsm = dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
914   if (gsm) 
915   {
916     fTPCResponse.SetVoltageMap(*gsm);
917     TString vals;
918     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
919     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
920     AliInfo(vals.Data());
921     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
922     AliInfo(vals.Data());
923     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
924     AliInfo(vals.Data());
925     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
926     AliInfo(vals.Data());
927   }
928   else AliInfo("no voltage map, ideal default assumed");
929 }
930
931 //______________________________________________________________________________
932 void AliPIDResponse::SetTRDPidResponseMaster()
933 {
934   //
935   // Load the TRD pid params and references from the OADB
936   //
937   if(fTRDPIDResponseObject) return;
938   AliOADBContainer contParams("contParams"); 
939
940   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
941   if(statusResponse){
942     AliError("Failed initializing PID Response Object from OADB");
943   } else {
944     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
945     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
946     if(!fTRDPIDResponseObject){
947       AliError(Form("TRD Response not found in run %d", fRun));
948     }
949   }
950 }
951
952 //______________________________________________________________________________
953 void AliPIDResponse::InitializeTRDResponse(){
954   //
955   // Set PID Params and references to the TRD PID response
956   // 
957     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
958 }
959
960 //______________________________________________________________________________
961 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
962
963     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
964         // backward compatibility for setting with 8 slices
965         TRDslicesForPID[0] = 0;
966         TRDslicesForPID[1] = 7;
967     }
968     else{
969         if(method==AliTRDPIDResponse::kLQ1D){
970             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
971             TRDslicesForPID[1] = 0;
972         }
973         if(method==AliTRDPIDResponse::kLQ2D){
974             TRDslicesForPID[0] = 1;
975             TRDslicesForPID[1] = 7;
976         }
977     }
978     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
979 }
980
981 //______________________________________________________________________________
982 void AliPIDResponse::SetTOFPidResponseMaster()
983 {
984   //
985   // Load the TOF pid params from the OADB
986   //
987
988   if (fTOFPIDParams) delete fTOFPIDParams;
989   fTOFPIDParams=NULL;
990
991   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
992   if (oadbf && oadbf->IsOpen()) {
993     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
994     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
995     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
996     oadbf->Close();
997     delete oadbc;
998   }
999   delete oadbf;
1000
1001   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1002 }
1003
1004 //______________________________________________________________________________
1005 void AliPIDResponse::InitializeTOFResponse(){
1006   //
1007   // Set PID Params to the TOF PID response
1008   //
1009
1010   AliInfo("TOF PID Params loaded from OADB");
1011   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1012   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1013   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1014                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1015   
1016   for (Int_t i=0;i<4;i++) {
1017     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1018   }
1019   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1020
1021   AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1022   Float_t t0Spread[4];
1023   for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1024   AliInfo(Form("  TZERO spreads from data: (A+C)/2 %f A %f C %f (A'-C')/2: %f",t0Spread[0],t0Spread[1],t0Spread[2],t0Spread[3]));
1025   Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1026   Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1027   if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1028     fResT0AC=t0Spread[3];
1029     fResT0A=TMath::Sqrt(a);
1030     fResT0C=TMath::Sqrt(c);
1031   } else {
1032     AliInfo("  TZERO spreads not present or inconsistent, loading default");
1033     fResT0A=75.;
1034     fResT0C=65.;
1035     fResT0AC=55.;
1036   }
1037   AliInfo(Form("  TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1038
1039 }
1040
1041
1042 //______________________________________________________________________________
1043 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1044   //
1045   // Check whether track is identified as electron under a given electron efficiency hypothesis
1046     //
1047
1048   Double_t probs[AliPID::kSPECIES];
1049   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1050
1051   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1052   // Take mean of the TRD momenta in the given tracklets
1053   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1054   Int_t nmomenta = 0;
1055   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1056     if(vtrack->GetTRDmomentum(iPl) > 0.){
1057       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
1058     }
1059   }
1060   p = TMath::Mean(nmomenta, trdmomenta);
1061
1062   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1063 }
1064
1065 //______________________________________________________________________________
1066 void AliPIDResponse::SetEMCALPidResponseMaster()
1067 {
1068   //
1069   // Load the EMCAL pid response functions from the OADB
1070   //
1071   TObjArray* fEMCALPIDParamsRun      = NULL;
1072   TObjArray* fEMCALPIDParamsPass     = NULL;
1073
1074   if(fEMCALPIDParams) return;
1075   AliOADBContainer contParams("contParams"); 
1076
1077   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1078   if(statusPars){
1079     AliError("Failed initializing PID Params from OADB");
1080   } 
1081   else {
1082     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1083
1084     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1085     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1086     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1087
1088     if(!fEMCALPIDParams){
1089       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1090       AliInfo("Will take the standard LHC11d instead ...");
1091
1092       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1093       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1094       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1095
1096       if(!fEMCALPIDParams){
1097         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));     
1098       }
1099     }
1100   }
1101 }
1102
1103 //______________________________________________________________________________
1104 void AliPIDResponse::InitializeEMCALResponse(){
1105   //
1106   // Set PID Params to the EMCAL PID response
1107   // 
1108   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1109
1110 }
1111
1112 //______________________________________________________________________________
1113 void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
1114 {
1115   //
1116   // create detector PID information and setup the transient pointer in the track
1117   //
1118   
1119   // check if detector number is inside accepted range
1120   if (detector == kNdetectors) return;
1121   
1122   // get detector pid
1123   AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1124   if (!detPID) {
1125     detPID=new AliDetectorPID;
1126     (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1127   }
1128   
1129   //check if values exist
1130   if (detPID->HasRawProbabilitiy(detector) && detPID->HasNumberOfSigmas(detector)) return;
1131   
1132   //TODO: which particles to include? See also the loops below...
1133   Double_t values[AliPID::kSPECIESC]={0};
1134
1135   //nsigmas
1136   for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1137     values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
1138   detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC);
1139   
1140   //probabilities
1141   EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1142   detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1143 }
1144
1145 //______________________________________________________________________________
1146 void AliPIDResponse::FillTrackDetectorPID()
1147 {
1148   //
1149   // create detector PID information and setup the transient pointer in the track
1150   //
1151
1152   if (!fCurrentEvent) return;
1153   
1154   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1155     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1156     if (!track) continue;
1157
1158     for (Int_t idet=0; idet<kNdetectors; ++idet){
1159       FillTrackDetectorPID(track, (EDetector)idet);
1160     }
1161   }
1162 }
1163
1164 //______________________________________________________________________________
1165 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1166   //
1167   // Set TOF response function
1168   // Input option for event_time used
1169   //
1170   
1171     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1172     if(t0spread < 10) t0spread = 80;
1173
1174     // T0 from TOF algorithm
1175
1176     Bool_t flagT0TOF=kFALSE;
1177     Bool_t flagT0T0=kFALSE;
1178     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1179     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1180     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1181
1182     // T0-TOF arrays
1183     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1184     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1185     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1186       estimatedT0event[i]=0.0;
1187       estimatedT0resolution[i]=0.0;
1188       startTimeMask[i] = 0;
1189     }
1190
1191     Float_t resT0A=fResT0A;
1192     Float_t resT0C=fResT0C;
1193     Float_t resT0AC=fResT0AC;
1194     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1195         flagT0T0=kTRUE;
1196     }
1197
1198
1199     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1200
1201     if (tofHeader) { // read global info and T0-TOF
1202       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1203       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1204       if(t0spread < 10) t0spread = 80;
1205
1206       flagT0TOF=kTRUE;
1207       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1208         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1209         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1210         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1211       }
1212
1213       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1214       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1215       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1216       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1217         Int_t icurrent = (Int_t)ibin->GetAt(j);
1218         startTime[icurrent]=t0Bin->GetAt(j);
1219         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1220         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1221       }
1222     }
1223
1224     // for cut of 3 sigma on t0 spread
1225     Float_t t0cut = 3 * t0spread;
1226     if(t0cut < 500) t0cut = 500;
1227
1228     if(option == kFILL_T0){ // T0-FILL is used
1229         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1230           estimatedT0event[i]=0.0;
1231           estimatedT0resolution[i]=t0spread;
1232         }
1233         fTOFResponse.SetT0event(estimatedT0event);
1234         fTOFResponse.SetT0resolution(estimatedT0resolution);
1235     }
1236
1237     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1238         if(flagT0TOF){
1239             fTOFResponse.SetT0event(startTime);
1240             fTOFResponse.SetT0resolution(startTimeRes);
1241             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1242               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1243               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1244             }
1245         }
1246         else{
1247             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1248               estimatedT0event[i]=0.0;
1249               estimatedT0resolution[i]=t0spread;
1250               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1251             }
1252             fTOFResponse.SetT0event(estimatedT0event);
1253             fTOFResponse.SetT0resolution(estimatedT0resolution);
1254         }
1255     }
1256     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1257         Float_t t0AC=-10000;
1258         Float_t t0A=-10000;
1259         Float_t t0C=-10000;
1260         if(flagT0T0){
1261             t0AC= vevent->GetT0TOF()[0];
1262             t0A= vevent->GetT0TOF()[1];
1263             t0C= vevent->GetT0TOF()[2];
1264         }
1265
1266         Float_t t0t0Best = 0;
1267         Float_t t0t0BestRes = 9999;
1268         Int_t t0used=0;
1269         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1270             t0t0Best = t0AC;
1271             t0t0BestRes = resT0AC;
1272             t0used=6;
1273         }
1274         else if(TMath::Abs(t0C) < t0cut){
1275             t0t0Best = t0C;
1276             t0t0BestRes = resT0C;
1277             t0used=4;
1278         }
1279         else if(TMath::Abs(t0A) < t0cut){
1280             t0t0Best = t0A;
1281             t0t0BestRes = resT0A;
1282             t0used=2;
1283         }
1284
1285         if(flagT0TOF){ // if T0-TOF info is available
1286             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1287                 if(t0t0BestRes < 999){
1288                   if(startTimeRes[i] < t0spread){
1289                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1290                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1291                     estimatedT0event[i]=t0best / wtot;
1292                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1293                     startTimeMask[i] = t0used+1;
1294                   }
1295                   else {
1296                     estimatedT0event[i]=t0t0Best;
1297                     estimatedT0resolution[i]=t0t0BestRes;
1298                     startTimeMask[i] = t0used;
1299                   }
1300                 }
1301                 else{
1302                   estimatedT0event[i]=startTime[i];
1303                   estimatedT0resolution[i]=startTimeRes[i];
1304                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1305                 }
1306                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1307             }
1308             fTOFResponse.SetT0event(estimatedT0event);
1309             fTOFResponse.SetT0resolution(estimatedT0resolution);
1310         }
1311         else{ // if no T0-TOF info is available
1312             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1313               fTOFResponse.SetT0binMask(i,t0used);
1314               if(t0t0BestRes < 999){
1315                 estimatedT0event[i]=t0t0Best;
1316                 estimatedT0resolution[i]=t0t0BestRes;
1317               }
1318               else{
1319                 estimatedT0event[i]=0.0;
1320                 estimatedT0resolution[i]=t0spread;
1321               }
1322             }
1323             fTOFResponse.SetT0event(estimatedT0event);
1324             fTOFResponse.SetT0resolution(estimatedT0resolution);
1325         }
1326     }
1327
1328     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1329         Float_t t0AC=-10000;
1330         Float_t t0A=-10000;
1331         Float_t t0C=-10000;
1332         if(flagT0T0){
1333             t0AC= vevent->GetT0TOF()[0];
1334             t0A= vevent->GetT0TOF()[1];
1335             t0C= vevent->GetT0TOF()[2];
1336         }
1337
1338         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1339             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1340               estimatedT0event[i]=t0AC;
1341               estimatedT0resolution[i]=resT0AC;
1342               fTOFResponse.SetT0binMask(i,6);
1343             }
1344         }
1345         else if(TMath::Abs(t0C) < t0cut){
1346             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1347               estimatedT0event[i]=t0C;
1348               estimatedT0resolution[i]=resT0C;
1349               fTOFResponse.SetT0binMask(i,4);
1350             }
1351         }
1352         else if(TMath::Abs(t0A) < t0cut){
1353             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1354               estimatedT0event[i]=t0A;
1355               estimatedT0resolution[i]=resT0A;
1356               fTOFResponse.SetT0binMask(i,2);
1357             }
1358         }
1359         else{
1360             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1361               estimatedT0event[i]=0.0;
1362               estimatedT0resolution[i]=t0spread;
1363               fTOFResponse.SetT0binMask(i,0);
1364             }
1365         }
1366         fTOFResponse.SetT0event(estimatedT0event);
1367         fTOFResponse.SetT0resolution(estimatedT0resolution);
1368     }
1369     delete [] startTime;
1370     delete [] startTimeRes;
1371     delete [] startTimeMask;
1372     delete [] estimatedT0event;
1373     delete [] estimatedT0resolution;
1374 }
1375
1376 //______________________________________________________________________________
1377 // private non cached versions of the PID calculation
1378 //
1379
1380
1381 //______________________________________________________________________________
1382 Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
1383 {
1384   //
1385   // NumberOfSigmas for 'detCode'
1386   //
1387   
1388   switch (detCode){
1389     case kITS:   return GetNumberOfSigmasITS(track, type); break;
1390     case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
1391     case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
1392     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
1393     default: return -999.;
1394   }
1395   
1396 }
1397
1398
1399
1400 //______________________________________________________________________________
1401 Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
1402 {
1403   //
1404   // Calculate the number of sigmas in the ITS
1405   //
1406   
1407   AliVTrack *track=(AliVTrack*)vtrack;
1408   
1409   Float_t dEdx=track->GetITSsignal();
1410   if (dEdx<=0) return -999.;
1411   
1412   UChar_t clumap=track->GetITSClusterMap();
1413   Int_t nPointsForPid=0;
1414   for(Int_t i=2; i<6; i++){
1415     if(clumap&(1<<i)) ++nPointsForPid;
1416   }
1417   Float_t mom=track->P();
1418   
1419   //check for ITS standalone tracks
1420   Bool_t isSA=kTRUE;
1421   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1422   
1423   //TODO: in case of the electron, use the SA parametrisation,
1424   //      this needs to be changed if ITS provides a parametrisation
1425   //      for electrons also for ITS+TPC tracks
1426   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
1427 }
1428
1429 //______________________________________________________________________________
1430 Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
1431 {
1432   //
1433   // Calculate the number of sigmas in the TPC
1434   //
1435   
1436   AliVTrack *track=(AliVTrack*)vtrack;
1437   
1438   Double_t mom  = track->GetTPCmomentum();
1439   Double_t sig  = track->GetTPCsignal();
1440   UInt_t   sigN = track->GetTPCsignalN();
1441   
1442   Double_t nSigma = -999.;
1443   if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
1444   
1445   return nSigma;
1446 }
1447
1448 //______________________________________________________________________________
1449 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1450 {
1451   //
1452   // Calculate the number of sigmas in the EMCAL
1453   //
1454   
1455   AliVTrack *track=(AliVTrack*)vtrack;
1456   
1457   AliVCluster *matchedClus = NULL;
1458   
1459   Double_t mom     = -1.;
1460   Double_t pt      = -1.;
1461   Double_t EovP    = -1.;
1462   Double_t fClsE   = -1.;
1463   
1464   Int_t nMatchClus = -1;
1465   Int_t charge     = 0;
1466   
1467   // Track matching
1468   nMatchClus = track->GetEMCALcluster();
1469   if(nMatchClus > -1){
1470     
1471     mom    = track->P();
1472     pt     = track->Pt();
1473     charge = track->Charge();
1474     
1475     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1476     
1477     if(matchedClus){
1478       
1479       // matched cluster is EMCAL
1480       if(matchedClus->IsEMCAL()){
1481         
1482         fClsE       = matchedClus->E();
1483         EovP        = fClsE/mom;
1484         
1485         
1486         // NSigma value really meaningful only for electrons!
1487         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1488       }
1489     }
1490   }
1491   
1492   return -999;
1493   
1494 }
1495
1496
1497 //______________________________________________________________________________
1498 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1499 {
1500   //
1501   // Compute PID response of 'detCode'
1502   //
1503
1504   switch (detCode){
1505     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1506     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1507     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1508     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1509     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1510     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1511     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1512     default: return kDetNoSignal;
1513   }
1514 }
1515
1516 //______________________________________________________________________________
1517 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1518 {
1519   //
1520   // Compute PID response for the ITS
1521   //
1522   
1523   // look for cached value first
1524   // only the non SA tracks are cached
1525   if (track->GetDetectorPID()){
1526     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
1527   }
1528   
1529   // set flat distribution (no decision)
1530   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1531   
1532   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
1533     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
1534   
1535   //check for ITS standalone tracks
1536   Bool_t isSA=kTRUE;
1537   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1538
1539   Double_t mom=track->P();
1540   Double_t dedx=track->GetITSsignal();
1541   Double_t momITS=mom;
1542   UChar_t clumap=track->GetITSClusterMap();
1543   Int_t nPointsForPid=0;
1544   for(Int_t i=2; i<6; i++){
1545     if(clumap&(1<<i)) ++nPointsForPid;
1546   }
1547
1548   if(nPointsForPid<3) { // track not to be used for combined PID purposes
1549     //       track->ResetStatus(AliVTrack::kITSpid);
1550     return kDetNoSignal;
1551   }
1552
1553   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1554   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
1555     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1556     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1557     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1558     //TODO: in case of the electron, use the SA parametrisation,
1559     //      this needs to be changed if ITS provides a parametrisation
1560     //      for electrons also for ITS+TPC tracks
1561     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1562     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1563       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1564     } else {
1565       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1566       mismatch=kFALSE;
1567     }
1568
1569     // Check for particles heavier than (AliPID::kSPECIES - 1)
1570     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
1571
1572   }
1573
1574   if (mismatch){
1575     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
1576     return kDetNoSignal;
1577   }
1578
1579
1580   return kDetPidOk;
1581 }
1582 //______________________________________________________________________________
1583 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1584 {
1585   //
1586   // Compute PID response for the TPC
1587   //
1588   
1589   // set flat distribution (no decision)
1590   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1591   
1592   // check quality of the track
1593   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
1594   
1595   Double_t mom = track->GetTPCmomentum();
1596   
1597   Double_t dedx=track->GetTPCsignal();
1598   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1599   
1600   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1601   
1602   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1603     AliPID::EParticleType type=AliPID::EParticleType(j);
1604     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
1605     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
1606     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1607       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1608     } else {
1609       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1610       mismatch=kFALSE;
1611     }
1612   }
1613   
1614   if (mismatch){
1615     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1616     return kDetNoSignal;
1617   }
1618   
1619   return kDetPidOk;
1620 }
1621 //______________________________________________________________________________
1622 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1623 {
1624   //
1625   // Compute PID response for the
1626   //
1627   
1628   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
1629   
1630   // set flat distribution (no decision)
1631   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1632   
1633   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
1634   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
1635   
1636   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
1637   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1638     AliPID::EParticleType type=AliPID::EParticleType(j);
1639     Double_t nsigmas=GetNumberOfSigmasTOF(track,type) + meanCorrFactor;
1640     
1641     Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
1642     Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
1643     if (TMath::Abs(nsigmas) > (fRange+2)) {
1644       if(nsigmas < fTOFtail)
1645         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
1646       else
1647         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
1648     } else{
1649       if(nsigmas < fTOFtail)
1650         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
1651       else
1652         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
1653     }
1654     
1655     if (TMath::Abs(nsigmas)<5.){
1656       Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
1657       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
1658     }
1659   }
1660   
1661   if (mismatch){
1662     return kDetMismatch;
1663   }
1664   
1665   // TODO: Light nuclei
1666   
1667   return kDetPidOk;
1668 }
1669 //______________________________________________________________________________
1670 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
1671 {
1672   //
1673   // Compute PID response for the
1674   //
1675   
1676   UInt_t TRDslicesForPID[2];
1677   SetTRDSlices(TRDslicesForPID,PIDmethod);
1678   // set flat distribution (no decision)
1679   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1680   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
1681   
1682   Float_t mom[6]={0.};
1683   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
1684   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
1685   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
1686   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
1687     mom[ilayer] = track->GetTRDmomentum(ilayer);
1688     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
1689       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
1690     }
1691   }
1692   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
1693   return kDetPidOk;
1694   
1695 }
1696 //______________________________________________________________________________
1697 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1698 {
1699   //
1700   // Compute PID response for the EMCAL
1701   //
1702   
1703   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1704   
1705   AliVCluster *matchedClus = NULL;
1706   
1707   Double_t mom     = -1.;
1708   Double_t pt      = -1.;
1709   Double_t EovP    = -1.;
1710   Double_t fClsE   = -1.;
1711   
1712   Int_t nMatchClus = -1;
1713   Int_t charge     = 0;
1714   
1715   // Track matching
1716   nMatchClus = track->GetEMCALcluster();
1717   
1718   if(nMatchClus > -1){
1719     
1720     mom    = track->P();
1721     pt     = track->Pt();
1722     charge = track->Charge();
1723     
1724     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1725     
1726     if(matchedClus){
1727       
1728       // matched cluster is EMCAL
1729       if(matchedClus->IsEMCAL()){
1730         
1731         fClsE       = matchedClus->E();
1732         EovP        = fClsE/mom;
1733         
1734         
1735         // compute the probabilities
1736         if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
1737           
1738           // in case everything is OK
1739           return kDetPidOk;
1740         }
1741       }
1742     }
1743   }
1744   
1745   // in all other cases set flat distribution (no decision)
1746   for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
1747   return kDetNoSignal;
1748   
1749 }
1750 //______________________________________________________________________________
1751 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
1752 {
1753   //
1754   // Compute PID response for the PHOS
1755   //
1756   
1757   // set flat distribution (no decision)
1758   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1759   return kDetNoSignal;
1760 }
1761 //______________________________________________________________________________
1762 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1763 {
1764   //
1765   // Compute PID response for the HMPID
1766   //
1767   // set flat distribution (no decision)
1768   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1769   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
1770   
1771   track->GetHMPIDpid(p);
1772   
1773   return kDetPidOk;
1774 }