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