New splines and first eta maps for the 2013 pPb runs (Benjamin)
[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  && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
613   if (fRun >= 194480) { fLHCperiod="LHC13B"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
614
615   //exception new pp MC productions from 2011
616   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
617   // exception for 11f1
618   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
619   // exception for 12f1a, 12f1b and 12i3
620   if (fCurrentFile.Contains("LHC12f1a/") || fCurrentFile.Contains("LHC12f1b/")
621       || fCurrentFile.Contains("LHC12i3/")) fMCperiodTPC="LHC12F1";
622   // exception for 12c4
623   if (fCurrentFile.Contains("LHC12c4/")) fMCperiodTPC="LHC12C4";
624 }
625
626 //______________________________________________________________________________
627 void AliPIDResponse::SetITSParametrisation()
628 {
629   //
630   // Set the ITS parametrisation
631   //
632 }
633
634  
635 //______________________________________________________________________________
636 void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
637 {
638   if (h->GetBinContent(binX, binY) <= 1e-4)
639     return; // Reject bins without content (within some numerical precision) or with strange content
640     
641   Double_t coord[2] = {0, 0};
642   coord[0] = h->GetXaxis()->GetBinCenter(binX);
643   coord[1] = h->GetYaxis()->GetBinCenter(binY);
644   Double_t binError = h->GetBinError(binX, binY);
645   if (binError <= 0) {
646     binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
647     printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
648   }
649   linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
650 }
651
652
653 //______________________________________________________________________________
654 TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
655 {
656   if (!h)
657     return 0x0;
658   
659   // Interpolate to finer map
660   TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
661   
662   Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
663   Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
664   Int_t nBinsX = 30;
665   // 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,
666   // scale the number of bins correspondingly
667   Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - 0.0016) * 40);
668   Int_t nBinsXrefined = nBinsX * refineFactorX;
669   Int_t nBinsYrefined = nBinsY * refineFactorY; 
670   
671   TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()),  Form("%s (refined)", h->GetTitle()),
672                             nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
673                             nBinsYrefined, lowerMapBoundY, upperMapBoundY);
674   
675   for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
676     for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
677       
678       hRefined->SetBinContent(binX, binY, 1); // Default value is 1
679       
680       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
681       Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
682       
683       /*OLD
684       linExtrapolation->ClearPoints();
685       
686       // For interpolation: Just take the corresponding bin from the old histo.
687       // For extrapolation: take the last available bin from the old histo.
688       // If the boundaries are to be skipped, also skip the corresponding bins
689       Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
690       if (oldBinX < 1)  
691         oldBinX = 1;
692       if (oldBinX > nBinsX)
693         oldBinX = nBinsX;
694       
695       Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
696       if (oldBinY < 1)  
697         oldBinY = 1;
698       if (oldBinY > nBinsY)
699         oldBinY = nBinsY;
700       
701       // Neighbours left column
702       if (oldBinX >= 2) {
703         if (oldBinY >= 2) {
704           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
705         }
706         
707         AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
708         
709         if (oldBinY < nBinsY) {
710           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
711         }
712       }
713       
714       // Neighbours (and point itself) same column
715       if (oldBinY >= 2) {
716         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
717       }
718         
719       AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
720         
721       if (oldBinY < nBinsY) {
722         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
723       }
724       
725       // Neighbours right column
726       if (oldBinX < nBinsX) {
727         if (oldBinY >= 2) {
728           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
729         }
730         
731         AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
732         
733         if (oldBinY < nBinsY) {
734           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
735         }
736       }
737       
738       
739       // Fit 2D-hyperplane
740       if (linExtrapolation->GetNpoints() <= 0)
741         continue;
742         
743       if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
744         continue;
745       
746       // Fill the bin of the refined histogram with the extrapolated value
747       Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
748                                  + linExtrapolation->GetParameter(2) * centerY;
749       */
750       Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;
751       hRefined->SetBinContent(binX, binY, interpolatedValue);      
752     }
753   } 
754   
755   
756   // Problem: Interpolation does not work before/beyond center of first/last bin (as the name suggests).
757   // Therefore, for each row in dEdx: Take last bin from old map and interpolate values from center and edge.
758   // Assume line through these points and extropolate to last bin of refined map
759   const Double_t firstOldXbinUpEdge = h->GetXaxis()->GetBinUpEdge(1);
760   const Double_t firstOldXbinCenter = h->GetXaxis()->GetBinCenter(1);
761   
762   const Double_t oldXbinHalfWidth = firstOldXbinUpEdge - firstOldXbinCenter;
763   
764   const Double_t lastOldXbinLowEdge = h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
765   const Double_t lastOldXbinCenter = h->GetXaxis()->GetBinCenter(h->GetNbinsX());
766   
767   for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
768     Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
769     
770     const Double_t interpolatedCenterFirstXbin = h->Interpolate(firstOldXbinCenter, centerY);
771     const Double_t interpolatedUpEdgeFirstXbin = h->Interpolate(firstOldXbinUpEdge, centerY);
772     
773     const Double_t extrapolationSlopeFirstXbin = (interpolatedUpEdgeFirstXbin - interpolatedCenterFirstXbin) / oldXbinHalfWidth;
774     const Double_t extrapolationOffsetFirstXbin = interpolatedCenterFirstXbin;
775     
776     
777     const Double_t interpolatedCenterLastXbin = h->Interpolate(lastOldXbinCenter, centerY);
778     const Double_t interpolatedLowEdgeLastXbin = h->Interpolate(lastOldXbinLowEdge, centerY);
779     
780     const Double_t extrapolationSlopeLastXbin = (interpolatedCenterLastXbin - interpolatedLowEdgeLastXbin) / oldXbinHalfWidth;
781     const Double_t extrapolationOffsetLastXbin = interpolatedCenterLastXbin;
782
783     for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
784       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
785      
786       if (centerX < firstOldXbinCenter) {
787         Double_t extrapolatedValue = extrapolationOffsetFirstXbin + (centerX - firstOldXbinCenter) * extrapolationSlopeFirstXbin;
788         hRefined->SetBinContent(binX, binY, extrapolatedValue);      
789       }
790       else if (centerX <= lastOldXbinCenter) {
791         continue;
792       }
793       else {
794         Double_t extrapolatedValue = extrapolationOffsetLastXbin + (centerX - lastOldXbinCenter) * extrapolationSlopeLastXbin;
795         hRefined->SetBinContent(binX, binY, extrapolatedValue);     
796       }
797     }
798   } 
799   
800   delete linExtrapolation;
801   
802   return hRefined;
803 }
804
805 //______________________________________________________________________________
806 void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
807                                    Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
808 {
809   //
810   // Load the TPC eta correction maps from the OADB
811   //
812   
813   if (fUseTPCEtaCorrection == kFALSE) {
814     // Disable eta correction via setting no maps
815     if (!fTPCResponse.SetEtaCorrMap(0x0))
816       AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled"); 
817     else
818       AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
819     
820     if (!fTPCResponse.SetSigmaParams(0x0, 0))
821       AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma"); 
822     else
823       AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
824     
825     return;
826   }
827   
828   TString dataType = "DATA";
829   TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
830   
831   if (fIsMC)  {
832     if (!fTuneMConData) {
833       period=fMCperiodTPC;
834       dataType="MC";
835     }
836     fRecoPass = 1;
837     
838     if (!fTuneMConData && fMCperiodTPC.IsNull()) {
839       AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
840       return;
841     }
842   }
843
844   Int_t recopass = fRecoPass;
845   if (fTuneMConData)
846     recopass = fRecoPassUser;
847   
848   TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), recopass);
849   
850   AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), recopass));
851   
852   // Invalidate old maps
853   fTPCResponse.SetEtaCorrMap(0x0);
854   fTPCResponse.SetSigmaParams(0x0, 0);
855   
856   // Load the eta correction maps
857   AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass)); 
858   
859   Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
860                                               Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
861   if (statusCont) {
862     AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
863   }
864   else {
865     AliInfo(Form("Loading TPC eta correction map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
866     
867     TH2D* etaMap = 0x0;
868     
869     if (fIsMC && !fTuneMConData) {
870       TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
871       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
872       if (!etaMap) {
873         // Try default object
874         etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(defaultObj.Data()));
875       }
876     }
877     else {
878       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetObject(fRun, defaultObj.Data()));
879     }
880     
881         
882     if (!etaMap) {
883       AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
884     }
885     else {
886       TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
887       
888       if (etaMapRefined) {
889         if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
890           AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
891           fTPCResponse.SetEtaCorrMap(0x0);
892         }
893         else {
894           AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s", 
895                        refineFactorMapX, refineFactorMapY, fOADBPath.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle()));
896         }
897         
898         delete etaMapRefined;
899       }
900       else {
901         AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun));
902       }
903     }
904   }
905   
906   // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
907   AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass)); 
908   
909   statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
910                                              Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
911   if (statusCont) {
912     AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
913   }
914   else {
915     AliInfo(Form("Loading TPC eta sigma map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
916     
917     TObjArray* etaSigmaPars = 0x0;
918     
919     if (fIsMC && !fTuneMConData) {
920       TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
921       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
922       if (!etaSigmaPars) {
923         // Try default object
924         etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(defaultObj.Data()));
925       }
926     }
927     else {
928       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetObject(fRun, defaultObj.Data()));
929     }
930     
931     if (!etaSigmaPars) {
932       AliError(Form("TPC eta sigma parametrisation not found for run %d -> Using old sigma parametrisation!!!", fRun));
933     }
934     else {
935       TH2D* etaSigmaPar1Map = dynamic_cast<TH2D *>(etaSigmaPars->FindObject("sigmaPar1Map"));
936       TNamed* sigmaPar0Info = dynamic_cast<TNamed *>(etaSigmaPars->FindObject("sigmaPar0"));
937       Double_t sigmaPar0 = 0.0;
938       
939       if (sigmaPar0Info) {
940         TString sigmaPar0String = sigmaPar0Info->GetTitle();
941         sigmaPar0 = sigmaPar0String.Atof();
942       }
943       else {
944         // Something is weired because the object for parameter 0 could not be loaded -> New sigma parametrisation can not be used!
945         etaSigmaPar1Map = 0x0;
946       }
947       
948       TH2D* etaSigmaPar1MapRefined = RefineHistoViaLinearInterpolation(etaSigmaPar1Map, refineFactorSigmaMapX, refineFactorSigmaMapY);
949       
950       
951       if (etaSigmaPar1MapRefined) {
952         if (!fTPCResponse.SetSigmaParams(etaSigmaPar1MapRefined, sigmaPar0)) {
953           AliError(Form("Failed to set TPC eta sigma map for run %d -> Using old sigma parametrisation!!!", fRun));
954           fTPCResponse.SetSigmaParams(0x0, 0);
955         }
956         else {
957           AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s", 
958                        refineFactorSigmaMapX, refineFactorSigmaMapY, fOADBPath.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle()));
959         }
960         
961         delete etaSigmaPar1MapRefined;
962       }
963       else {
964         AliError(Form("Failed to set TPC eta sigma map for run %d (map was loaded, but couldn't be refined) -> Using old sigma parametrisation!!!",
965                       fRun));
966       }
967     }
968   }
969 }
970
971 //______________________________________________________________________________
972 void AliPIDResponse::SetTPCPidResponseMaster()
973 {
974   //
975   // Load the TPC pid response functions from the OADB
976   // Load the TPC voltage maps from OADB
977   //
978   //don't load twice for the moment
979    if (fArrPidResponseMaster) return;
980  
981
982   //reset the PID response functions
983   delete fArrPidResponseMaster;
984   fArrPidResponseMaster=NULL;
985   
986   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
987   TFile *f=NULL;
988   if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
989   
990   TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
991   f=TFile::Open(fileNamePIDresponse.Data());
992   if (f && f->IsOpen() && !f->IsZombie()){
993     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
994   }
995   delete f;
996
997   TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
998   f=TFile::Open(fileNameVoltageMaps.Data());
999   if (f && f->IsOpen() && !f->IsZombie()){
1000     fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
1001   }
1002   delete f;
1003   
1004   if (!fArrPidResponseMaster){
1005     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
1006     return;
1007   }
1008   fArrPidResponseMaster->SetOwner();
1009
1010   if (!fOADBvoltageMaps)
1011   {
1012     AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
1013   }
1014   fArrPidResponseMaster->SetOwner();
1015 }
1016
1017 //______________________________________________________________________________
1018 void AliPIDResponse::SetTPCParametrisation()
1019 {
1020   //
1021   // Change BB parametrisation for current run
1022   //
1023   
1024   //
1025   //reset old splines
1026   //
1027   fTPCResponse.ResetSplines();
1028   
1029   if (fLHCperiod.IsNull()) {
1030     AliError("No period set, not changing parametrisation");
1031     return;
1032   }
1033   
1034   //
1035   // Set default parametrisations for data and MC
1036   //
1037   
1038   //data type
1039   TString datatype="DATA";
1040   //in case of mc fRecoPass is per default 1
1041   if (fIsMC) {
1042       if(!fTuneMConData) datatype="MC";
1043       fRecoPass=1;
1044   }
1045
1046   // period
1047   TString period=fLHCperiod;
1048   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
1049
1050   Int_t recopass = fRecoPass;
1051   if(fTuneMConData) recopass = fRecoPassUser;
1052     
1053   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1054   Bool_t found=kFALSE;
1055   //
1056   //set the new PID splines
1057   //
1058   if (fArrPidResponseMaster){
1059     //for MC don't use period information
1060     //if (fIsMC) period="[A-Z0-9]*";
1061     //for MC use MC period information
1062     //pattern for the default entry (valid for all particles)
1063     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1064
1065     //find particle id and gain scenario
1066     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
1067     {
1068       TObject *grAll=NULL;
1069       TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
1070       gainScenario.ToUpper();
1071       //loop over entries and filter them
1072       for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
1073       {
1074         TObject *responseFunction=fArrPidResponseMaster->At(iresp);
1075         if (responseFunction==NULL) continue;
1076         TString responseName=responseFunction->GetName();
1077          
1078         if (!reg.MatchB(responseName)) continue;
1079
1080         TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
1081         TObject* tmp=NULL;
1082         tmp=arr->At(1); if (!tmp) continue;
1083         TString particleName=tmp->GetName();
1084         tmp=arr->At(3); if (!tmp) continue;
1085         TString gainScenarioName=tmp->GetName();
1086         delete arr;
1087         if (particleName.IsNull()) continue;
1088         if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
1089         else 
1090         {
1091           for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1092           {
1093             TString particle=AliPID::ParticleName(ispec);
1094             particle.ToUpper();
1095             //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
1096             if ( particle == particleName && gainScenario == gainScenarioName )
1097             {
1098               fTPCResponse.SetResponseFunction( responseFunction,
1099                                                 (AliPID::EParticleType)ispec,
1100                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1101               fTPCResponse.SetUseDatabase(kTRUE);
1102               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
1103               found=kTRUE;
1104               break;
1105             }
1106           }
1107         }
1108       }
1109       
1110       // Retrieve responsefunction for pions - will (if available) be used for muons if there are no dedicated muon splines.
1111       // For light nuclei, try to set the proton spline, if no dedicated splines are available.
1112       // In both cases: Use default splines, if no dedicated splines and no pion/proton splines are available.
1113       TObject* responseFunctionPion = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kPion,                             
1114                                                                         (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1115       TObject* responseFunctionProton = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kProton,                             
1116                                                                           (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1117       
1118       for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1119       {
1120         if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
1121           (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
1122         {
1123           if (ispec == AliPID::kMuon) { // Muons
1124             if (responseFunctionPion) {
1125               fTPCResponse.SetResponseFunction( responseFunctionPion,
1126                                                 (AliPID::EParticleType)ispec,
1127                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1128               fTPCResponse.SetUseDatabase(kTRUE);
1129               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionPion->GetName()));
1130               found=kTRUE;  
1131             }
1132             else if (grAll) {
1133               fTPCResponse.SetResponseFunction( grAll,
1134                                                 (AliPID::EParticleType)ispec,
1135                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1136               fTPCResponse.SetUseDatabase(kTRUE);
1137               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1138               found=kTRUE;
1139             }
1140             //else
1141             //  AliError(Form("No splines found for muons (also no pion splines and no default splines) for gain scenario %d!", igainScenario));
1142           }
1143           else if (ispec >= AliPID::kSPECIES) { // Light nuclei
1144             if (responseFunctionProton) {
1145               fTPCResponse.SetResponseFunction( responseFunctionProton,
1146                                                 (AliPID::EParticleType)ispec,
1147                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1148               fTPCResponse.SetUseDatabase(kTRUE);
1149               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunctionProton->GetName()));
1150               found=kTRUE;  
1151             }
1152             else if (grAll) {
1153               fTPCResponse.SetResponseFunction( grAll,
1154                                                 (AliPID::EParticleType)ispec,
1155                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1156               fTPCResponse.SetUseDatabase(kTRUE);
1157               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1158               found=kTRUE;
1159             }
1160             //else
1161             //  AliError(Form("No splines found for species %d (also no proton splines and no default splines) for gain scenario %d!",
1162             //                ispec, igainScenario));
1163           }
1164         }
1165       }
1166     }
1167   }
1168   else AliInfo("no fArrPidResponseMaster");
1169
1170   if (!found){
1171     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1172   }
1173
1174   //
1175   // Setup resolution parametrisation
1176   //
1177   
1178   //default
1179   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1180   
1181   if (fRun>=122195){
1182     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1183   }
1184
1185   if (fRun>=186636){
1186 //   if (fRun>=188356){
1187     fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1188   }
1189   
1190   if (fArrPidResponseMaster)
1191   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1192   
1193   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
1194
1195   //read in the voltage map
1196   TVectorF* gsm = 0x0;
1197   if (fOADBvoltageMaps) gsm=dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
1198   if (gsm) 
1199   {
1200     fTPCResponse.SetVoltageMap(*gsm);
1201     TString vals;
1202     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1203     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1204     AliInfo(vals.Data());
1205     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1206     AliInfo(vals.Data());
1207     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1208     AliInfo(vals.Data());
1209     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1210     AliInfo(vals.Data());
1211   }
1212   else AliInfo("no voltage map, ideal default assumed");
1213 }
1214
1215 //______________________________________________________________________________
1216 void AliPIDResponse::SetTRDPidResponseMaster()
1217 {
1218   //
1219   // Load the TRD pid params and references from the OADB
1220   //
1221   if(fTRDPIDResponseObject) return;
1222   AliOADBContainer contParams("contParams"); 
1223
1224   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1225   if(statusResponse){
1226     AliError("Failed initializing PID Response Object from OADB");
1227   } else {
1228     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1229     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1230     if(!fTRDPIDResponseObject){
1231       AliError(Form("TRD Response not found in run %d", fRun));
1232     }
1233   }
1234 }
1235
1236 //______________________________________________________________________________
1237 void AliPIDResponse::InitializeTRDResponse(){
1238   //
1239   // Set PID Params and references to the TRD PID response
1240   // 
1241     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1242 }
1243
1244 //______________________________________________________________________________
1245 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1246
1247     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1248         // backward compatibility for setting with 8 slices
1249         TRDslicesForPID[0] = 0;
1250         TRDslicesForPID[1] = 7;
1251     }
1252     else{
1253         if(method==AliTRDPIDResponse::kLQ1D){
1254             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1255             TRDslicesForPID[1] = 0;
1256         }
1257         if(method==AliTRDPIDResponse::kLQ2D){
1258             TRDslicesForPID[0] = 1;
1259             TRDslicesForPID[1] = 7;
1260         }
1261     }
1262     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
1263 }
1264
1265 //______________________________________________________________________________
1266 void AliPIDResponse::SetTOFPidResponseMaster()
1267 {
1268   //
1269   // Load the TOF pid params from the OADB
1270   //
1271
1272   if (fTOFPIDParams) delete fTOFPIDParams;
1273   fTOFPIDParams=NULL;
1274
1275   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1276   if (oadbf && oadbf->IsOpen()) {
1277     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1278     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1279     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
1280     oadbf->Close();
1281     delete oadbc;
1282   }
1283   delete oadbf;
1284
1285   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1286 }
1287
1288 //______________________________________________________________________________
1289 void AliPIDResponse::InitializeTOFResponse(){
1290   //
1291   // Set PID Params to the TOF PID response
1292   //
1293
1294   AliInfo("TOF PID Params loaded from OADB");
1295   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1296   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1297   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1298                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1299   
1300   for (Int_t i=0;i<4;i++) {
1301     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1302   }
1303   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1304
1305   AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1306   Float_t t0Spread[4];
1307   for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1308   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]));
1309   Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1310   Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1311   if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1312     fResT0AC=t0Spread[3];
1313     fResT0A=TMath::Sqrt(a);
1314     fResT0C=TMath::Sqrt(c);
1315   } else {
1316     AliInfo("  TZERO spreads not present or inconsistent, loading default");
1317     fResT0A=75.;
1318     fResT0C=65.;
1319     fResT0AC=55.;
1320   }
1321   AliInfo(Form("  TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1322
1323 }
1324
1325
1326 //______________________________________________________________________________
1327 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1328   //
1329   // Check whether track is identified as electron under a given electron efficiency hypothesis
1330     //
1331
1332   Double_t probs[AliPID::kSPECIES];
1333   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1334
1335   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1336   // Take mean of the TRD momenta in the given tracklets
1337   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1338   Int_t nmomenta = 0;
1339   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1340     if(vtrack->GetTRDmomentum(iPl) > 0.){
1341       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
1342     }
1343   }
1344   p = TMath::Mean(nmomenta, trdmomenta);
1345
1346   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1347 }
1348
1349 //______________________________________________________________________________
1350 void AliPIDResponse::SetEMCALPidResponseMaster()
1351 {
1352   //
1353   // Load the EMCAL pid response functions from the OADB
1354   //
1355   TObjArray* fEMCALPIDParamsRun      = NULL;
1356   TObjArray* fEMCALPIDParamsPass     = NULL;
1357
1358   if(fEMCALPIDParams) return;
1359   AliOADBContainer contParams("contParams"); 
1360
1361   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1362   if(statusPars){
1363     AliError("Failed initializing PID Params from OADB");
1364   } 
1365   else {
1366     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1367
1368     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1369     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1370     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1371
1372     if(!fEMCALPIDParams){
1373       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1374       AliInfo("Will take the standard LHC11d instead ...");
1375
1376       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1377       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1378       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1379
1380       if(!fEMCALPIDParams){
1381         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));     
1382       }
1383     }
1384   }
1385 }
1386
1387 //______________________________________________________________________________
1388 void AliPIDResponse::InitializeEMCALResponse(){
1389   //
1390   // Set PID Params to the EMCAL PID response
1391   // 
1392   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1393
1394 }
1395
1396 //______________________________________________________________________________
1397 void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
1398 {
1399   //
1400   // create detector PID information and setup the transient pointer in the track
1401   //
1402   
1403   // check if detector number is inside accepted range
1404   if (detector == kNdetectors) return;
1405   
1406   // get detector pid
1407   AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1408   if (!detPID) {
1409     detPID=new AliDetectorPID;
1410     (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1411   }
1412   
1413   //check if values exist
1414   if (detPID->HasRawProbability(detector) && detPID->HasNumberOfSigmas(detector)) return;
1415   
1416   //TODO: which particles to include? See also the loops below...
1417   Double_t values[AliPID::kSPECIESC]={0};
1418
1419   //probabilities
1420   EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1421   detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1422   
1423   //nsigmas
1424   for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1425     values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
1426   // the pid status is the same for probabilities and nSigmas, so it is
1427   // fine to use the one from the probabilities also here
1428   detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC, status);
1429   
1430 }
1431
1432 //______________________________________________________________________________
1433 void AliPIDResponse::FillTrackDetectorPID()
1434 {
1435   //
1436   // create detector PID information and setup the transient pointer in the track
1437   //
1438
1439   if (!fCurrentEvent) return;
1440   
1441   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1442     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1443     if (!track) continue;
1444
1445     for (Int_t idet=0; idet<kNdetectors; ++idet){
1446       FillTrackDetectorPID(track, (EDetector)idet);
1447     }
1448   }
1449 }
1450
1451 //______________________________________________________________________________
1452 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1453   //
1454   // Set TOF response function
1455   // Input option for event_time used
1456   //
1457   
1458     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1459     if(t0spread < 10) t0spread = 80;
1460
1461     // T0 from TOF algorithm
1462
1463     Bool_t flagT0TOF=kFALSE;
1464     Bool_t flagT0T0=kFALSE;
1465     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1466     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1467     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1468
1469     // T0-TOF arrays
1470     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1471     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1472     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1473       estimatedT0event[i]=0.0;
1474       estimatedT0resolution[i]=0.0;
1475       startTimeMask[i] = 0;
1476     }
1477
1478     Float_t resT0A=fResT0A;
1479     Float_t resT0C=fResT0C;
1480     Float_t resT0AC=fResT0AC;
1481     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1482         flagT0T0=kTRUE;
1483     }
1484
1485
1486     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1487
1488     if (tofHeader) { // read global info and T0-TOF
1489       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1490       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1491       if(t0spread < 10) t0spread = 80;
1492
1493       flagT0TOF=kTRUE;
1494       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1495         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1496         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1497         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1498       }
1499
1500       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1501       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1502       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1503       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1504         Int_t icurrent = (Int_t)ibin->GetAt(j);
1505         startTime[icurrent]=t0Bin->GetAt(j);
1506         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1507         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1508       }
1509     }
1510
1511     // for cut of 3 sigma on t0 spread
1512     Float_t t0cut = 3 * t0spread;
1513     if(t0cut < 500) t0cut = 500;
1514
1515     if(option == kFILL_T0){ // T0-FILL is used
1516         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1517           estimatedT0event[i]=0.0;
1518           estimatedT0resolution[i]=t0spread;
1519         }
1520         fTOFResponse.SetT0event(estimatedT0event);
1521         fTOFResponse.SetT0resolution(estimatedT0resolution);
1522     }
1523
1524     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1525         if(flagT0TOF){
1526             fTOFResponse.SetT0event(startTime);
1527             fTOFResponse.SetT0resolution(startTimeRes);
1528             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1529               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1530               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1531             }
1532         }
1533         else{
1534             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1535               estimatedT0event[i]=0.0;
1536               estimatedT0resolution[i]=t0spread;
1537               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1538             }
1539             fTOFResponse.SetT0event(estimatedT0event);
1540             fTOFResponse.SetT0resolution(estimatedT0resolution);
1541         }
1542     }
1543     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1544         Float_t t0AC=-10000;
1545         Float_t t0A=-10000;
1546         Float_t t0C=-10000;
1547         if(flagT0T0){
1548             t0A= vevent->GetT0TOF()[1];
1549             t0C= vevent->GetT0TOF()[2];
1550         //      t0AC= vevent->GetT0TOF()[0];
1551         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1552         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1553         t0AC /= resT0AC*resT0AC;
1554         }
1555
1556         Float_t t0t0Best = 0;
1557         Float_t t0t0BestRes = 9999;
1558         Int_t t0used=0;
1559         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1560             t0t0Best = t0AC;
1561             t0t0BestRes = resT0AC;
1562             t0used=6;
1563         }
1564         else if(TMath::Abs(t0C) < t0cut){
1565             t0t0Best = t0C;
1566             t0t0BestRes = resT0C;
1567             t0used=4;
1568         }
1569         else if(TMath::Abs(t0A) < t0cut){
1570             t0t0Best = t0A;
1571             t0t0BestRes = resT0A;
1572             t0used=2;
1573         }
1574
1575         if(flagT0TOF){ // if T0-TOF info is available
1576             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1577                 if(t0t0BestRes < 999){
1578                   if(startTimeRes[i] < t0spread){
1579                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1580                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1581                     estimatedT0event[i]=t0best / wtot;
1582                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1583                     startTimeMask[i] = t0used+1;
1584                   }
1585                   else {
1586                     estimatedT0event[i]=t0t0Best;
1587                     estimatedT0resolution[i]=t0t0BestRes;
1588                     startTimeMask[i] = t0used;
1589                   }
1590                 }
1591                 else{
1592                   estimatedT0event[i]=startTime[i];
1593                   estimatedT0resolution[i]=startTimeRes[i];
1594                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1595                 }
1596                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1597             }
1598             fTOFResponse.SetT0event(estimatedT0event);
1599             fTOFResponse.SetT0resolution(estimatedT0resolution);
1600         }
1601         else{ // if no T0-TOF info is available
1602             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1603               fTOFResponse.SetT0binMask(i,t0used);
1604               if(t0t0BestRes < 999){
1605                 estimatedT0event[i]=t0t0Best;
1606                 estimatedT0resolution[i]=t0t0BestRes;
1607               }
1608               else{
1609                 estimatedT0event[i]=0.0;
1610                 estimatedT0resolution[i]=t0spread;
1611               }
1612             }
1613             fTOFResponse.SetT0event(estimatedT0event);
1614             fTOFResponse.SetT0resolution(estimatedT0resolution);
1615         }
1616     }
1617
1618     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1619         Float_t t0AC=-10000;
1620         Float_t t0A=-10000;
1621         Float_t t0C=-10000;
1622         if(flagT0T0){
1623             t0A= vevent->GetT0TOF()[1];
1624             t0C= vevent->GetT0TOF()[2];
1625         //      t0AC= vevent->GetT0TOF()[0];
1626         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1627         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1628         t0AC /= resT0AC*resT0AC;
1629         }
1630
1631         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1632             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1633               estimatedT0event[i]=t0AC;
1634               estimatedT0resolution[i]=resT0AC;
1635               fTOFResponse.SetT0binMask(i,6);
1636             }
1637         }
1638         else if(TMath::Abs(t0C) < t0cut){
1639             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1640               estimatedT0event[i]=t0C;
1641               estimatedT0resolution[i]=resT0C;
1642               fTOFResponse.SetT0binMask(i,4);
1643             }
1644         }
1645         else if(TMath::Abs(t0A) < t0cut){
1646             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1647               estimatedT0event[i]=t0A;
1648               estimatedT0resolution[i]=resT0A;
1649               fTOFResponse.SetT0binMask(i,2);
1650             }
1651         }
1652         else{
1653             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1654               estimatedT0event[i]=0.0;
1655               estimatedT0resolution[i]=t0spread;
1656               fTOFResponse.SetT0binMask(i,0);
1657             }
1658         }
1659         fTOFResponse.SetT0event(estimatedT0event);
1660         fTOFResponse.SetT0resolution(estimatedT0resolution);
1661     }
1662     delete [] startTime;
1663     delete [] startTimeRes;
1664     delete [] startTimeMask;
1665     delete [] estimatedT0event;
1666     delete [] estimatedT0resolution;
1667 }
1668
1669 //______________________________________________________________________________
1670 // private non cached versions of the PID calculation
1671 //
1672
1673
1674 //______________________________________________________________________________
1675 Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detector, const AliVParticle *vtrack, AliPID::EParticleType type) const
1676 {
1677   //
1678   // NumberOfSigmas for 'detCode'
1679   //
1680
1681   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
1682   
1683   switch (detector){
1684     case kITS:   return GetNumberOfSigmasITS(track, type); break;
1685     case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
1686     case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
1687     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
1688     default: return -999.;
1689   }
1690
1691   return -999.;
1692 }
1693
1694 //______________________________________________________________________________
1695 Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
1696 {
1697   //
1698   // Calculate the number of sigmas in the ITS
1699   //
1700   
1701   AliVTrack *track=(AliVTrack*)vtrack;
1702
1703   const EDetPidStatus pidStatus=GetITSPIDStatus(track);
1704   if (pidStatus!=kDetPidOk) return -999.;
1705     
1706   UChar_t clumap=track->GetITSClusterMap();
1707   Int_t nPointsForPid=0;
1708   for(Int_t i=2; i<6; i++){
1709     if(clumap&(1<<i)) ++nPointsForPid;
1710   }
1711   Float_t mom=track->P();
1712   
1713   //check for ITS standalone tracks
1714   Bool_t isSA=kTRUE;
1715   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1716   
1717   const Float_t dEdx=track->GetITSsignal();
1718
1719   //TODO: in case of the electron, use the SA parametrisation,
1720   //      this needs to be changed if ITS provides a parametrisation
1721   //      for electrons also for ITS+TPC tracks
1722   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
1723 }
1724
1725 //______________________________________________________________________________
1726 Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
1727 {
1728   //
1729   // Calculate the number of sigmas in the TPC
1730   //
1731   
1732   AliVTrack *track=(AliVTrack*)vtrack;
1733
1734   const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
1735   if (pidStatus!=kDetPidOk) return -999.;
1736   
1737   Double_t nSigma = -999.;
1738   
1739   if (fTuneMConData)
1740     this->GetTPCsignalTunedOnData(track);
1741   
1742   nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1743   
1744   return nSigma;
1745 }
1746
1747 //______________________________________________________________________________
1748 Float_t AliPIDResponse::GetNumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
1749 {
1750   //
1751   // Calculate the number of sigmas in the TOF
1752   //
1753   
1754   AliVTrack *track=(AliVTrack*)vtrack;
1755
1756   const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
1757   if (pidStatus!=kDetPidOk) return -999.;
1758
1759   
1760   return GetNumberOfSigmasTOFold(vtrack, type);
1761 }
1762
1763 //______________________________________________________________________________
1764 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1765 {
1766   //
1767   // Calculate the number of sigmas in the EMCAL
1768   //
1769   
1770   AliVTrack *track=(AliVTrack*)vtrack;
1771
1772   const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
1773   if (pidStatus!=kDetPidOk) return -999.;
1774
1775   const Int_t nMatchClus = track->GetEMCALcluster();
1776   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1777   
1778   const Double_t mom    = track->P();
1779   const Double_t pt     = track->Pt();
1780   const Int_t    charge = track->Charge();
1781   const Double_t fClsE  = matchedClus->E();
1782   const Double_t EovP   = fClsE/mom;
1783   
1784   return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1785 }
1786
1787
1788 //______________________________________________________________________________
1789 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1790 {
1791   //
1792   // Compute PID response of 'detCode'
1793   //
1794
1795   switch (detCode){
1796     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1797     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1798     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1799     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1800     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1801     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1802     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1803     default: return kDetNoSignal;
1804   }
1805 }
1806
1807 //______________________________________________________________________________
1808 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1809 {
1810   //
1811   // Compute PID response for the ITS
1812   //
1813   
1814   // set flat distribution (no decision)
1815   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1816   
1817   const EDetPidStatus pidStatus=GetITSPIDStatus(track);
1818   if (pidStatus!=kDetPidOk) return pidStatus;
1819   
1820   if (track->GetDetectorPID()){
1821     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
1822   }
1823   
1824   //check for ITS standalone tracks
1825   Bool_t isSA=kTRUE;
1826   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1827
1828   Double_t mom=track->P();
1829   Double_t dedx=track->GetITSsignal();
1830   Double_t momITS=mom;
1831   UChar_t clumap=track->GetITSClusterMap();
1832   Int_t nPointsForPid=0;
1833   for(Int_t i=2; i<6; i++){
1834     if(clumap&(1<<i)) ++nPointsForPid;
1835   }
1836
1837   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1838   for (Int_t j=0; j<nSpecies; j++) {
1839     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1840     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1841     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1842     //TODO: in case of the electron, use the SA parametrisation,
1843     //      this needs to be changed if ITS provides a parametrisation
1844     //      for electrons also for ITS+TPC tracks
1845     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1846     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1847       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1848     } else {
1849       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1850       mismatch=kFALSE;
1851     }
1852   }
1853
1854   if (mismatch){
1855     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1856   }
1857
1858   return kDetPidOk;
1859 }
1860 //______________________________________________________________________________
1861 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1862 {
1863   //
1864   // Compute PID response for the TPC
1865   //
1866   
1867   // set flat distribution (no decision)
1868   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1869   
1870   const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
1871   if (pidStatus!=kDetPidOk) return pidStatus;
1872   
1873   Double_t dedx=track->GetTPCsignal();
1874   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1875   
1876   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1877   
1878   Double_t bethe = 0.;
1879   Double_t sigma = 0.;
1880   
1881   for (Int_t j=0; j<nSpecies; j++) {
1882     AliPID::EParticleType type=AliPID::EParticleType(j);
1883     
1884     bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1885     sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1886     
1887     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1888       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1889     } else {
1890       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1891       mismatch=kFALSE;
1892     }
1893   }
1894   
1895   if (mismatch){
1896     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1897   }
1898   
1899   return kDetPidOk;
1900 }
1901 //______________________________________________________________________________
1902 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1903 {
1904   //
1905   // Compute PID probabilities for TOF
1906   //
1907   
1908   // set flat distribution (no decision)
1909   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1910   
1911   const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
1912   if (pidStatus!=kDetPidOk) return pidStatus;
1913
1914   const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
1915   
1916   for (Int_t j=0; j<nSpecies; j++) {
1917     AliPID::EParticleType type=AliPID::EParticleType(j);
1918     const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
1919     
1920     const Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
1921     const Double_t sig     = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
1922     if (TMath::Abs(nsigmas) > (fRange+2)) {
1923       if(nsigmas < fTOFtail)
1924         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
1925       else
1926         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
1927     } else{
1928       if(nsigmas < fTOFtail)
1929         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
1930       else
1931         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
1932     }    
1933   }
1934   
1935   return kDetPidOk;
1936 }
1937 //______________________________________________________________________________
1938 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
1939 {
1940   //
1941   // Compute PID probabilities for the TRD
1942   //
1943   
1944   // set flat distribution (no decision)
1945   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1946   
1947   const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
1948   if (pidStatus!=kDetPidOk) return pidStatus;
1949
1950   UInt_t TRDslicesForPID[2];
1951   SetTRDSlices(TRDslicesForPID,PIDmethod);
1952   
1953   Float_t mom[6]={0.};
1954   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
1955   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
1956   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
1957   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
1958     mom[ilayer] = track->GetTRDmomentum(ilayer);
1959     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
1960       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
1961     }
1962   }
1963   
1964   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
1965   return kDetPidOk;
1966 }
1967
1968 //______________________________________________________________________________
1969 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1970 {
1971   //
1972   // Compute PID response for the EMCAL
1973   //
1974   
1975   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1976
1977   const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
1978   if (pidStatus!=kDetPidOk) return pidStatus;
1979
1980   const Int_t nMatchClus = track->GetEMCALcluster();
1981   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1982   
1983   const Double_t mom    = track->P();
1984   const Double_t pt     = track->Pt();
1985   const Int_t    charge = track->Charge();
1986   const Double_t fClsE  = matchedClus->E();
1987   const Double_t EovP   = fClsE/mom;
1988   
1989   // compute the probabilities
1990   fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p);
1991   return kDetPidOk;
1992 }
1993
1994 //______________________________________________________________________________
1995 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
1996 {
1997   //
1998   // Compute PID response for the PHOS
1999   //
2000   
2001   // set flat distribution (no decision)
2002   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2003   return kDetNoSignal;
2004 }
2005
2006 //______________________________________________________________________________
2007 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2008 {
2009   //
2010   // Compute PID response for the HMPID
2011   //
2012   
2013   // set flat distribution (no decision)
2014   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2015   
2016   const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
2017   if (pidStatus!=kDetPidOk) return pidStatus;
2018   
2019   track->GetHMPIDpid(p);
2020   
2021   return kDetPidOk;
2022 }
2023
2024 //______________________________________________________________________________
2025 AliPIDResponse::EDetPidStatus AliPIDResponse::GetITSPIDStatus(const AliVTrack *track) const
2026 {
2027   // compute ITS pid status
2028
2029   // check status bits
2030   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
2031     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
2032
2033   const Float_t dEdx=track->GetITSsignal();
2034   if (dEdx<=0) return kDetNoSignal;
2035   
2036   // requite at least 3 pid clusters
2037   const UChar_t clumap=track->GetITSClusterMap();
2038   Int_t nPointsForPid=0;
2039   for(Int_t i=2; i<6; i++){
2040     if(clumap&(1<<i)) ++nPointsForPid;
2041   }
2042   
2043   if(nPointsForPid<3) { 
2044     return kDetNoSignal;
2045   }
2046   
2047   return kDetPidOk;
2048 }
2049
2050 //______________________________________________________________________________
2051 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetTPCPIDStatus(const AliVTrack *track) const
2052 {
2053   // compute TPC pid status
2054   
2055   // check quality of the track
2056   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
2057
2058   // check pid values
2059   const Double_t dedx=track->GetTPCsignal();
2060   const UShort_t signalN=track->GetTPCsignalN();
2061   if (signalN<10 || dedx<10) return kDetNoSignal;
2062
2063   if (!(fArrPidResponseMaster && fArrPidResponseMaster->At(AliPID::kPion))) return kDetNoParams;
2064   
2065   return kDetPidOk;
2066 }
2067
2068 //______________________________________________________________________________
2069 AliPIDResponse::EDetPidStatus AliPIDResponse::GetTRDPIDStatus(const AliVTrack *track) const
2070 {
2071   // compute TRD pid status
2072
2073   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
2074   return kDetPidOk;
2075 }
2076
2077 //______________________________________________________________________________
2078 AliPIDResponse::EDetPidStatus AliPIDResponse::GetTOFPIDStatus(const AliVTrack *track) const
2079 {
2080   // compute TOF pid status
2081
2082   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
2083   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
2084
2085   return kDetPidOk;
2086 }
2087
2088 //______________________________________________________________________________
2089 Float_t AliPIDResponse::GetTOFMismatchProbability(const AliVTrack *track) const
2090 {
2091   // compute mismatch probability cross-checking at 5 sigmas with TPC
2092   // currently just implemented as a 5 sigma compatibility cut
2093
2094   // check pid status
2095   const EDetPidStatus tofStatus=GetTOFPIDStatus(track);
2096   if (tofStatus!=kDetPidOk) return 0.;
2097
2098   //mismatch
2099   const EDetPidStatus tpcStatus=GetTPCPIDStatus(track);
2100   if (tpcStatus!=kDetPidOk) return 0.;
2101   
2102   const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2103   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
2104   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
2105     AliPID::EParticleType type=AliPID::EParticleType(j);
2106     const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
2107     
2108     if (TMath::Abs(nsigmas)<5.){
2109       const Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
2110       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
2111     }
2112   }
2113   
2114   if (mismatch){
2115     return 1.;
2116   }
2117   
2118   return 0.;
2119 }
2120
2121
2122
2123 //______________________________________________________________________________
2124 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetHMPIDPIDStatus(const AliVTrack *track) const
2125 {
2126   // compute HMPID pid status
2127   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
2128   return kDetPidOk;
2129 }
2130
2131 //______________________________________________________________________________
2132 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetPHOSPIDStatus(const AliVTrack */*track*/) const
2133 {
2134   // compute PHOS pid status
2135   return kDetNoSignal;  
2136 }
2137
2138 //______________________________________________________________________________
2139 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetEMCALPIDStatus(const AliVTrack *track) const
2140 {
2141   // compute EMCAL pid status
2142
2143
2144   // Track matching
2145   const Int_t nMatchClus = track->GetEMCALcluster();
2146   if (nMatchClus<0) return kDetNoSignal;
2147
2148   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2149
2150   if (!(matchedClus && matchedClus->IsEMCAL())) return kDetNoSignal;
2151
2152   const Int_t charge = track->Charge();
2153   if (TMath::Abs(charge)!=1) return kDetNoSignal;
2154
2155   if (!(fEMCALPIDParams && fEMCALPIDParams->At(AliPID::kElectron))) return kDetNoParams;
2156   
2157   return kDetPidOk;
2158
2159 }
2160
2161 //______________________________________________________________________________
2162 AliPIDResponse::EDetPidStatus AliPIDResponse::GetPIDStatus(EDetector detector, const AliVTrack *track) const
2163 {
2164   //
2165   // check pid status for a track
2166   //
2167
2168   switch (detector){
2169     case kITS:   return GetITSPIDStatus(track);   break;
2170     case kTPC:   return GetTPCPIDStatus(track);   break;
2171     case kTRD:   return GetTRDPIDStatus(track);   break;
2172     case kTOF:   return GetTOFPIDStatus(track);   break;
2173     case kPHOS:  return GetPHOSPIDStatus(track);  break;
2174     case kEMCAL: return GetEMCALPIDStatus(track); break;
2175     case kHMPID: return GetHMPIDPIDStatus(track); break;
2176     default: return kDetNoSignal;
2177   }
2178   return kDetNoSignal;
2179   
2180 }