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