4660716b63c0bdc0550327a4b8e92f537903ef34
[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 <TH2D.h>
31 #include <TSpline.h>
32 #include <TFile.h>
33 #include <TArrayI.h>
34 #include <TArrayF.h>
35 #include <TLinearFitter.h>
36
37 #include <AliVEvent.h>
38 #include <AliVTrack.h>
39 #include <AliLog.h>
40 #include <AliPID.h>
41 #include <AliOADBContainer.h>
42 #include <AliTRDPIDResponseObject.h>
43 #include <AliTOFPIDParams.h>
44
45 #include "AliPIDResponse.h"
46 #include "AliDetectorPID.h"
47
48 #include "AliCentrality.h"
49
50 ClassImp(AliPIDResponse);
51
52 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
53 TNamed("PIDResponse","PIDResponse"),
54 fITSResponse(isMC),
55 fTPCResponse(),
56 fTRDResponse(),
57 fTOFResponse(),
58 fEMCALResponse(),
59 fRange(5.),
60 fITSPIDmethod(kITSTruncMean),
61 fIsMC(isMC),
62 fCachePID(kTRUE),
63 fOADBPath(),
64 fCustomTPCpidResponse(),
65 fBeamType("PP"),
66 fLHCperiod(),
67 fMCperiodTPC(),
68 fMCperiodUser(),
69 fCurrentFile(),
70 fRecoPass(0),
71 fRecoPassUser(-1),
72 fRun(0),
73 fOldRun(0),
74 fResT0A(75.),
75 fResT0C(65.),
76 fResT0AC(55.),
77 fArrPidResponseMaster(NULL),
78 fResolutionCorrection(NULL),
79 fOADBvoltageMaps(NULL),
80 fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
81 fTRDPIDResponseObject(NULL),
82 fTOFtail(1.1),
83 fTOFPIDParams(NULL),
84 fEMCALPIDParams(NULL),
85 fCurrentEvent(NULL),
86 fCurrCentrality(0.0),
87 fTuneMConData(kFALSE)
88 {
89   //
90   // default ctor
91   //
92   AliLog::SetClassDebugLevel("AliPIDResponse",0);
93   AliLog::SetClassDebugLevel("AliESDpid",0);
94   AliLog::SetClassDebugLevel("AliAODpidUtil",0);
95
96 }
97
98 //______________________________________________________________________________
99 AliPIDResponse::~AliPIDResponse()
100 {
101   //
102   // dtor
103   //
104   delete fArrPidResponseMaster;
105   delete fTRDPIDResponseObject;
106   delete fTOFPIDParams;
107 }
108
109 //______________________________________________________________________________
110 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
111 TNamed(other),
112 fITSResponse(other.fITSResponse),
113 fTPCResponse(other.fTPCResponse),
114 fTRDResponse(other.fTRDResponse),
115 fTOFResponse(other.fTOFResponse),
116 fEMCALResponse(other.fEMCALResponse),
117 fRange(other.fRange),
118 fITSPIDmethod(other.fITSPIDmethod),
119 fIsMC(other.fIsMC),
120 fCachePID(other.fCachePID),
121 fOADBPath(other.fOADBPath),
122 fCustomTPCpidResponse(other.fCustomTPCpidResponse),
123 fBeamType("PP"),
124 fLHCperiod(),
125 fMCperiodTPC(),
126 fMCperiodUser(other.fMCperiodUser),
127 fCurrentFile(),
128 fRecoPass(0),
129 fRecoPassUser(other.fRecoPassUser),
130 fRun(0),
131 fOldRun(0),
132 fResT0A(75.),
133 fResT0C(65.),
134 fResT0AC(55.),
135 fArrPidResponseMaster(NULL),
136 fResolutionCorrection(NULL),
137 fOADBvoltageMaps(NULL),
138 fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
139 fTRDPIDResponseObject(NULL),
140 fTOFtail(1.1),
141 fTOFPIDParams(NULL),
142 fEMCALPIDParams(NULL),
143 fCurrentEvent(NULL),
144 fCurrCentrality(0.0),
145 fTuneMConData(kFALSE)
146 {
147   //
148   // copy ctor
149   //
150 }
151
152 //______________________________________________________________________________
153 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
154 {
155   //
156   // copy ctor
157   //
158   if(this!=&other) {
159     delete fArrPidResponseMaster;
160     TNamed::operator=(other);
161     fITSResponse=other.fITSResponse;
162     fTPCResponse=other.fTPCResponse;
163     fTRDResponse=other.fTRDResponse;
164     fTOFResponse=other.fTOFResponse;
165     fEMCALResponse=other.fEMCALResponse;
166     fRange=other.fRange;
167     fITSPIDmethod=other.fITSPIDmethod;
168     fOADBPath=other.fOADBPath;
169     fCustomTPCpidResponse=other.fCustomTPCpidResponse;
170     fIsMC=other.fIsMC;
171     fCachePID=other.fCachePID;
172     fBeamType="PP";
173     fLHCperiod="";
174     fMCperiodTPC="";
175     fMCperiodUser=other.fMCperiodUser;
176     fCurrentFile="";
177     fRecoPass=0;
178     fRecoPassUser=other.fRecoPassUser;
179     fRun=0;
180     fOldRun=0;
181     fResT0A=75.;
182     fResT0C=65.;
183     fResT0AC=55.;
184     fArrPidResponseMaster=NULL;
185     fResolutionCorrection=NULL;
186     fOADBvoltageMaps=NULL;
187         fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
188     fTRDPIDResponseObject=NULL;
189     fEMCALPIDParams=NULL;
190     fTOFtail=1.1;
191     fTOFPIDParams=NULL;
192     fCurrentEvent=other.fCurrentEvent;
193
194   }
195   return *this;
196 }
197
198 //______________________________________________________________________________
199 Float_t AliPIDResponse::NumberOfSigmas(EDetector detector, const AliVParticle *vtrack, AliPID::EParticleType type) const
200 {
201   //
202   // NumberOfSigmas for 'detCode'
203   //
204   
205   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
206   // look for cached value first
207   const AliDetectorPID *detPID=track->GetDetectorPID();
208   
209   if ( detPID && detPID->HasNumberOfSigmas(detector)){
210     return detPID->GetNumberOfSigmas(detector, type);
211   } else if (fCachePID) {
212     FillTrackDetectorPID(track, detector);
213     detPID=track->GetDetectorPID();
214     return detPID->GetNumberOfSigmas(detector, type);
215   }
216   
217   return GetNumberOfSigmas(detector, track, type);
218 }
219
220 //______________________________________________________________________________
221 AliPIDResponse::EDetPidStatus AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track,
222                                                              AliPID::EParticleType type, Double_t &val) const
223 {
224   //
225   // NumberOfSigmas with detector status as return value
226   //
227   
228   val=NumberOfSigmas(detCode, track, type);
229   return CheckPIDStatus(detCode, (AliVTrack*)track);
230 }
231
232 //______________________________________________________________________________
233 // public buffered versions of the PID calculation
234 //
235
236 //______________________________________________________________________________
237 Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
238 {
239   //
240   // Calculate the number of sigmas in the ITS
241   //
242   
243   return NumberOfSigmas(kITS, vtrack, type);
244 }
245
246 //______________________________________________________________________________
247 Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
248 {
249   //
250   // Calculate the number of sigmas in the TPC
251   //
252   
253   return NumberOfSigmas(kTPC, vtrack, type);
254 }
255
256 //______________________________________________________________________________
257 Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack, 
258                                            AliPID::EParticleType type,
259                                            AliTPCPIDResponse::ETPCdEdxSource dedxSource) const
260 {
261   //get number of sigmas according the selected TPC gain configuration scenario
262   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
263
264 //   return 0.;
265   Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection);
266
267   return nSigma;
268 }
269
270 //______________________________________________________________________________
271 Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
272 {
273   //
274   // Calculate the number of sigmas in the TOF
275   //
276   
277   return NumberOfSigmas(kTOF, vtrack, type);
278 }
279
280 //______________________________________________________________________________
281 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
282 {
283   //
284   // Calculate the number of sigmas in the EMCAL
285   //
286   
287   return NumberOfSigmas(kEMCAL, vtrack, type);
288 }
289
290 //______________________________________________________________________________
291 Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4])  const
292 {
293   //
294   // emcal nsigma with eop and showershape
295   //
296   AliVTrack *track=(AliVTrack*)vtrack;
297   
298   AliVCluster *matchedClus = NULL;
299
300   Double_t mom     = -1.; 
301   Double_t pt      = -1.; 
302   Double_t EovP    = -1.;
303   Double_t fClsE   = -1.;
304
305   // initialize eop and shower shape parameters
306   eop = -1.;
307   for(Int_t i = 0; i < 4; i++){
308     showershape[i] = -1.;
309   }
310   
311   Int_t nMatchClus = -1;
312   Int_t charge     = 0;
313   
314   // Track matching
315   nMatchClus = track->GetEMCALcluster();
316   if(nMatchClus > -1){
317
318     mom    = track->P();
319     pt     = track->Pt();
320     charge = track->Charge();
321     
322     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
323     
324     if(matchedClus){
325       
326       // matched cluster is EMCAL
327       if(matchedClus->IsEMCAL()){
328         
329         fClsE       = matchedClus->E();
330         EovP        = fClsE/mom;
331         
332         // fill used EMCAL variables here
333         eop            = EovP; // E/p
334         showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
335         showershape[1] = matchedClus->GetM02(); // long axis
336         showershape[2] = matchedClus->GetM20(); // short axis
337         showershape[3] = matchedClus->GetDispersion(); // dispersion
338
339         // look for cached value first
340         const AliDetectorPID *detPID=track->GetDetectorPID();
341         const EDetector detector=kEMCAL;
342         
343         if ( detPID && detPID->HasNumberOfSigmas(detector)){
344           return detPID->GetNumberOfSigmas(detector, type);
345         } else if (fCachePID) {
346           FillTrackDetectorPID(track, detector);
347           detPID=track->GetDetectorPID();
348           return detPID->GetNumberOfSigmas(detector, type);
349         }
350         
351         // NSigma value really meaningful only for electrons!
352         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
353       }
354     }
355   }
356   return -999;
357 }
358
359 //______________________________________________________________________________
360 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode  detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
361 {
362   // Compute PID response of 'detCode'
363   
364   // find detector code from detector bit mask
365   Int_t detector=-1;
366   for (Int_t idet=0; idet<kNdetectors; ++idet) if ( (detCode&(1<<idet)) ) { detector=idet; break; }
367   if (detector==-1) return kDetNoSignal;
368
369   return ComputePIDProbability((EDetector)detector, track, nSpecies, p);
370 }
371
372 //______________________________________________________________________________
373 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetector detector,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
374 {
375   //
376   // Compute PID response of 'detector'
377   //
378
379   const AliDetectorPID *detPID=track->GetDetectorPID();
380
381   if ( detPID && detPID->HasRawProbability(detector)){
382     return detPID->GetRawProbability(detector, p, nSpecies);
383   } else if (fCachePID) {
384     FillTrackDetectorPID(track, detector);
385     detPID=track->GetDetectorPID();
386     return detPID->GetRawProbability(detector, p, nSpecies);
387   }
388   
389   //if no caching return values calculated from scratch
390   return GetComputePIDProbability(detector, track, nSpecies, p);
391 }
392
393 //______________________________________________________________________________
394 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
395 {
396   // Compute PID response for the ITS
397   return ComputePIDProbability(kITS, track, nSpecies, p);
398 }
399
400 //______________________________________________________________________________
401 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
402 {
403   // Compute PID response for the TPC
404   return ComputePIDProbability(kTPC, track, nSpecies, p);
405 }
406
407 //______________________________________________________________________________
408 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
409 {
410   // Compute PID response for the
411   return ComputePIDProbability(kTOF, track, nSpecies, p);
412 }
413
414 //______________________________________________________________________________
415 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
416 {
417   // Compute PID response for the
418   return ComputePIDProbability(kTRD, track, nSpecies, p);
419 }
420
421 //______________________________________________________________________________
422 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
423 {
424   // Compute PID response for the EMCAL
425   return ComputePIDProbability(kEMCAL, track, nSpecies, p);
426 }
427 //______________________________________________________________________________
428 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
429 {
430   // Compute PID response for the PHOS
431   
432   // set flat distribution (no decision)
433   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
434   return kDetNoSignal;
435 }
436
437 //______________________________________________________________________________
438 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
439 {
440   // Compute PID response for the HMPID
441   return ComputePIDProbability(kHMPID, track, nSpecies, p);
442 }
443
444 //______________________________________________________________________________
445 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
446 {
447   // Compute PID response for the
448   return GetComputeTRDProbability(track, nSpecies, p, PIDmethod);
449 }
450
451 //______________________________________________________________________________
452 AliPIDResponse::EDetPidStatus AliPIDResponse::CheckPIDStatus(EDetector detector, const AliVTrack *track) const
453 {
454   // calculate detector pid status
455   
456   const Int_t iDetCode=(Int_t)detector;
457   if (iDetCode<0||iDetCode>=kNdetectors) return kDetNoSignal;
458   const AliDetectorPID *detPID=track->GetDetectorPID();
459   
460   if ( detPID ){
461     return detPID->GetPIDStatus(detector);
462   } else if (fCachePID) {
463     FillTrackDetectorPID(track, detector);
464     detPID=track->GetDetectorPID();
465     return detPID->GetPIDStatus(detector);
466   }
467   
468   // if not buffered and no buffering is requested
469   return GetPIDStatus(detector, track);
470 }
471
472 //______________________________________________________________________________
473 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
474 {
475   //
476   // Apply settings for the current event
477   //
478   fRecoPass=pass;
479   
480
481   fCurrentEvent=NULL;
482   if (!event) return;
483   fCurrentEvent=event;
484   if (run>0) fRun=run;
485   else fRun=event->GetRunNumber();
486   
487   if (fRun!=fOldRun){
488     ExecNewRun();
489     fOldRun=fRun;
490   }
491   
492   //TPC resolution parametrisation PbPb
493   if ( fResolutionCorrection ){
494     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
495     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
496   }
497   
498   //TOF resolution
499   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
500
501
502   // Get and set centrality
503   AliCentrality *centrality = event->GetCentrality();
504   if(centrality){
505     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
506   }
507   else{
508     fCurrCentrality = -1;
509   }
510 }
511
512 //______________________________________________________________________________
513 void AliPIDResponse::ExecNewRun()
514 {
515   //
516   // Things to Execute upon a new run
517   //
518   SetRecoInfo();
519   
520   SetITSParametrisation();
521   
522   SetTPCPidResponseMaster();
523   SetTPCParametrisation();
524   SetTPCEtaMaps();
525
526   SetTRDPidResponseMaster(); 
527   InitializeTRDResponse();
528
529   SetEMCALPidResponseMaster(); 
530   InitializeEMCALResponse();
531   
532   SetTOFPidResponseMaster();
533   InitializeTOFResponse();
534
535   if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
536 }
537
538 //______________________________________________________________________________
539 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
540 {
541   //
542   // Get TPC multiplicity in bins of 150
543   //
544   
545   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
546   Double_t tpcMulti=0.;
547   if(vertexTPC){
548     Double_t vertexContribTPC=vertexTPC->GetNContributors();
549     tpcMulti=vertexContribTPC/150.;
550     if (tpcMulti>20.) tpcMulti=20.;
551   }
552   
553   return tpcMulti;
554 }
555
556 //______________________________________________________________________________
557 void AliPIDResponse::SetRecoInfo()
558 {
559   //
560   // Set reconstruction information
561   //
562   
563   //reset information
564   fLHCperiod="";
565   fMCperiodTPC="";
566   
567   fBeamType="";
568     
569   fBeamType="PP";
570   
571   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
572   TPRegexp reg12a17("LHC1[2-3][a-z]");
573
574   //find the period by run number (UGLY, but not stored in ESD and AOD... )
575   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
576   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
577   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
578   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
579   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
580   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
581   else if (fRun>=136851&&fRun<=139846) {
582     fLHCperiod="LHC10H";
583     fMCperiodTPC="LHC10H8";
584     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
585     fBeamType="PBPB";
586   }
587   else if (fRun>=139847&&fRun<=146974) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
588   //TODO: periods 11B (146975-150721), 11C (150722-155837) are not yet treated assume 11d for the moment
589   else if (fRun>=146975&&fRun<=155837) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
590   else if (fRun>=155838&&fRun<=159649) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
591   // also for 11e (159650-162750),f(162751-165771) use 11d
592   else if (fRun>=159650&&fRun<=162750) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
593   else if (fRun>=162751&&fRun<=165771) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
594   
595   else if (fRun>=165772 && fRun<=170718) {
596     fLHCperiod="LHC11H";
597     fMCperiodTPC="LHC11A10";
598     fBeamType="PBPB";
599     if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
600   }
601   if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
602   // for the moment use LHC12b parameters up to LHC12e
603   if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
604 //   if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
605 //   if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
606 //   if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
607
608 //   if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
609 //   if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
610 //   if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
611 // for the moment use 12g parametrisation for all full gain runs (LHC12f+)
612   if (fRun >= 186636  ) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
613
614   //exception new pp MC productions from 2011
615   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
616   // exception for 11f1
617   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
618   // exception for 12f1a, 12f1b and 12i3
619   if (fCurrentFile.Contains("LHC12f1a/") || fCurrentFile.Contains("LHC12f1b/")
620       || fCurrentFile.Contains("LHC12i3/")) fMCperiodTPC="LHC12F1";
621   // exception for 12c4
622   if (fCurrentFile.Contains("LHC12c4/")) fMCperiodTPC="LHC12C4";
623 }
624
625 //______________________________________________________________________________
626 void AliPIDResponse::SetITSParametrisation()
627 {
628   //
629   // Set the ITS parametrisation
630   //
631 }
632
633  
634 //______________________________________________________________________________
635 void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
636 {
637   if (h->GetBinContent(binX, binY) <= 1e-4)
638     return; // Reject bins without content (within some numerical precision) or with strange content
639     
640   Double_t coord[2] = {0, 0};
641   coord[0] = h->GetXaxis()->GetBinCenter(binX);
642   coord[1] = h->GetYaxis()->GetBinCenter(binY);
643   Double_t binError = h->GetBinError(binX, binY);
644   if (binError <= 0) {
645     binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
646     printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
647   }
648   linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
649 }
650
651
652 //______________________________________________________________________________
653 TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
654 {
655   if (!h)
656     return 0x0;
657   
658   // Interpolate to finer map
659   TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
660   
661   Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
662   Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
663   Int_t nBinsX = 30;
664   // Binning was find to yield good results, if 40 bins are chosen for the range 0.0016 to 0.02. For the new variable range,
665   // scale the number of bins correspondingly
666   Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - 0.0016) * 40);
667   Int_t nBinsXrefined = nBinsX * refineFactorX;
668   Int_t nBinsYrefined = nBinsY * refineFactorY; 
669   
670   TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()),  Form("%s (refined)", h->GetTitle()),
671                             nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
672                             nBinsYrefined, lowerMapBoundY, upperMapBoundY);
673   
674   for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
675     for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
676       
677       hRefined->SetBinContent(binX, binY, 1); // Default value is 1
678       
679       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
680       Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
681       
682       /*OLD
683       linExtrapolation->ClearPoints();
684       
685       // For interpolation: Just take the corresponding bin from the old histo.
686       // For extrapolation: take the last available bin from the old histo.
687       // If the boundaries are to be skipped, also skip the corresponding bins
688       Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
689       if (oldBinX < 1)  
690         oldBinX = 1;
691       if (oldBinX > nBinsX)
692         oldBinX = nBinsX;
693       
694       Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
695       if (oldBinY < 1)  
696         oldBinY = 1;
697       if (oldBinY > nBinsY)
698         oldBinY = nBinsY;
699       
700       // Neighbours left column
701       if (oldBinX >= 2) {
702         if (oldBinY >= 2) {
703           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
704         }
705         
706         AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
707         
708         if (oldBinY < nBinsY) {
709           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
710         }
711       }
712       
713       // Neighbours (and point itself) same column
714       if (oldBinY >= 2) {
715         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
716       }
717         
718       AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
719         
720       if (oldBinY < nBinsY) {
721         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
722       }
723       
724       // Neighbours right column
725       if (oldBinX < nBinsX) {
726         if (oldBinY >= 2) {
727           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
728         }
729         
730         AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
731         
732         if (oldBinY < nBinsY) {
733           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
734         }
735       }
736       
737       
738       // Fit 2D-hyperplane
739       if (linExtrapolation->GetNpoints() <= 0)
740         continue;
741         
742       if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
743         continue;
744       
745       // Fill the bin of the refined histogram with the extrapolated value
746       Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
747                                  + linExtrapolation->GetParameter(2) * centerY;
748       */
749       Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;
750       hRefined->SetBinContent(binX, binY, interpolatedValue);      
751     }
752   } 
753   
754   
755   // Problem: Interpolation does not work before/beyond center of first/last bin (as the name suggests).
756   // Therefore, for each row in dEdx: Take last bin from old map and interpolate values from center and edge.
757   // Assume line through these points and extropolate to last bin of refined map
758   const Double_t firstOldXbinUpEdge = h->GetXaxis()->GetBinUpEdge(1);
759   const Double_t firstOldXbinCenter = h->GetXaxis()->GetBinCenter(1);
760   
761   const Double_t oldXbinHalfWidth = firstOldXbinUpEdge - firstOldXbinCenter;
762   
763   const Double_t lastOldXbinLowEdge = h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
764   const Double_t lastOldXbinCenter = h->GetXaxis()->GetBinCenter(h->GetNbinsX());
765   
766   for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
767     Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
768     
769     const Double_t interpolatedCenterFirstXbin = h->Interpolate(firstOldXbinCenter, centerY);
770     const Double_t interpolatedUpEdgeFirstXbin = h->Interpolate(firstOldXbinUpEdge, centerY);
771     
772     const Double_t extrapolationSlopeFirstXbin = (interpolatedUpEdgeFirstXbin - interpolatedCenterFirstXbin) / oldXbinHalfWidth;
773     const Double_t extrapolationOffsetFirstXbin = interpolatedCenterFirstXbin;
774     
775     
776     const Double_t interpolatedCenterLastXbin = h->Interpolate(lastOldXbinCenter, centerY);
777     const Double_t interpolatedLowEdgeLastXbin = h->Interpolate(lastOldXbinLowEdge, centerY);
778     
779     const Double_t extrapolationSlopeLastXbin = (interpolatedCenterLastXbin - interpolatedLowEdgeLastXbin) / oldXbinHalfWidth;
780     const Double_t extrapolationOffsetLastXbin = interpolatedCenterLastXbin;
781
782     for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
783       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
784      
785       if (centerX < firstOldXbinCenter) {
786         Double_t extrapolatedValue = extrapolationOffsetFirstXbin + (centerX - firstOldXbinCenter) * extrapolationSlopeFirstXbin;
787         hRefined->SetBinContent(binX, binY, extrapolatedValue);      
788       }
789       else if (centerX <= lastOldXbinCenter) {
790         continue;
791       }
792       else {
793         Double_t extrapolatedValue = extrapolationOffsetLastXbin + (centerX - lastOldXbinCenter) * extrapolationSlopeLastXbin;
794         hRefined->SetBinContent(binX, binY, extrapolatedValue);     
795       }
796     }
797   } 
798   
799   delete linExtrapolation;
800   
801   return hRefined;
802 }
803
804 //______________________________________________________________________________
805 void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
806                                    Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
807 {
808   //
809   // Load the TPC eta correction maps from the OADB
810   //
811   
812   if (fUseTPCEtaCorrection == kFALSE) {
813     // Disable eta correction via setting no maps
814     if (!fTPCResponse.SetEtaCorrMap(0x0))
815       AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled"); 
816     else
817       AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
818     
819     if (!fTPCResponse.SetSigmaParams(0x0, 0))
820       AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma"); 
821     else
822       AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
823     
824     return;
825   }
826   
827   TString dataType = "DATA";
828   TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
829   
830   if (fIsMC)  {
831     if (!fTuneMConData) {
832       period=fMCperiodTPC;
833       dataType="MC";
834     }
835     fRecoPass = 1;
836     
837     if (!fTuneMConData && fMCperiodTPC.IsNull()) {
838       AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
839       return;
840     }
841   }
842
843   Int_t recopass = fRecoPass;
844   if (fTuneMConData)
845     recopass = fRecoPassUser;
846   
847   TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), recopass);
848   
849   AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), recopass));
850   
851   // Invalidate old maps
852   fTPCResponse.SetEtaCorrMap(0x0);
853   fTPCResponse.SetSigmaParams(0x0, 0);
854   
855   // Load the eta correction maps
856   AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass)); 
857   
858   Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
859                                               Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
860   if (statusCont) {
861     AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
862   }
863   else {
864     AliInfo(Form("Loading TPC eta correction map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
865     
866     TH2D* etaMap = 0x0;
867     
868     if (fIsMC && !fTuneMConData) {
869       TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
870       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
871       if (!etaMap) {
872         // Try default object
873         etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(defaultObj.Data()));
874       }
875     }
876     else {
877       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetObject(fRun, defaultObj.Data()));
878     }
879     
880         
881     if (!etaMap) {
882       AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
883     }
884     else {
885       TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
886       
887       if (etaMapRefined) {
888         if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
889           AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
890           fTPCResponse.SetEtaCorrMap(0x0);
891         }
892         else {
893           AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s", 
894                        refineFactorMapX, refineFactorMapY, fOADBPath.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle()));
895         }
896         
897         delete etaMapRefined;
898       }
899       else {
900         AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun));
901       }
902     }
903   }
904   
905   // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
906   AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass)); 
907   
908   statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
909                                              Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
910   if (statusCont) {
911     AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
912   }
913   else {
914     AliInfo(Form("Loading TPC eta sigma map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
915     
916     TObjArray* etaSigmaPars = 0x0;
917     
918     if (fIsMC && !fTuneMConData) {
919       TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
920       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
921       if (!etaSigmaPars) {
922         // Try default object
923         etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(defaultObj.Data()));
924       }
925     }
926     else {
927       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetObject(fRun, defaultObj.Data()));
928     }
929     
930     if (!etaSigmaPars) {
931       AliError(Form("TPC eta sigma parametrisation not found for run %d -> Using old sigma parametrisation!!!", fRun));
932     }
933     else {
934       TH2D* etaSigmaPar1Map = dynamic_cast<TH2D *>(etaSigmaPars->FindObject("sigmaPar1Map"));
935       TNamed* sigmaPar0Info = dynamic_cast<TNamed *>(etaSigmaPars->FindObject("sigmaPar0"));
936       Double_t sigmaPar0 = 0.0;
937       
938       if (sigmaPar0Info) {
939         TString sigmaPar0String = sigmaPar0Info->GetTitle();
940         sigmaPar0 = sigmaPar0String.Atof();
941       }
942       else {
943         // Something is weired because the object for parameter 0 could not be loaded -> New sigma parametrisation can not be used!
944         etaSigmaPar1Map = 0x0;
945       }
946       
947       TH2D* etaSigmaPar1MapRefined = RefineHistoViaLinearInterpolation(etaSigmaPar1Map, refineFactorSigmaMapX, refineFactorSigmaMapY);
948       
949       
950       if (etaSigmaPar1MapRefined) {
951         if (!fTPCResponse.SetSigmaParams(etaSigmaPar1MapRefined, sigmaPar0)) {
952           AliError(Form("Failed to set TPC eta sigma map for run %d -> Using old sigma parametrisation!!!", fRun));
953           fTPCResponse.SetSigmaParams(0x0, 0);
954         }
955         else {
956           AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s", 
957                        refineFactorSigmaMapX, refineFactorSigmaMapY, fOADBPath.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle()));
958         }
959         
960         delete etaSigmaPar1MapRefined;
961       }
962       else {
963         AliError(Form("Failed to set TPC eta sigma map for run %d (map was loaded, but couldn't be refined) -> Using old sigma parametrisation!!!",
964                       fRun));
965       }
966     }
967   }
968 }
969
970 //______________________________________________________________________________
971 void AliPIDResponse::SetTPCPidResponseMaster()
972 {
973   //
974   // Load the TPC pid response functions from the OADB
975   // Load the TPC voltage maps from OADB
976   //
977   //don't load twice for the moment
978    if (fArrPidResponseMaster) return;
979  
980
981   //reset the PID response functions
982   delete fArrPidResponseMaster;
983   fArrPidResponseMaster=NULL;
984   
985   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
986   TFile *f=NULL;
987   if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
988   
989   TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
990   f=TFile::Open(fileNamePIDresponse.Data());
991   if (f && f->IsOpen() && !f->IsZombie()){
992     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
993   }
994   delete f;
995
996   TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
997   f=TFile::Open(fileNameVoltageMaps.Data());
998   if (f && f->IsOpen() && !f->IsZombie()){
999     fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
1000   }
1001   delete f;
1002   
1003   if (!fArrPidResponseMaster){
1004     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
1005     return;
1006   }
1007   fArrPidResponseMaster->SetOwner();
1008
1009   if (!fOADBvoltageMaps)
1010   {
1011     AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
1012   }
1013   fArrPidResponseMaster->SetOwner();
1014 }
1015
1016 //______________________________________________________________________________
1017 void AliPIDResponse::SetTPCParametrisation()
1018 {
1019   //
1020   // Change BB parametrisation for current run
1021   //
1022   
1023   //
1024   //reset old splines
1025   //
1026   fTPCResponse.ResetSplines();
1027   
1028   if (fLHCperiod.IsNull()) {
1029     AliError("No period set, not changing parametrisation");
1030     return;
1031   }
1032   
1033   //
1034   // Set default parametrisations for data and MC
1035   //
1036   
1037   //data type
1038   TString datatype="DATA";
1039   //in case of mc fRecoPass is per default 1
1040   if (fIsMC) {
1041       if(!fTuneMConData) datatype="MC";
1042       fRecoPass=1;
1043   }
1044
1045   // period
1046   TString period=fLHCperiod;
1047   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
1048
1049   Int_t recopass = fRecoPass;
1050   if(fTuneMConData) recopass = fRecoPassUser;
1051     
1052   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1053   Bool_t found=kFALSE;
1054   //
1055   //set the new PID splines
1056   //
1057   if (fArrPidResponseMaster){
1058     //for MC don't use period information
1059     //if (fIsMC) period="[A-Z0-9]*";
1060     //for MC use MC period information
1061     //pattern for the default entry (valid for all particles)
1062     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1063
1064     //find particle id and gain scenario
1065     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
1066     {
1067       TObject *grAll=NULL;
1068       TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
1069       gainScenario.ToUpper();
1070       //loop over entries and filter them
1071       for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
1072       {
1073         TObject *responseFunction=fArrPidResponseMaster->At(iresp);
1074         if (responseFunction==NULL) continue;
1075         TString responseName=responseFunction->GetName();
1076          
1077         if (!reg.MatchB(responseName)) continue;
1078
1079         TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
1080         TObject* tmp=NULL;
1081         tmp=arr->At(1); if (!tmp) continue;
1082         TString particleName=tmp->GetName();
1083         tmp=arr->At(3); if (!tmp) continue;
1084         TString gainScenarioName=tmp->GetName();
1085         delete arr;
1086         if (particleName.IsNull()) continue;
1087         if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
1088         else 
1089         {
1090           for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1091           {
1092             TString particle=AliPID::ParticleName(ispec);
1093             particle.ToUpper();
1094             //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
1095             if ( particle == particleName && gainScenario == gainScenarioName )
1096             {
1097               fTPCResponse.SetResponseFunction( responseFunction,
1098                                                 (AliPID::EParticleType)ispec,
1099                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1100               fTPCResponse.SetUseDatabase(kTRUE);
1101               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
1102               found=kTRUE;
1103               break;
1104             }
1105           }
1106         }
1107       }
1108       
1109       // Retrieve responsefunction for pions - will (if available) be used for muons if there are no dedicated muon splines.
1110       // For light nuclei, try to set the proton spline, if no dedicated splines are available.
1111       // In both cases: Use default splines, if no dedicated splines and no pion/proton splines are available.
1112       TObject* responseFunctionPion = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kPion,                             
1113                                                                         (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1114       TObject* responseFunctionProton = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kProton,                             
1115                                                                           (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1116       
1117       for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1118       {
1119         if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
1120           (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
1121         {
1122           if (ispec == AliPID::kMuon) { // Muons
1123             if (responseFunctionPion) {
1124               fTPCResponse.SetResponseFunction( responseFunctionPion,
1125                                                 (AliPID::EParticleType)ispec,
1126                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1127               fTPCResponse.SetUseDatabase(kTRUE);
1128               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionPion->GetName()));
1129               found=kTRUE;  
1130             }
1131             else if (grAll) {
1132               fTPCResponse.SetResponseFunction( grAll,
1133                                                 (AliPID::EParticleType)ispec,
1134                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1135               fTPCResponse.SetUseDatabase(kTRUE);
1136               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1137               found=kTRUE;
1138             }
1139             //else
1140             //  AliError(Form("No splines found for muons (also no pion splines and no default splines) for gain scenario %d!", igainScenario));
1141           }
1142           else if (ispec >= AliPID::kSPECIES) { // Light nuclei
1143             if (responseFunctionProton) {
1144               fTPCResponse.SetResponseFunction( responseFunctionProton,
1145                                                 (AliPID::EParticleType)ispec,
1146                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1147               fTPCResponse.SetUseDatabase(kTRUE);
1148               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionProton->GetName()));
1149               found=kTRUE;  
1150             }
1151             else if (grAll) {
1152               fTPCResponse.SetResponseFunction( grAll,
1153                                                 (AliPID::EParticleType)ispec,
1154                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1155               fTPCResponse.SetUseDatabase(kTRUE);
1156               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1157               found=kTRUE;
1158             }
1159             //else
1160             //  AliError(Form("No splines found for species %d (also no proton splines and no default splines) for gain scenario %d!",
1161             //                ispec, igainScenario));
1162           }
1163         }
1164       }
1165     }
1166   }
1167   else AliInfo("no fArrPidResponseMaster");
1168
1169   if (!found){
1170     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1171   }
1172
1173   //
1174   // Setup resolution parametrisation
1175   //
1176   
1177   //default
1178   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1179   
1180   if (fRun>=122195){
1181     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1182   }
1183
1184   if (fRun>=186636){
1185 //   if (fRun>=188356){
1186     fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1187   }
1188   
1189   if (fArrPidResponseMaster)
1190   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1191   
1192   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
1193
1194   //read in the voltage map
1195   TVectorF* gsm = 0x0;
1196   if (fOADBvoltageMaps) gsm=dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
1197   if (gsm) 
1198   {
1199     fTPCResponse.SetVoltageMap(*gsm);
1200     TString vals;
1201     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1202     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1203     AliInfo(vals.Data());
1204     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1205     AliInfo(vals.Data());
1206     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1207     AliInfo(vals.Data());
1208     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1209     AliInfo(vals.Data());
1210   }
1211   else AliInfo("no voltage map, ideal default assumed");
1212 }
1213
1214 //______________________________________________________________________________
1215 void AliPIDResponse::SetTRDPidResponseMaster()
1216 {
1217   //
1218   // Load the TRD pid params and references from the OADB
1219   //
1220   if(fTRDPIDResponseObject) return;
1221   AliOADBContainer contParams("contParams"); 
1222
1223   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1224   if(statusResponse){
1225     AliError("Failed initializing PID Response Object from OADB");
1226   } else {
1227     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1228     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1229     if(!fTRDPIDResponseObject){
1230       AliError(Form("TRD Response not found in run %d", fRun));
1231     }
1232   }
1233 }
1234
1235 //______________________________________________________________________________
1236 void AliPIDResponse::InitializeTRDResponse(){
1237   //
1238   // Set PID Params and references to the TRD PID response
1239   // 
1240     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1241 }
1242
1243 //______________________________________________________________________________
1244 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1245
1246     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1247         // backward compatibility for setting with 8 slices
1248         TRDslicesForPID[0] = 0;
1249         TRDslicesForPID[1] = 7;
1250     }
1251     else{
1252         if(method==AliTRDPIDResponse::kLQ1D){
1253             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1254             TRDslicesForPID[1] = 0;
1255         }
1256         if(method==AliTRDPIDResponse::kLQ2D){
1257             TRDslicesForPID[0] = 1;
1258             TRDslicesForPID[1] = 7;
1259         }
1260     }
1261     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
1262 }
1263
1264 //______________________________________________________________________________
1265 void AliPIDResponse::SetTOFPidResponseMaster()
1266 {
1267   //
1268   // Load the TOF pid params from the OADB
1269   //
1270
1271   if (fTOFPIDParams) delete fTOFPIDParams;
1272   fTOFPIDParams=NULL;
1273
1274   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1275   if (oadbf && oadbf->IsOpen()) {
1276     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1277     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1278     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
1279     oadbf->Close();
1280     delete oadbc;
1281   }
1282   delete oadbf;
1283
1284   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1285 }
1286
1287 //______________________________________________________________________________
1288 void AliPIDResponse::InitializeTOFResponse(){
1289   //
1290   // Set PID Params to the TOF PID response
1291   //
1292
1293   AliInfo("TOF PID Params loaded from OADB");
1294   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1295   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1296   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1297                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1298   
1299   for (Int_t i=0;i<4;i++) {
1300     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1301   }
1302   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1303
1304   AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1305   Float_t t0Spread[4];
1306   for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1307   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]));
1308   Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1309   Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1310   if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1311     fResT0AC=t0Spread[3];
1312     fResT0A=TMath::Sqrt(a);
1313     fResT0C=TMath::Sqrt(c);
1314   } else {
1315     AliInfo("  TZERO spreads not present or inconsistent, loading default");
1316     fResT0A=75.;
1317     fResT0C=65.;
1318     fResT0AC=55.;
1319   }
1320   AliInfo(Form("  TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1321
1322 }
1323
1324
1325 //______________________________________________________________________________
1326 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1327   //
1328   // Check whether track is identified as electron under a given electron efficiency hypothesis
1329     //
1330
1331   Double_t probs[AliPID::kSPECIES];
1332   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1333
1334   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1335   // Take mean of the TRD momenta in the given tracklets
1336   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1337   Int_t nmomenta = 0;
1338   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1339     if(vtrack->GetTRDmomentum(iPl) > 0.){
1340       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
1341     }
1342   }
1343   p = TMath::Mean(nmomenta, trdmomenta);
1344
1345   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1346 }
1347
1348 //______________________________________________________________________________
1349 void AliPIDResponse::SetEMCALPidResponseMaster()
1350 {
1351   //
1352   // Load the EMCAL pid response functions from the OADB
1353   //
1354   TObjArray* fEMCALPIDParamsRun      = NULL;
1355   TObjArray* fEMCALPIDParamsPass     = NULL;
1356
1357   if(fEMCALPIDParams) return;
1358   AliOADBContainer contParams("contParams"); 
1359
1360   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1361   if(statusPars){
1362     AliError("Failed initializing PID Params from OADB");
1363   } 
1364   else {
1365     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1366
1367     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1368     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1369     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1370
1371     if(!fEMCALPIDParams){
1372       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1373       AliInfo("Will take the standard LHC11d instead ...");
1374
1375       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1376       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1377       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1378
1379       if(!fEMCALPIDParams){
1380         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));     
1381       }
1382     }
1383   }
1384 }
1385
1386 //______________________________________________________________________________
1387 void AliPIDResponse::InitializeEMCALResponse(){
1388   //
1389   // Set PID Params to the EMCAL PID response
1390   // 
1391   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1392
1393 }
1394
1395 //______________________________________________________________________________
1396 void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
1397 {
1398   //
1399   // create detector PID information and setup the transient pointer in the track
1400   //
1401   
1402   // check if detector number is inside accepted range
1403   if (detector == kNdetectors) return;
1404   
1405   // get detector pid
1406   AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1407   if (!detPID) {
1408     detPID=new AliDetectorPID;
1409     (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1410   }
1411   
1412   //check if values exist
1413   if (detPID->HasRawProbability(detector) && detPID->HasNumberOfSigmas(detector)) return;
1414   
1415   //TODO: which particles to include? See also the loops below...
1416   Double_t values[AliPID::kSPECIESC]={0};
1417
1418   //probabilities
1419   EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1420   detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1421   
1422   //nsigmas
1423   for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1424     values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
1425   // the pid status is the same for probabilities and nSigmas, so it is
1426   // fine to use the one from the probabilities also here
1427   detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC, status);
1428   
1429 }
1430
1431 //______________________________________________________________________________
1432 void AliPIDResponse::FillTrackDetectorPID()
1433 {
1434   //
1435   // create detector PID information and setup the transient pointer in the track
1436   //
1437
1438   if (!fCurrentEvent) return;
1439   
1440   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1441     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1442     if (!track) continue;
1443
1444     for (Int_t idet=0; idet<kNdetectors; ++idet){
1445       FillTrackDetectorPID(track, (EDetector)idet);
1446     }
1447   }
1448 }
1449
1450 //______________________________________________________________________________
1451 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1452   //
1453   // Set TOF response function
1454   // Input option for event_time used
1455   //
1456   
1457     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1458     if(t0spread < 10) t0spread = 80;
1459
1460     // T0 from TOF algorithm
1461
1462     Bool_t flagT0TOF=kFALSE;
1463     Bool_t flagT0T0=kFALSE;
1464     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1465     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1466     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1467
1468     // T0-TOF arrays
1469     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1470     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1471     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1472       estimatedT0event[i]=0.0;
1473       estimatedT0resolution[i]=0.0;
1474       startTimeMask[i] = 0;
1475     }
1476
1477     Float_t resT0A=fResT0A;
1478     Float_t resT0C=fResT0C;
1479     Float_t resT0AC=fResT0AC;
1480     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1481         flagT0T0=kTRUE;
1482     }
1483
1484
1485     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1486
1487     if (tofHeader) { // read global info and T0-TOF
1488       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1489       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1490       if(t0spread < 10) t0spread = 80;
1491
1492       flagT0TOF=kTRUE;
1493       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1494         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1495         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1496         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1497       }
1498
1499       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1500       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1501       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1502       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1503         Int_t icurrent = (Int_t)ibin->GetAt(j);
1504         startTime[icurrent]=t0Bin->GetAt(j);
1505         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1506         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1507       }
1508     }
1509
1510     // for cut of 3 sigma on t0 spread
1511     Float_t t0cut = 3 * t0spread;
1512     if(t0cut < 500) t0cut = 500;
1513
1514     if(option == kFILL_T0){ // T0-FILL is used
1515         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1516           estimatedT0event[i]=0.0;
1517           estimatedT0resolution[i]=t0spread;
1518         }
1519         fTOFResponse.SetT0event(estimatedT0event);
1520         fTOFResponse.SetT0resolution(estimatedT0resolution);
1521     }
1522
1523     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1524         if(flagT0TOF){
1525             fTOFResponse.SetT0event(startTime);
1526             fTOFResponse.SetT0resolution(startTimeRes);
1527             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1528               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1529               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1530             }
1531         }
1532         else{
1533             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1534               estimatedT0event[i]=0.0;
1535               estimatedT0resolution[i]=t0spread;
1536               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1537             }
1538             fTOFResponse.SetT0event(estimatedT0event);
1539             fTOFResponse.SetT0resolution(estimatedT0resolution);
1540         }
1541     }
1542     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1543         Float_t t0AC=-10000;
1544         Float_t t0A=-10000;
1545         Float_t t0C=-10000;
1546         if(flagT0T0){
1547             t0A= vevent->GetT0TOF()[1];
1548             t0C= vevent->GetT0TOF()[2];
1549         //      t0AC= vevent->GetT0TOF()[0];
1550         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1551         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1552         t0AC /= resT0AC*resT0AC;
1553         }
1554
1555         Float_t t0t0Best = 0;
1556         Float_t t0t0BestRes = 9999;
1557         Int_t t0used=0;
1558         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1559             t0t0Best = t0AC;
1560             t0t0BestRes = resT0AC;
1561             t0used=6;
1562         }
1563         else if(TMath::Abs(t0C) < t0cut){
1564             t0t0Best = t0C;
1565             t0t0BestRes = resT0C;
1566             t0used=4;
1567         }
1568         else if(TMath::Abs(t0A) < t0cut){
1569             t0t0Best = t0A;
1570             t0t0BestRes = resT0A;
1571             t0used=2;
1572         }
1573
1574         if(flagT0TOF){ // if T0-TOF info is available
1575             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1576                 if(t0t0BestRes < 999){
1577                   if(startTimeRes[i] < t0spread){
1578                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1579                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1580                     estimatedT0event[i]=t0best / wtot;
1581                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1582                     startTimeMask[i] = t0used+1;
1583                   }
1584                   else {
1585                     estimatedT0event[i]=t0t0Best;
1586                     estimatedT0resolution[i]=t0t0BestRes;
1587                     startTimeMask[i] = t0used;
1588                   }
1589                 }
1590                 else{
1591                   estimatedT0event[i]=startTime[i];
1592                   estimatedT0resolution[i]=startTimeRes[i];
1593                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1594                 }
1595                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1596             }
1597             fTOFResponse.SetT0event(estimatedT0event);
1598             fTOFResponse.SetT0resolution(estimatedT0resolution);
1599         }
1600         else{ // if no T0-TOF info is available
1601             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1602               fTOFResponse.SetT0binMask(i,t0used);
1603               if(t0t0BestRes < 999){
1604                 estimatedT0event[i]=t0t0Best;
1605                 estimatedT0resolution[i]=t0t0BestRes;
1606               }
1607               else{
1608                 estimatedT0event[i]=0.0;
1609                 estimatedT0resolution[i]=t0spread;
1610               }
1611             }
1612             fTOFResponse.SetT0event(estimatedT0event);
1613             fTOFResponse.SetT0resolution(estimatedT0resolution);
1614         }
1615     }
1616
1617     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1618         Float_t t0AC=-10000;
1619         Float_t t0A=-10000;
1620         Float_t t0C=-10000;
1621         if(flagT0T0){
1622             t0A= vevent->GetT0TOF()[1];
1623             t0C= vevent->GetT0TOF()[2];
1624         //      t0AC= vevent->GetT0TOF()[0];
1625         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1626         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1627         t0AC /= resT0AC*resT0AC;
1628         }
1629
1630         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1631             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1632               estimatedT0event[i]=t0AC;
1633               estimatedT0resolution[i]=resT0AC;
1634               fTOFResponse.SetT0binMask(i,6);
1635             }
1636         }
1637         else if(TMath::Abs(t0C) < t0cut){
1638             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1639               estimatedT0event[i]=t0C;
1640               estimatedT0resolution[i]=resT0C;
1641               fTOFResponse.SetT0binMask(i,4);
1642             }
1643         }
1644         else if(TMath::Abs(t0A) < t0cut){
1645             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1646               estimatedT0event[i]=t0A;
1647               estimatedT0resolution[i]=resT0A;
1648               fTOFResponse.SetT0binMask(i,2);
1649             }
1650         }
1651         else{
1652             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1653               estimatedT0event[i]=0.0;
1654               estimatedT0resolution[i]=t0spread;
1655               fTOFResponse.SetT0binMask(i,0);
1656             }
1657         }
1658         fTOFResponse.SetT0event(estimatedT0event);
1659         fTOFResponse.SetT0resolution(estimatedT0resolution);
1660     }
1661     delete [] startTime;
1662     delete [] startTimeRes;
1663     delete [] startTimeMask;
1664     delete [] estimatedT0event;
1665     delete [] estimatedT0resolution;
1666 }
1667
1668 //______________________________________________________________________________
1669 // private non cached versions of the PID calculation
1670 //
1671
1672
1673 //______________________________________________________________________________
1674 Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detector, const AliVParticle *vtrack, AliPID::EParticleType type) const
1675 {
1676   //
1677   // NumberOfSigmas for 'detCode'
1678   //
1679
1680   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
1681   
1682   switch (detector){
1683     case kITS:   return GetNumberOfSigmasITS(track, type); break;
1684     case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
1685     case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
1686     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
1687     default: return -999.;
1688   }
1689
1690   return -999.;
1691 }
1692
1693 //______________________________________________________________________________
1694 Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
1695 {
1696   //
1697   // Calculate the number of sigmas in the ITS
1698   //
1699   
1700   AliVTrack *track=(AliVTrack*)vtrack;
1701
1702   const EDetPidStatus pidStatus=GetITSPIDStatus(track);
1703   if (pidStatus!=kDetPidOk) return -999.;
1704     
1705   UChar_t clumap=track->GetITSClusterMap();
1706   Int_t nPointsForPid=0;
1707   for(Int_t i=2; i<6; i++){
1708     if(clumap&(1<<i)) ++nPointsForPid;
1709   }
1710   Float_t mom=track->P();
1711   
1712   //check for ITS standalone tracks
1713   Bool_t isSA=kTRUE;
1714   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1715   
1716   const Float_t dEdx=track->GetITSsignal();
1717
1718   //TODO: in case of the electron, use the SA parametrisation,
1719   //      this needs to be changed if ITS provides a parametrisation
1720   //      for electrons also for ITS+TPC tracks
1721   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
1722 }
1723
1724 //______________________________________________________________________________
1725 Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
1726 {
1727   //
1728   // Calculate the number of sigmas in the TPC
1729   //
1730   
1731   AliVTrack *track=(AliVTrack*)vtrack;
1732
1733   const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
1734   if (pidStatus!=kDetPidOk) return -999.;
1735   
1736   Double_t nSigma = -999.;
1737   
1738   if (fTuneMConData)
1739     this->GetTPCsignalTunedOnData(track);
1740   
1741   nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1742   
1743   return nSigma;
1744 }
1745
1746 //______________________________________________________________________________
1747 Float_t AliPIDResponse::GetNumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
1748 {
1749   //
1750   // Calculate the number of sigmas in the TOF
1751   //
1752   
1753   AliVTrack *track=(AliVTrack*)vtrack;
1754
1755   const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
1756   if (pidStatus!=kDetPidOk) return -999.;
1757
1758   
1759   return GetNumberOfSigmasTOFold(vtrack, type);
1760 }
1761
1762 //______________________________________________________________________________
1763 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1764 {
1765   //
1766   // Calculate the number of sigmas in the EMCAL
1767   //
1768   
1769   AliVTrack *track=(AliVTrack*)vtrack;
1770
1771   const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
1772   if (pidStatus!=kDetPidOk) return -999.;
1773
1774   const Int_t nMatchClus = track->GetEMCALcluster();
1775   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1776   
1777   const Double_t mom    = track->P();
1778   const Double_t pt     = track->Pt();
1779   const Int_t    charge = track->Charge();
1780   const Double_t fClsE  = matchedClus->E();
1781   const Double_t EovP   = fClsE/mom;
1782   
1783   return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1784 }
1785
1786
1787 //______________________________________________________________________________
1788 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1789 {
1790   //
1791   // Compute PID response of 'detCode'
1792   //
1793
1794   switch (detCode){
1795     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1796     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1797     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1798     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1799     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1800     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1801     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1802     default: return kDetNoSignal;
1803   }
1804 }
1805
1806 //______________________________________________________________________________
1807 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1808 {
1809   //
1810   // Compute PID response for the ITS
1811   //
1812   
1813   // set flat distribution (no decision)
1814   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1815   
1816   const EDetPidStatus pidStatus=GetITSPIDStatus(track);
1817   if (pidStatus!=kDetPidOk) return pidStatus;
1818   
1819   if (track->GetDetectorPID()){
1820     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
1821   }
1822   
1823   //check for ITS standalone tracks
1824   Bool_t isSA=kTRUE;
1825   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1826
1827   Double_t mom=track->P();
1828   Double_t dedx=track->GetITSsignal();
1829   Double_t momITS=mom;
1830   UChar_t clumap=track->GetITSClusterMap();
1831   Int_t nPointsForPid=0;
1832   for(Int_t i=2; i<6; i++){
1833     if(clumap&(1<<i)) ++nPointsForPid;
1834   }
1835
1836   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1837   for (Int_t j=0; j<nSpecies; j++) {
1838     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1839     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1840     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1841     //TODO: in case of the electron, use the SA parametrisation,
1842     //      this needs to be changed if ITS provides a parametrisation
1843     //      for electrons also for ITS+TPC tracks
1844     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1845     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1846       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1847     } else {
1848       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1849       mismatch=kFALSE;
1850     }
1851   }
1852
1853   if (mismatch){
1854     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1855   }
1856
1857   return kDetPidOk;
1858 }
1859 //______________________________________________________________________________
1860 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1861 {
1862   //
1863   // Compute PID response for the TPC
1864   //
1865   
1866   // set flat distribution (no decision)
1867   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1868   
1869   const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
1870   if (pidStatus!=kDetPidOk) return pidStatus;
1871   
1872   Double_t dedx=track->GetTPCsignal();
1873   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1874   
1875   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1876   
1877   Double_t bethe = 0.;
1878   Double_t sigma = 0.;
1879   
1880   for (Int_t j=0; j<nSpecies; j++) {
1881     AliPID::EParticleType type=AliPID::EParticleType(j);
1882     
1883     bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1884     sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1885     
1886     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1887       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1888     } else {
1889       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1890       mismatch=kFALSE;
1891     }
1892   }
1893   
1894   if (mismatch){
1895     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1896   }
1897   
1898   return kDetPidOk;
1899 }
1900 //______________________________________________________________________________
1901 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1902 {
1903   //
1904   // Compute PID probabilities for TOF
1905   //
1906   
1907   // set flat distribution (no decision)
1908   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1909   
1910   const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
1911   if (pidStatus!=kDetPidOk) return pidStatus;
1912
1913   const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
1914   
1915   for (Int_t j=0; j<nSpecies; j++) {
1916     AliPID::EParticleType type=AliPID::EParticleType(j);
1917     const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
1918     
1919     const Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
1920     const Double_t sig     = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
1921     if (TMath::Abs(nsigmas) > (fRange+2)) {
1922       if(nsigmas < fTOFtail)
1923         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
1924       else
1925         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
1926     } else{
1927       if(nsigmas < fTOFtail)
1928         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
1929       else
1930         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
1931     }    
1932   }
1933   
1934   return kDetPidOk;
1935 }
1936 //______________________________________________________________________________
1937 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
1938 {
1939   //
1940   // Compute PID probabilities for the TRD
1941   //
1942   
1943   // set flat distribution (no decision)
1944   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1945   
1946   const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
1947   if (pidStatus!=kDetPidOk) return pidStatus;
1948
1949   UInt_t TRDslicesForPID[2];
1950   SetTRDSlices(TRDslicesForPID,PIDmethod);
1951   
1952   Float_t mom[6]={0.};
1953   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
1954   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
1955   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
1956   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
1957     mom[ilayer] = track->GetTRDmomentum(ilayer);
1958     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
1959       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
1960     }
1961   }
1962   
1963   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
1964   return kDetPidOk;
1965 }
1966
1967 //______________________________________________________________________________
1968 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1969 {
1970   //
1971   // Compute PID response for the EMCAL
1972   //
1973   
1974   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1975
1976   const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
1977   if (pidStatus!=kDetPidOk) return pidStatus;
1978
1979   const Int_t nMatchClus = track->GetEMCALcluster();
1980   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1981   
1982   const Double_t mom    = track->P();
1983   const Double_t pt     = track->Pt();
1984   const Int_t    charge = track->Charge();
1985   const Double_t fClsE  = matchedClus->E();
1986   const Double_t EovP   = fClsE/mom;
1987   
1988   // compute the probabilities
1989   fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p);
1990   return kDetPidOk;
1991 }
1992
1993 //______________________________________________________________________________
1994 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
1995 {
1996   //
1997   // Compute PID response for the PHOS
1998   //
1999   
2000   // set flat distribution (no decision)
2001   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2002   return kDetNoSignal;
2003 }
2004
2005 //______________________________________________________________________________
2006 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2007 {
2008   //
2009   // Compute PID response for the HMPID
2010   //
2011   
2012   // set flat distribution (no decision)
2013   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2014   
2015   const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
2016   if (pidStatus!=kDetPidOk) return pidStatus;
2017   
2018   track->GetHMPIDpid(p);
2019   
2020   return kDetPidOk;
2021 }
2022
2023 //______________________________________________________________________________
2024 AliPIDResponse::EDetPidStatus AliPIDResponse::GetITSPIDStatus(const AliVTrack *track) const
2025 {
2026   // compute ITS pid status
2027
2028   // check status bits
2029   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
2030     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
2031
2032   const Float_t dEdx=track->GetITSsignal();
2033   if (dEdx<=0) return kDetNoSignal;
2034   
2035   // requite at least 3 pid clusters
2036   const UChar_t clumap=track->GetITSClusterMap();
2037   Int_t nPointsForPid=0;
2038   for(Int_t i=2; i<6; i++){
2039     if(clumap&(1<<i)) ++nPointsForPid;
2040   }
2041   
2042   if(nPointsForPid<3) { 
2043     return kDetNoSignal;
2044   }
2045   
2046   return kDetPidOk;
2047 }
2048
2049 //______________________________________________________________________________
2050 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetTPCPIDStatus(const AliVTrack *track) const
2051 {
2052   // compute TPC pid status
2053   
2054   // check quality of the track
2055   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
2056
2057   // check pid values
2058   const Double_t dedx=track->GetTPCsignal();
2059   const UShort_t signalN=track->GetTPCsignalN();
2060   if (signalN<10 || dedx<10) return kDetNoSignal;
2061
2062   if (!(fArrPidResponseMaster && fArrPidResponseMaster->At(AliPID::kPion))) return kDetNoParams;
2063   
2064   return kDetPidOk;
2065 }
2066
2067 //______________________________________________________________________________
2068 AliPIDResponse::EDetPidStatus AliPIDResponse::GetTRDPIDStatus(const AliVTrack *track) const
2069 {
2070   // compute TRD pid status
2071
2072   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
2073   return kDetPidOk;
2074 }
2075
2076 //______________________________________________________________________________
2077 AliPIDResponse::EDetPidStatus AliPIDResponse::GetTOFPIDStatus(const AliVTrack *track) const
2078 {
2079   // compute TOF pid status
2080
2081   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
2082   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
2083
2084   return kDetPidOk;
2085 }
2086
2087 //______________________________________________________________________________
2088 Float_t AliPIDResponse::GetTOFMismatchProbability(const AliVTrack *track) const
2089 {
2090   // compute mismatch probability cross-checking at 5 sigmas with TPC
2091   // currently just implemented as a 5 sigma compatibility cut
2092
2093   // check pid status
2094   const EDetPidStatus tofStatus=GetTOFPIDStatus(track);
2095   if (tofStatus!=kDetPidOk) return 0.;
2096
2097   //mismatch
2098   const EDetPidStatus tpcStatus=GetTPCPIDStatus(track);
2099   if (tpcStatus!=kDetPidOk) return 0.;
2100   
2101   const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2102   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
2103   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
2104     AliPID::EParticleType type=AliPID::EParticleType(j);
2105     const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
2106     
2107     if (TMath::Abs(nsigmas)<5.){
2108       const Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
2109       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
2110     }
2111   }
2112   
2113   if (mismatch){
2114     return 1.;
2115   }
2116   
2117   return 0.;
2118 }
2119
2120
2121
2122 //______________________________________________________________________________
2123 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetHMPIDPIDStatus(const AliVTrack *track) const
2124 {
2125   // compute HMPID pid status
2126   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
2127   return kDetPidOk;
2128 }
2129
2130 //______________________________________________________________________________
2131 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetPHOSPIDStatus(const AliVTrack */*track*/) const
2132 {
2133   // compute PHOS pid status
2134   return kDetNoSignal;  
2135 }
2136
2137 //______________________________________________________________________________
2138 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetEMCALPIDStatus(const AliVTrack *track) const
2139 {
2140   // compute EMCAL pid status
2141
2142
2143   // Track matching
2144   const Int_t nMatchClus = track->GetEMCALcluster();
2145   if (nMatchClus<0) return kDetNoSignal;
2146
2147   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2148
2149   if (!(matchedClus && matchedClus->IsEMCAL())) return kDetNoSignal;
2150
2151   const Int_t charge = track->Charge();
2152   if (TMath::Abs(charge)!=1) return kDetNoSignal;
2153
2154   if (!(fEMCALPIDParams && fEMCALPIDParams->At(AliPID::kElectron))) return kDetNoParams;
2155   
2156   return kDetPidOk;
2157
2158 }
2159
2160 //______________________________________________________________________________
2161 AliPIDResponse::EDetPidStatus AliPIDResponse::GetPIDStatus(EDetector detector, const AliVTrack *track) const
2162 {
2163   //
2164   // check pid status for a track
2165   //
2166
2167   switch (detector){
2168     case kITS:   return GetITSPIDStatus(track);   break;
2169     case kTPC:   return GetTPCPIDStatus(track);   break;
2170     case kTRD:   return GetTRDPIDStatus(track);   break;
2171     case kTOF:   return GetTOFPIDStatus(track);   break;
2172     case kPHOS:  return GetPHOSPIDStatus(track);  break;
2173     case kEMCAL: return GetEMCALPIDStatus(track); break;
2174     case kHMPID: return GetHMPIDPIDStatus(track); break;
2175     default: return kDetNoSignal;
2176   }
2177   return kDetNoSignal;
2178   
2179 }