- Code cleanup
[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) const
378 {
379   //
380   //
381   //
382   val=-9999.;
383   switch (detector){
384     case kITS:   return GetSignalDeltaITS(track,type,val); break;
385     case kTPC:   return GetSignalDeltaTPC(track,type,val); break;
386     case kTOF:   return GetSignalDeltaTOF(track,type,val); break;
387     case kHMPID: return GetSignalDeltaHMPID(track,type,val); 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) const
395 {
396   //
397   //
398   //
399   Double_t val=-9999.;
400   EDetPidStatus stat=GetSignalDelta(detCode, track, type, val);
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   return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1803 }
1804
1805 //______________________________________________________________________________
1806 Float_t AliPIDResponse::GetNumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
1807 {
1808   //
1809   // Calculate the number of sigmas in the TOF
1810   //
1811   
1812   AliVTrack *track=(AliVTrack*)vtrack;
1813
1814   const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
1815   if (pidStatus!=kDetPidOk) return -999.;
1816   
1817   return GetNumberOfSigmasTOFold(vtrack, type);
1818 }
1819 //______________________________________________________________________________
1820
1821 Float_t AliPIDResponse::GetNumberOfSigmasHMPID(const AliVParticle *vtrack, AliPID::EParticleType type) const
1822 {
1823   //
1824   // Calculate the number of sigmas in the HMPID
1825   //  
1826   AliVTrack *track=(AliVTrack*)vtrack;
1827     
1828   const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
1829   if (pidStatus!=kDetPidOk) return -999.; 
1830   
1831   return fHMPIDResponse.GetNumberOfSigmas(track, type);
1832 }
1833
1834 //______________________________________________________________________________
1835 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1836 {
1837   //
1838   // Calculate the number of sigmas in the EMCAL
1839   //
1840   
1841   AliVTrack *track=(AliVTrack*)vtrack;
1842
1843   const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
1844   if (pidStatus!=kDetPidOk) return -999.;
1845
1846   const Int_t nMatchClus = track->GetEMCALcluster();
1847   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1848   
1849   const Double_t mom    = track->P();
1850   const Double_t pt     = track->Pt();
1851   const Int_t    charge = track->Charge();
1852   const Double_t fClsE  = matchedClus->E();
1853   const Double_t EovP   = fClsE/mom;
1854   
1855   return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1856 }
1857
1858 //______________________________________________________________________________
1859 AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaITS(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
1860 {
1861   //
1862   // Signal minus expected Signal for ITS
1863   //
1864   AliVTrack *track=(AliVTrack*)vtrack;
1865   val=fITSResponse.GetSignalDelta(track,type);
1866   
1867   return GetITSPIDStatus(track);
1868 }
1869
1870 //______________________________________________________________________________
1871 AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
1872 {
1873   //
1874   // Signal minus expected Signal for TPC
1875   //
1876   AliVTrack *track=(AliVTrack*)vtrack;
1877   val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1878   
1879   return GetTPCPIDStatus(track);
1880 }
1881
1882 //______________________________________________________________________________
1883 AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTOF(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
1884 {
1885   //
1886   // Signal minus expected Signal for TOF
1887   //
1888   AliVTrack *track=(AliVTrack*)vtrack;
1889   val=GetSignalDeltaTOFold(track, type);
1890   
1891   return GetTOFPIDStatus(track);
1892 }
1893
1894 //______________________________________________________________________________
1895 AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaHMPID(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val) const
1896 {
1897   //
1898   // Signal minus expected Signal for HMPID
1899   //
1900   AliVTrack *track=(AliVTrack*)vtrack;
1901   val=fHMPIDResponse.GetSignalDelta(track, type);
1902   
1903   return GetHMPIDPIDStatus(track);
1904 }
1905
1906 //______________________________________________________________________________
1907 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1908 {
1909   //
1910   // Compute PID response of 'detCode'
1911   //
1912
1913   switch (detCode){
1914     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1915     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1916     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1917     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1918     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1919     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1920     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1921     default: return kDetNoSignal;
1922   }
1923 }
1924
1925 //______________________________________________________________________________
1926 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1927 {
1928   //
1929   // Compute PID response for the ITS
1930   //
1931   
1932   // set flat distribution (no decision)
1933   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1934   
1935   const EDetPidStatus pidStatus=GetITSPIDStatus(track);
1936   if (pidStatus!=kDetPidOk) return pidStatus;
1937   
1938   if (track->GetDetectorPID()){
1939     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
1940   }
1941   
1942   //check for ITS standalone tracks
1943   Bool_t isSA=kTRUE;
1944   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1945
1946   Double_t mom=track->P();
1947   Double_t dedx=track->GetITSsignal();
1948   Double_t momITS=mom;
1949   UChar_t clumap=track->GetITSClusterMap();
1950   Int_t nPointsForPid=0;
1951   for(Int_t i=2; i<6; i++){
1952     if(clumap&(1<<i)) ++nPointsForPid;
1953   }
1954
1955   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1956   for (Int_t j=0; j<nSpecies; j++) {
1957     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1958     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1959     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1960     //TODO: in case of the electron, use the SA parametrisation,
1961     //      this needs to be changed if ITS provides a parametrisation
1962     //      for electrons also for ITS+TPC tracks
1963     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1964     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1965       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1966     } else {
1967       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1968       mismatch=kFALSE;
1969     }
1970   }
1971
1972   if (mismatch){
1973     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1974   }
1975
1976   return kDetPidOk;
1977 }
1978 //______________________________________________________________________________
1979 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1980 {
1981   //
1982   // Compute PID response for the TPC
1983   //
1984   
1985   // set flat distribution (no decision)
1986   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1987   
1988   const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
1989   if (pidStatus!=kDetPidOk) return pidStatus;
1990   
1991   Double_t dedx=track->GetTPCsignal();
1992   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1993   
1994   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1995   
1996   Double_t bethe = 0.;
1997   Double_t sigma = 0.;
1998   
1999   for (Int_t j=0; j<nSpecies; j++) {
2000     AliPID::EParticleType type=AliPID::EParticleType(j);
2001     
2002     bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
2003     sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
2004     
2005     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
2006       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
2007     } else {
2008       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
2009       mismatch=kFALSE;
2010     }
2011   }
2012   
2013   if (mismatch){
2014     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2015   }
2016   
2017   return kDetPidOk;
2018 }
2019 //______________________________________________________________________________
2020 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2021 {
2022   //
2023   // Compute PID probabilities for TOF
2024   //
2025   
2026   // set flat distribution (no decision)
2027   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2028   
2029   const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
2030   if (pidStatus!=kDetPidOk) return pidStatus;
2031
2032   const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2033   
2034   for (Int_t j=0; j<nSpecies; j++) {
2035     AliPID::EParticleType type=AliPID::EParticleType(j);
2036     const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
2037     
2038     const Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
2039     const Double_t sig     = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
2040     if (TMath::Abs(nsigmas) > (fRange+2)) {
2041       if(nsigmas < fTOFtail)
2042         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
2043       else
2044         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
2045     } else{
2046       if(nsigmas < fTOFtail)
2047         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
2048       else
2049         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
2050     }    
2051   }
2052   
2053   return kDetPidOk;
2054 }
2055 //______________________________________________________________________________
2056 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
2057 {
2058   //
2059   // Compute PID probabilities for the TRD
2060   //
2061   
2062   // set flat distribution (no decision)
2063   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2064   
2065   const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
2066   if (pidStatus!=kDetPidOk) return pidStatus;
2067
2068   UInt_t TRDslicesForPID[2];
2069   SetTRDSlices(TRDslicesForPID,PIDmethod);
2070   
2071   Float_t mom[6]={0.};
2072   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
2073   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
2074   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
2075   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
2076     mom[ilayer] = track->GetTRDmomentum(ilayer);
2077     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
2078       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
2079     }
2080   }
2081   
2082   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
2083   return kDetPidOk;
2084 }
2085
2086 //______________________________________________________________________________
2087 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2088 {
2089   //
2090   // Compute PID response for the EMCAL
2091   //
2092   
2093   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2094
2095   const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
2096   if (pidStatus!=kDetPidOk) return pidStatus;
2097
2098   const Int_t nMatchClus = track->GetEMCALcluster();
2099   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2100   
2101   const Double_t mom    = track->P();
2102   const Double_t pt     = track->Pt();
2103   const Int_t    charge = track->Charge();
2104   const Double_t fClsE  = matchedClus->E();
2105   const Double_t EovP   = fClsE/mom;
2106   
2107   // compute the probabilities
2108   fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p);
2109   return kDetPidOk;
2110 }
2111
2112 //______________________________________________________________________________
2113 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
2114 {
2115   //
2116   // Compute PID response for the PHOS
2117   //
2118   
2119   // set flat distribution (no decision)
2120   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2121   return kDetNoSignal;
2122 }
2123
2124 //______________________________________________________________________________
2125 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2126 {
2127   //
2128   // Compute PID response for the HMPID
2129   //
2130   
2131   // set flat distribution (no decision)
2132   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2133   
2134   const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
2135   if (pidStatus!=kDetPidOk) return pidStatus;
2136   
2137   fHMPIDResponse.GetProbability(track,nSpecies,p);
2138     
2139   return kDetPidOk;
2140 }
2141
2142 //______________________________________________________________________________
2143 AliPIDResponse::EDetPidStatus AliPIDResponse::GetITSPIDStatus(const AliVTrack *track) const
2144 {
2145   // compute ITS pid status
2146
2147   // check status bits
2148   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
2149     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
2150
2151   const Float_t dEdx=track->GetITSsignal();
2152   if (dEdx<=0) return kDetNoSignal;
2153   
2154   // requite at least 3 pid clusters
2155   const UChar_t clumap=track->GetITSClusterMap();
2156   Int_t nPointsForPid=0;
2157   for(Int_t i=2; i<6; i++){
2158     if(clumap&(1<<i)) ++nPointsForPid;
2159   }
2160   
2161   if(nPointsForPid<3) { 
2162     return kDetNoSignal;
2163   }
2164   
2165   return kDetPidOk;
2166 }
2167
2168 //______________________________________________________________________________
2169 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetTPCPIDStatus(const AliVTrack *track) const
2170 {
2171   // compute TPC pid status
2172   
2173   // check quality of the track
2174   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
2175
2176   // check pid values
2177   const Double_t dedx=track->GetTPCsignal();
2178   const UShort_t signalN=track->GetTPCsignalN();
2179   if (signalN<10 || dedx<10) return kDetNoSignal;
2180
2181   if (!(fArrPidResponseMaster && fArrPidResponseMaster->At(AliPID::kPion))) return kDetNoParams;
2182   
2183   return kDetPidOk;
2184 }
2185
2186 //______________________________________________________________________________
2187 AliPIDResponse::EDetPidStatus AliPIDResponse::GetTRDPIDStatus(const AliVTrack *track) const
2188 {
2189   // compute TRD pid status
2190
2191   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
2192   return kDetPidOk;
2193 }
2194
2195 //______________________________________________________________________________
2196 AliPIDResponse::EDetPidStatus AliPIDResponse::GetTOFPIDStatus(const AliVTrack *track) const
2197 {
2198   // compute TOF pid status
2199
2200   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
2201   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
2202
2203   return kDetPidOk;
2204 }
2205
2206 //______________________________________________________________________________
2207 Float_t AliPIDResponse::GetTOFMismatchProbability(const AliVTrack *track) const
2208 {
2209   // compute mismatch probability cross-checking at 5 sigmas with TPC
2210   // currently just implemented as a 5 sigma compatibility cut
2211
2212   // check pid status
2213   const EDetPidStatus tofStatus=GetTOFPIDStatus(track);
2214   if (tofStatus!=kDetPidOk) return 0.;
2215
2216   //mismatch
2217   const EDetPidStatus tpcStatus=GetTPCPIDStatus(track);
2218   if (tpcStatus!=kDetPidOk) return 0.;
2219   
2220   const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2221   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
2222   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
2223     AliPID::EParticleType type=AliPID::EParticleType(j);
2224     const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
2225     
2226     if (TMath::Abs(nsigmas)<5.){
2227       const Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
2228       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
2229     }
2230   }
2231   
2232   if (mismatch){
2233     return 1.;
2234   }
2235   
2236   return 0.;
2237 }
2238
2239 //______________________________________________________________________________
2240 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetHMPIDPIDStatus(const AliVTrack *track) const
2241 {
2242   // compute HMPID pid status
2243   
2244   Int_t ch = track->GetHMPIDcluIdx()/1000000;
2245   Double_t HMPIDsignal = track->GetHMPIDsignal(); 
2246   
2247   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0 || ch<0 || ch>6 || HMPIDsignal<0) return kDetNoSignal;
2248   
2249   return kDetPidOk;
2250 }
2251
2252 //______________________________________________________________________________
2253 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetPHOSPIDStatus(const AliVTrack */*track*/) const
2254 {
2255   // compute PHOS pid status
2256   return kDetNoSignal;  
2257 }
2258
2259 //______________________________________________________________________________
2260 AliPIDResponse::EDetPidStatus AliPIDResponse:: GetEMCALPIDStatus(const AliVTrack *track) const
2261 {
2262   // compute EMCAL pid status
2263
2264
2265   // Track matching
2266   const Int_t nMatchClus = track->GetEMCALcluster();
2267   if (nMatchClus<0) return kDetNoSignal;
2268
2269   AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2270
2271   if (!(matchedClus && matchedClus->IsEMCAL())) return kDetNoSignal;
2272
2273   const Int_t charge = track->Charge();
2274   if (TMath::Abs(charge)!=1) return kDetNoSignal;
2275
2276   if (!(fEMCALPIDParams && fEMCALPIDParams->At(AliPID::kElectron))) return kDetNoParams;
2277   
2278   return kDetPidOk;
2279
2280 }
2281
2282 //______________________________________________________________________________
2283 AliPIDResponse::EDetPidStatus AliPIDResponse::GetPIDStatus(EDetector detector, const AliVTrack *track) const
2284 {
2285   //
2286   // check pid status for a track
2287   //
2288
2289   switch (detector){
2290     case kITS:   return GetITSPIDStatus(track);   break;
2291     case kTPC:   return GetTPCPIDStatus(track);   break;
2292     case kTRD:   return GetTRDPIDStatus(track);   break;
2293     case kTOF:   return GetTOFPIDStatus(track);   break;
2294     case kPHOS:  return GetPHOSPIDStatus(track);  break;
2295     case kEMCAL: return GetEMCALPIDStatus(track); break;
2296     case kHMPID: return GetHMPIDPIDStatus(track); break;
2297     default: return kDetNoSignal;
2298   }
2299   return kDetNoSignal;
2300   
2301 }