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