TPC PID update:
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDResponse.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */
17
18 //-----------------------------------------------------------------
19 //        Base class for handling the pid response               //
20 //        functions of all detectors                             //
21 //        and give access to the nsigmas                         //
22 //                                                               //
23 //   Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch //
24 //-----------------------------------------------------------------
25
26 #include <TList.h>
27 #include <TObjArray.h>
28 #include <TPRegexp.h>
29 #include <TF1.h>
30 #include <TH2D.h>
31 #include <TSpline.h>
32 #include <TFile.h>
33 #include <TArrayI.h>
34 #include <TArrayF.h>
35 #include <TLinearFitter.h>
36
37 #include <AliVEvent.h>
38 #include <AliVTrack.h>
39 #include <AliLog.h>
40 #include <AliPID.h>
41 #include <AliOADBContainer.h>
42 #include <AliTRDPIDResponseObject.h>
43 #include <AliTOFPIDParams.h>
44
45 #include "AliPIDResponse.h"
46 #include "AliDetectorPID.h"
47
48 #include "AliCentrality.h"
49
50 ClassImp(AliPIDResponse);
51
52 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
53 TNamed("PIDResponse","PIDResponse"),
54 fITSResponse(isMC),
55 fTPCResponse(),
56 fTRDResponse(),
57 fTOFResponse(),
58 fEMCALResponse(),
59 fRange(5.),
60 fITSPIDmethod(kITSTruncMean),
61 fIsMC(isMC),
62 fCachePID(kTRUE),
63 fOADBPath(),
64 fCustomTPCpidResponse(),
65 fBeamType("PP"),
66 fLHCperiod(),
67 fMCperiodTPC(),
68 fMCperiodUser(),
69 fCurrentFile(),
70 fRecoPass(0),
71 fRecoPassUser(-1),
72 fRun(0),
73 fOldRun(0),
74 fResT0A(75.),
75 fResT0C(65.),
76 fResT0AC(55.),
77 fArrPidResponseMaster(NULL),
78 fResolutionCorrection(NULL),
79 fOADBvoltageMaps(NULL),
80 fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
81 fTRDPIDResponseObject(NULL),
82 fTOFtail(1.1),
83 fTOFPIDParams(NULL),
84 fEMCALPIDParams(NULL),
85 fCurrentEvent(NULL),
86 fCurrCentrality(0.0),
87 fTuneMConData(kFALSE)
88 {
89   //
90   // default ctor
91   //
92   AliLog::SetClassDebugLevel("AliPIDResponse",0);
93   AliLog::SetClassDebugLevel("AliESDpid",0);
94   AliLog::SetClassDebugLevel("AliAODpidUtil",0);
95
96 }
97
98 //______________________________________________________________________________
99 AliPIDResponse::~AliPIDResponse()
100 {
101   //
102   // dtor
103   //
104   delete fArrPidResponseMaster;
105   delete fTRDPIDResponseObject;
106   delete fTOFPIDParams;
107 }
108
109 //______________________________________________________________________________
110 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
111 TNamed(other),
112 fITSResponse(other.fITSResponse),
113 fTPCResponse(other.fTPCResponse),
114 fTRDResponse(other.fTRDResponse),
115 fTOFResponse(other.fTOFResponse),
116 fEMCALResponse(other.fEMCALResponse),
117 fRange(other.fRange),
118 fITSPIDmethod(other.fITSPIDmethod),
119 fIsMC(other.fIsMC),
120 fCachePID(other.fCachePID),
121 fOADBPath(other.fOADBPath),
122 fCustomTPCpidResponse(other.fCustomTPCpidResponse),
123 fBeamType("PP"),
124 fLHCperiod(),
125 fMCperiodTPC(),
126 fMCperiodUser(other.fMCperiodUser),
127 fCurrentFile(),
128 fRecoPass(0),
129 fRecoPassUser(other.fRecoPassUser),
130 fRun(0),
131 fOldRun(0),
132 fResT0A(75.),
133 fResT0C(65.),
134 fResT0AC(55.),
135 fArrPidResponseMaster(NULL),
136 fResolutionCorrection(NULL),
137 fOADBvoltageMaps(NULL),
138 fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
139 fTRDPIDResponseObject(NULL),
140 fTOFtail(1.1),
141 fTOFPIDParams(NULL),
142 fEMCALPIDParams(NULL),
143 fCurrentEvent(NULL),
144 fCurrCentrality(0.0),
145 fTuneMConData(kFALSE)
146 {
147   //
148   // copy ctor
149   //
150 }
151
152 //______________________________________________________________________________
153 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
154 {
155   //
156   // copy ctor
157   //
158   if(this!=&other) {
159     delete fArrPidResponseMaster;
160     TNamed::operator=(other);
161     fITSResponse=other.fITSResponse;
162     fTPCResponse=other.fTPCResponse;
163     fTRDResponse=other.fTRDResponse;
164     fTOFResponse=other.fTOFResponse;
165     fEMCALResponse=other.fEMCALResponse;
166     fRange=other.fRange;
167     fITSPIDmethod=other.fITSPIDmethod;
168     fOADBPath=other.fOADBPath;
169     fCustomTPCpidResponse=other.fCustomTPCpidResponse;
170     fIsMC=other.fIsMC;
171     fCachePID=other.fCachePID;
172     fBeamType="PP";
173     fLHCperiod="";
174     fMCperiodTPC="";
175     fMCperiodUser=other.fMCperiodUser;
176     fCurrentFile="";
177     fRecoPass=0;
178     fRecoPassUser=other.fRecoPassUser;
179     fRun=0;
180     fOldRun=0;
181     fResT0A=75.;
182     fResT0C=65.;
183     fResT0AC=55.;
184     fArrPidResponseMaster=NULL;
185     fResolutionCorrection=NULL;
186     fOADBvoltageMaps=NULL;
187         fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
188     fTRDPIDResponseObject=NULL;
189     fEMCALPIDParams=NULL;
190     fTOFtail=1.1;
191     fTOFPIDParams=NULL;
192     fCurrentEvent=other.fCurrentEvent;
193
194   }
195   return *this;
196 }
197
198 //______________________________________________________________________________
199 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
200 {
201   //
202   // NumberOfSigmas for 'detCode'
203   //
204
205   switch (detCode){
206     case kDetITS: return NumberOfSigmasITS(track, type); break;
207     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
208     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
209     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
210     default: return -999.;
211   }
212
213 }
214
215 //______________________________________________________________________________
216 Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
217 {
218   //
219   // NumberOfSigmas for 'detCode'
220   //
221   return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
222 }
223
224 //______________________________________________________________________________
225 // public buffered versions of the PID calculation
226 //
227
228 //______________________________________________________________________________
229 Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
230 {
231   //
232   // Calculate the number of sigmas in the ITS
233   //
234   
235   AliVTrack *track=(AliVTrack*)vtrack;
236   
237   // look for cached value first
238   const AliDetectorPID *detPID=track->GetDetectorPID();
239   const EDetector detector=kITS;
240
241   if ( detPID && detPID->HasNumberOfSigmas(detector)){
242     return detPID->GetNumberOfSigmas(detector, type);
243   } else if (fCachePID) {
244     FillTrackDetectorPID(track, detector);
245     detPID=track->GetDetectorPID();
246     return detPID->GetNumberOfSigmas(detector, type);
247   }
248
249   return GetNumberOfSigmasITS(track, type);
250 }
251
252 //______________________________________________________________________________
253 Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
254 {
255   //
256   // Calculate the number of sigmas in the TPC
257   //
258   
259   AliVTrack *track=(AliVTrack*)vtrack;
260   
261   // look for cached value first
262   const AliDetectorPID *detPID=track->GetDetectorPID();
263   const EDetector detector=kTPC;
264   
265   if ( detPID && detPID->HasNumberOfSigmas(detector)){
266     return detPID->GetNumberOfSigmas(detector, type);
267   } else if (fCachePID) {
268     FillTrackDetectorPID(track, detector);
269     detPID=track->GetDetectorPID();
270     return detPID->GetNumberOfSigmas(detector, type);
271   }
272   
273   return GetNumberOfSigmasTPC(track, type);
274 }
275
276 //______________________________________________________________________________
277 Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack, 
278                                            AliPID::EParticleType type,
279                                            AliTPCPIDResponse::ETPCdEdxSource dedxSource) const
280 {
281   //get number of sigmas according the selected TPC gain configuration scenario
282   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
283
284   Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection);
285
286   return nSigma;
287 }
288
289 //______________________________________________________________________________
290 Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
291 {
292   //
293   // Calculate the number of sigmas in the TOF
294   //
295   
296   AliVTrack *track=(AliVTrack*)vtrack;
297   
298   // look for cached value first
299   const AliDetectorPID *detPID=track->GetDetectorPID();
300   const EDetector detector=kTOF;
301   
302   if ( detPID && detPID->HasNumberOfSigmas(detector)){
303     return detPID->GetNumberOfSigmas(detector, type);
304   } else if (fCachePID) {
305     FillTrackDetectorPID(track, detector);
306     detPID=track->GetDetectorPID();
307     return detPID->GetNumberOfSigmas(detector, type);
308   }
309   
310   return GetNumberOfSigmasTOF(track, type);
311 }
312
313 //______________________________________________________________________________
314 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
315 {
316   //
317   // Calculate the number of sigmas in the EMCAL
318   //
319   
320   AliVTrack *track=(AliVTrack*)vtrack;
321
322   // look for cached value first
323   const AliDetectorPID *detPID=track->GetDetectorPID();
324   const EDetector detector=kEMCAL;
325   
326   if ( detPID && detPID->HasNumberOfSigmas(detector)){
327     return detPID->GetNumberOfSigmas(detector, type);
328   } else if (fCachePID) {
329     FillTrackDetectorPID(track, detector);
330     detPID=track->GetDetectorPID();
331     return detPID->GetNumberOfSigmas(detector, type);
332   }
333   
334   return GetNumberOfSigmasEMCAL(track, type);
335 }
336
337 //______________________________________________________________________________
338 Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4])  const
339 {
340   //
341   // emcal nsigma with eop and showershape
342   //
343   AliVTrack *track=(AliVTrack*)vtrack;
344   
345   AliVCluster *matchedClus = NULL;
346
347   Double_t mom     = -1.; 
348   Double_t pt      = -1.; 
349   Double_t EovP    = -1.;
350   Double_t fClsE   = -1.;
351
352   // initialize eop and shower shape parameters
353   eop = -1.;
354   for(Int_t i = 0; i < 4; i++){
355     showershape[i] = -1.;
356   }
357   
358   Int_t nMatchClus = -1;
359   Int_t charge     = 0;
360   
361   // Track matching
362   nMatchClus = track->GetEMCALcluster();
363   if(nMatchClus > -1){
364
365     mom    = track->P();
366     pt     = track->Pt();
367     charge = track->Charge();
368     
369     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
370     
371     if(matchedClus){
372       
373       // matched cluster is EMCAL
374       if(matchedClus->IsEMCAL()){
375         
376         fClsE       = matchedClus->E();
377         EovP        = fClsE/mom;
378         
379         // fill used EMCAL variables here
380         eop            = EovP; // E/p
381         showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
382         showershape[1] = matchedClus->GetM02(); // long axis
383         showershape[2] = matchedClus->GetM20(); // short axis
384         showershape[3] = matchedClus->GetDispersion(); // dispersion
385
386         // look for cached value first
387         const AliDetectorPID *detPID=track->GetDetectorPID();
388         const EDetector detector=kEMCAL;
389         
390         if ( detPID && detPID->HasNumberOfSigmas(detector)){
391           return detPID->GetNumberOfSigmas(detector, type);
392         } else if (fCachePID) {
393           FillTrackDetectorPID(track, detector);
394           detPID=track->GetDetectorPID();
395           return detPID->GetNumberOfSigmas(detector, type);
396         }
397         
398         // NSigma value really meaningful only for electrons!
399         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
400       }
401     }
402   }
403   return -999;
404 }
405
406
407 //______________________________________________________________________________
408 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
409 {
410   //
411   // Compute PID response of 'detCode'
412   //
413
414   switch (detCode){
415     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
416     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
417     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
418     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
419     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
420     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
421     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
422     default: return kDetNoSignal;
423   }
424 }
425
426 //______________________________________________________________________________
427 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
428 {
429   //
430   // Compute PID response of 'detCode'
431   //
432
433   return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
434 }
435
436 //______________________________________________________________________________
437 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
438 {
439   //
440   // Compute PID response for the ITS
441   //
442
443   // look for cached value first
444   const AliDetectorPID *detPID=track->GetDetectorPID();
445   const EDetector detector=kITS;
446   
447   if ( detPID && detPID->HasRawProbabilitiy(detector)){
448     return detPID->GetRawProbability(detector, p, nSpecies);
449   } else if (fCachePID) {
450     FillTrackDetectorPID(track, detector);
451     detPID=track->GetDetectorPID();
452     return detPID->GetRawProbability(detector, p, nSpecies);
453   }
454   
455   return GetComputeITSProbability(track, nSpecies, p);
456 }
457 //______________________________________________________________________________
458 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
459 {
460   //
461   // Compute PID response for the TPC
462   //
463   
464   // look for cached value first
465   const AliDetectorPID *detPID=track->GetDetectorPID();
466   const EDetector detector=kTPC;
467   
468   if ( detPID && detPID->HasRawProbabilitiy(detector)){
469     return detPID->GetRawProbability(detector, p, nSpecies);
470   } else if (fCachePID) {
471     FillTrackDetectorPID(track, detector);
472     detPID=track->GetDetectorPID();
473     return detPID->GetRawProbability(detector, p, nSpecies);
474   }
475   
476   return GetComputeTPCProbability(track, nSpecies, p);
477 }
478 //______________________________________________________________________________
479 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
480 {
481   //
482   // Compute PID response for the
483   //
484   
485   const AliDetectorPID *detPID=track->GetDetectorPID();
486   const EDetector detector=kTOF;
487   
488   if ( detPID && detPID->HasRawProbabilitiy(detector)){
489     return detPID->GetRawProbability(detector, p, nSpecies);
490   } else if (fCachePID) {
491     FillTrackDetectorPID(track, detector);
492     detPID=track->GetDetectorPID();
493     return detPID->GetRawProbability(detector, p, nSpecies);
494   }
495   
496   return GetComputeTOFProbability(track, nSpecies, p);
497 }
498 //______________________________________________________________________________
499 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
500 {
501   //
502   // Compute PID response for the
503   //
504   
505   const AliDetectorPID *detPID=track->GetDetectorPID();
506   const EDetector detector=kTRD;
507
508   // chacke only for the default method (1d at the moment)
509   if (PIDmethod!=AliTRDPIDResponse::kLQ1D) return GetComputeTRDProbability(track, nSpecies, p, PIDmethod);
510   
511   if ( detPID && detPID->HasRawProbabilitiy(detector)){
512     return detPID->GetRawProbability(detector, p, nSpecies);
513   } else if (fCachePID) {
514     FillTrackDetectorPID(track, detector);
515     detPID=track->GetDetectorPID();
516     return detPID->GetRawProbability(detector, p, nSpecies);
517   }
518   
519   return GetComputeTRDProbability(track, nSpecies, p);
520 }
521 //______________________________________________________________________________
522 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
523 {
524   //
525   // Compute PID response for the EMCAL
526   //
527   
528   const AliDetectorPID *detPID=track->GetDetectorPID();
529   const EDetector detector=kEMCAL;
530   
531   if ( detPID && detPID->HasRawProbabilitiy(detector)){
532     return detPID->GetRawProbability(detector, p, nSpecies);
533   } else if (fCachePID) {
534     FillTrackDetectorPID(track, detector);
535     detPID=track->GetDetectorPID();
536     return detPID->GetRawProbability(detector, p, nSpecies);
537   }
538   
539   return GetComputeEMCALProbability(track, nSpecies, p);
540 }
541 //______________________________________________________________________________
542 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
543 {
544   //
545   // Compute PID response for the PHOS
546   //
547   
548   // look for cached value first
549 //   if (track->GetDetectorPID()){
550 //     return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies);
551 //   }
552   
553   // set flat distribution (no decision)
554   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
555   return kDetNoSignal;
556 }
557 //______________________________________________________________________________
558 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
559 {
560   //
561   // Compute PID response for the HMPID
562   //
563
564   const AliDetectorPID *detPID=track->GetDetectorPID();
565   const EDetector detector=kHMPID;
566
567   if ( detPID && detPID->HasRawProbabilitiy(detector)){
568     return detPID->GetRawProbability(detector, p, nSpecies);
569   } else if (fCachePID) {
570     FillTrackDetectorPID(track, detector);
571     detPID=track->GetDetectorPID();
572     return detPID->GetRawProbability(detector, p, nSpecies);
573   }
574
575   return GetComputeHMPIDProbability(track, nSpecies, p);
576 }
577
578 //______________________________________________________________________________
579 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
580 {
581   //
582   // Apply settings for the current event
583   //
584   fRecoPass=pass;
585   
586
587   fCurrentEvent=NULL;
588   if (!event) return;
589   fCurrentEvent=event;
590   if (run>0) fRun=run;
591   else fRun=event->GetRunNumber();
592   
593   if (fRun!=fOldRun){
594     ExecNewRun();
595     fOldRun=fRun;
596   }
597   
598   //TPC resolution parametrisation PbPb
599   if ( fResolutionCorrection ){
600     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
601     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
602   }
603   
604   //TOF resolution
605   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
606
607
608   // Get and set centrality
609   AliCentrality *centrality = event->GetCentrality();
610   if(centrality){
611     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
612   }
613   else{
614     fCurrCentrality = -1;
615   }
616 }
617
618 //______________________________________________________________________________
619 void AliPIDResponse::ExecNewRun()
620 {
621   //
622   // Things to Execute upon a new run
623   //
624   SetRecoInfo();
625   
626   SetITSParametrisation();
627   
628   SetTPCPidResponseMaster();
629   SetTPCParametrisation();
630   SetTPCEtaMaps();
631
632   SetTRDPidResponseMaster(); 
633   InitializeTRDResponse();
634
635   SetEMCALPidResponseMaster(); 
636   InitializeEMCALResponse();
637   
638   SetTOFPidResponseMaster();
639   InitializeTOFResponse();
640
641   if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
642 }
643
644 //______________________________________________________________________________
645 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
646 {
647   //
648   // Get TPC multiplicity in bins of 150
649   //
650   
651   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
652   Double_t tpcMulti=0.;
653   if(vertexTPC){
654     Double_t vertexContribTPC=vertexTPC->GetNContributors();
655     tpcMulti=vertexContribTPC/150.;
656     if (tpcMulti>20.) tpcMulti=20.;
657   }
658   
659   return tpcMulti;
660 }
661
662 //______________________________________________________________________________
663 void AliPIDResponse::SetRecoInfo()
664 {
665   //
666   // Set reconstruction information
667   //
668   
669   //reset information
670   fLHCperiod="";
671   fMCperiodTPC="";
672   
673   fBeamType="";
674     
675   fBeamType="PP";
676   
677   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
678   TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
679
680   //find the period by run number (UGLY, but not stored in ESD and AOD... )
681   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
682   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
683   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
684   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
685   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
686   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
687   else if (fRun>=136851&&fRun<=139846) {
688     fLHCperiod="LHC10H";
689     fMCperiodTPC="LHC10H8";
690     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
691     fBeamType="PBPB";
692   }
693   else if (fRun>=139847&&fRun<=146974) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
694   //TODO: periods 11B (146975-150721), 11C (150722-155837) are not yet treated assume 11d for the moment
695   else if (fRun>=146975&&fRun<=155837) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
696   else if (fRun>=155838&&fRun<=159649) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
697   // also for 11e (159650-162750),f(162751-165771) use 11d
698   else if (fRun>=159650&&fRun<=162750) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
699   else if (fRun>=162751&&fRun<=165771) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
700   
701   else if (fRun>=165772 && fRun<=170718) {
702     fLHCperiod="LHC11H";
703     fMCperiodTPC="LHC11A10";
704     fBeamType="PBPB";
705     if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
706   }
707   if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
708   // for the moment use LHC12b parameters up to LHC12e
709   if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
710 //   if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
711 //   if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
712 //   if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
713
714 //   if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
715 //   if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
716 //   if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
717 // for the moment use 12g parametrisation for all full gain runs (LHC12f+)
718   if (fRun >= 186636  ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
719
720   //exception new pp MC productions from 2011
721   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
722   // exception for 11f1
723   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
724 }
725
726 //______________________________________________________________________________
727 void AliPIDResponse::SetITSParametrisation()
728 {
729   //
730   // Set the ITS parametrisation
731   //
732 }
733
734  
735 //______________________________________________________________________________
736 void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
737 {
738   if (h->GetBinContent(binX, binY) <= 1e-4)
739     return; // Reject bins without content (within some numerical precision) or with strange content
740     
741   Double_t coord[2] = {0, 0};
742   coord[0] = h->GetXaxis()->GetBinCenter(binX);
743   coord[1] = h->GetYaxis()->GetBinCenter(binY);
744   Double_t binError = h->GetBinError(binX, binY);
745   if (binError <= 0) {
746     binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
747     printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
748   }
749   linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
750 }
751
752
753 //______________________________________________________________________________
754 TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
755 {
756   if (!h)
757     return 0x0;
758   
759   // Interpolate to finer map
760   TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
761   
762   Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
763   Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
764   Int_t nBinsX = 20;
765   // 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,
766   // scale the number of bins correspondingly
767   Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - lowerMapBoundY) * 40);
768   Int_t nBinsXrefined = nBinsX * refineFactorX;
769   Int_t nBinsYrefined = nBinsY * refineFactorY; 
770   
771   TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()),  Form("%s (refined)", h->GetTitle()),
772                             nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
773                             nBinsYrefined, lowerMapBoundY, upperMapBoundY);
774   
775   for (Int_t binX = 1; binX <= nBinsXrefined; binX++)  {
776     for (Int_t binY = 1; binY <= nBinsYrefined; binY++)  {
777       
778       hRefined->SetBinContent(binX, binY, 1); // Default value is 1
779       
780       Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
781       Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
782       
783       /*TODO NOW NOW
784       linExtrapolation->ClearPoints();
785       
786       // For interpolation: Just take the corresponding bin from the old histo.
787       // For extrapolation: take the last available bin from the old histo.
788       // If the boundaries are to be skipped, also skip the corresponding bins
789       Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
790       if (oldBinX < 1)  
791         oldBinX = 1;
792       if (oldBinX > nBinsX)
793         oldBinX = nBinsX;
794       
795       Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
796       if (oldBinY < 1)  
797         oldBinY = 1;
798       if (oldBinY > nBinsY)
799         oldBinY = nBinsY;
800       
801       // Neighbours left column
802       if (oldBinX >= 2) {
803         if (oldBinY >= 2) {
804           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
805         }
806         
807         AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
808         
809         if (oldBinY < nBinsY) {
810           AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
811         }
812       }
813       
814       // Neighbours (and point itself) same column
815       if (oldBinY >= 2) {
816         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
817       }
818         
819       AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
820         
821       if (oldBinY < nBinsY) {
822         AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
823       }
824       
825       // Neighbours right column
826       if (oldBinX < nBinsX) {
827         if (oldBinY >= 2) {
828           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
829         }
830         
831         AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
832         
833         if (oldBinY < nBinsY) {
834           AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
835         }
836       }
837       
838       
839       // Fit 2D-hyperplane
840       if (linExtrapolation->GetNpoints() <= 0)
841         continue;
842         
843       if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
844         continue;
845       
846       // Fill the bin of the refined histogram with the extrapolated value
847       Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
848                                  + linExtrapolation->GetParameter(2) * centerY;
849       */
850       Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;//TODO NOW NOW NOW
851       hRefined->SetBinContent(binX, binY, interpolatedValue);      
852     }
853   } 
854   
855   delete linExtrapolation;
856   
857   return hRefined;
858 }
859
860 //______________________________________________________________________________
861 void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
862                                    Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
863 {
864   //
865   // Load the TPC eta correction maps from the OADB
866   //
867   
868   TString dataType = "DATA";
869   TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
870   
871   if (fIsMC)  {
872     dataType = "MC";
873     fRecoPass = 1;
874     
875     if (fMCperiodTPC.IsNull()) {
876       AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
877       return;
878     }
879   
880     period = fMCperiodTPC;
881   }
882   
883   TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), fRecoPass);
884   
885   AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), fRecoPass));
886   if (fUseTPCEtaCorrection == kFALSE) {
887     // Disable eta correction via setting no maps
888     if (!fTPCResponse.SetEtaCorrMap(0x0))
889       AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled"); 
890     else
891       AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
892     
893     if (!fTPCResponse.SetSigmaParams(0x0, 0))
894       AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma"); 
895     else
896       AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
897     
898     return;
899   }
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(), fRecoPass)); 
907   
908   Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
909                                               Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass));
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) {
919       TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
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(), fRecoPass)); 
957   
958   statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
959                                              Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass));
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) {
969       TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
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   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1100   Bool_t found=kFALSE;
1101   //
1102   //set the new PID splines
1103   //
1104   if (fArrPidResponseMaster){
1105     Int_t recopass = fRecoPass;
1106     if(fTuneMConData) recopass = fRecoPassUser;
1107     //for MC don't use period information
1108     //if (fIsMC) period="[A-Z0-9]*";
1109     //for MC use MC period information
1110     //pattern for the default entry (valid for all particles)
1111     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1112
1113     //find particle id ang gain scenario
1114     for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
1115     {
1116       TObject *grAll=NULL;
1117       TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
1118       gainScenario.ToUpper();
1119       //loop over entries and filter them
1120       for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
1121       {
1122         TObject *responseFunction=fArrPidResponseMaster->At(iresp);
1123         if (responseFunction==NULL) continue;
1124         TString responseName=responseFunction->GetName();
1125          
1126         if (!reg.MatchB(responseName)) continue;
1127
1128         TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
1129         TObject* tmp=NULL;
1130         tmp=arr->At(1); if (!tmp) continue;
1131         TString particleName=tmp->GetName();
1132         tmp=arr->At(3); if (!tmp) continue;
1133         TString gainScenarioName=tmp->GetName();
1134         delete arr;
1135         if (particleName.IsNull()) continue;
1136         if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
1137         else 
1138         {
1139           for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1140           {
1141             TString particle=AliPID::ParticleName(ispec);
1142             particle.ToUpper();
1143             //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
1144             if ( particle == particleName && gainScenario == gainScenarioName )
1145             {
1146               fTPCResponse.SetResponseFunction( responseFunction,
1147                                                 (AliPID::EParticleType)ispec,
1148                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1149               fTPCResponse.SetUseDatabase(kTRUE);
1150               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
1151               found=kTRUE;
1152               // overwrite default with proton spline (for light nuclei)
1153               if (ispec==AliPID::kProton) grAll=responseFunction;
1154               break;
1155             }
1156           }
1157         }
1158       }
1159       if (grAll)
1160       {
1161         for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1162         {
1163           if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
1164                                                  (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
1165           {
1166               fTPCResponse.SetResponseFunction( grAll,
1167                                                 (AliPID::EParticleType)ispec,
1168                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1169               fTPCResponse.SetUseDatabase(kTRUE);
1170               AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
1171               found=kTRUE;
1172           }
1173         }
1174       }
1175     }
1176   }
1177   else AliInfo("no fArrPidResponseMaster");
1178
1179   if (!found){
1180     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1181   }
1182
1183   //
1184   // Setup resolution parametrisation
1185   //
1186   
1187   //default
1188   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1189   
1190   if (fRun>=122195){
1191     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1192   }
1193
1194   if (fRun>=186636){
1195 //   if (fRun>=188356){
1196     fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1197   }
1198   
1199   if (fArrPidResponseMaster)
1200   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
1201   
1202   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
1203
1204   //read in the voltage map
1205   TVectorF* gsm = 0x0;
1206   if (fOADBvoltageMaps) gsm=dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
1207   if (gsm) 
1208   {
1209     fTPCResponse.SetVoltageMap(*gsm);
1210     TString vals;
1211     AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1212     vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1213     AliInfo(vals.Data());
1214     vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1215     AliInfo(vals.Data());
1216     vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1217     AliInfo(vals.Data());
1218     vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1219     AliInfo(vals.Data());
1220   }
1221   else AliInfo("no voltage map, ideal default assumed");
1222 }
1223
1224 //______________________________________________________________________________
1225 void AliPIDResponse::SetTRDPidResponseMaster()
1226 {
1227   //
1228   // Load the TRD pid params and references from the OADB
1229   //
1230   if(fTRDPIDResponseObject) return;
1231   AliOADBContainer contParams("contParams"); 
1232
1233   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1234   if(statusResponse){
1235     AliError("Failed initializing PID Response Object from OADB");
1236   } else {
1237     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1238     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1239     if(!fTRDPIDResponseObject){
1240       AliError(Form("TRD Response not found in run %d", fRun));
1241     }
1242   }
1243 }
1244
1245 //______________________________________________________________________________
1246 void AliPIDResponse::InitializeTRDResponse(){
1247   //
1248   // Set PID Params and references to the TRD PID response
1249   // 
1250     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1251 }
1252
1253 //______________________________________________________________________________
1254 void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1255
1256     if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1257         // backward compatibility for setting with 8 slices
1258         TRDslicesForPID[0] = 0;
1259         TRDslicesForPID[1] = 7;
1260     }
1261     else{
1262         if(method==AliTRDPIDResponse::kLQ1D){
1263             TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1264             TRDslicesForPID[1] = 0;
1265         }
1266         if(method==AliTRDPIDResponse::kLQ2D){
1267             TRDslicesForPID[0] = 1;
1268             TRDslicesForPID[1] = 7;
1269         }
1270     }
1271     AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
1272 }
1273
1274 //______________________________________________________________________________
1275 void AliPIDResponse::SetTOFPidResponseMaster()
1276 {
1277   //
1278   // Load the TOF pid params from the OADB
1279   //
1280
1281   if (fTOFPIDParams) delete fTOFPIDParams;
1282   fTOFPIDParams=NULL;
1283
1284   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1285   if (oadbf && oadbf->IsOpen()) {
1286     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1287     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1288     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
1289     oadbf->Close();
1290     delete oadbc;
1291   }
1292   delete oadbf;
1293
1294   if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1295 }
1296
1297 //______________________________________________________________________________
1298 void AliPIDResponse::InitializeTOFResponse(){
1299   //
1300   // Set PID Params to the TOF PID response
1301   //
1302
1303   AliInfo("TOF PID Params loaded from OADB");
1304   AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1305   AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1306   AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1307                fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1308   
1309   for (Int_t i=0;i<4;i++) {
1310     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1311   }
1312   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1313
1314   AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1315   Float_t t0Spread[4];
1316   for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1317   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]));
1318   Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1319   Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1320   if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1321     fResT0AC=t0Spread[3];
1322     fResT0A=TMath::Sqrt(a);
1323     fResT0C=TMath::Sqrt(c);
1324   } else {
1325     AliInfo("  TZERO spreads not present or inconsistent, loading default");
1326     fResT0A=75.;
1327     fResT0C=65.;
1328     fResT0AC=55.;
1329   }
1330   AliInfo(Form("  TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1331
1332 }
1333
1334
1335 //______________________________________________________________________________
1336 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1337   //
1338   // Check whether track is identified as electron under a given electron efficiency hypothesis
1339     //
1340
1341   Double_t probs[AliPID::kSPECIES];
1342   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod);
1343
1344   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1345   // Take mean of the TRD momenta in the given tracklets
1346   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1347   Int_t nmomenta = 0;
1348   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1349     if(vtrack->GetTRDmomentum(iPl) > 0.){
1350       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
1351     }
1352   }
1353   p = TMath::Mean(nmomenta, trdmomenta);
1354
1355   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1356 }
1357
1358 //______________________________________________________________________________
1359 void AliPIDResponse::SetEMCALPidResponseMaster()
1360 {
1361   //
1362   // Load the EMCAL pid response functions from the OADB
1363   //
1364   TObjArray* fEMCALPIDParamsRun      = NULL;
1365   TObjArray* fEMCALPIDParamsPass     = NULL;
1366
1367   if(fEMCALPIDParams) return;
1368   AliOADBContainer contParams("contParams"); 
1369
1370   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1371   if(statusPars){
1372     AliError("Failed initializing PID Params from OADB");
1373   } 
1374   else {
1375     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1376
1377     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1378     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1379     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1380
1381     if(!fEMCALPIDParams){
1382       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1383       AliInfo("Will take the standard LHC11d instead ...");
1384
1385       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1386       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
1387       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1388
1389       if(!fEMCALPIDParams){
1390         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));     
1391       }
1392     }
1393   }
1394 }
1395
1396 //______________________________________________________________________________
1397 void AliPIDResponse::InitializeEMCALResponse(){
1398   //
1399   // Set PID Params to the EMCAL PID response
1400   // 
1401   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1402
1403 }
1404
1405 //______________________________________________________________________________
1406 void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
1407 {
1408   //
1409   // create detector PID information and setup the transient pointer in the track
1410   //
1411   
1412   // check if detector number is inside accepted range
1413   if (detector == kNdetectors) return;
1414   
1415   // get detector pid
1416   AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1417   if (!detPID) {
1418     detPID=new AliDetectorPID;
1419     (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1420   }
1421   
1422   //check if values exist
1423   if (detPID->HasRawProbabilitiy(detector) && detPID->HasNumberOfSigmas(detector)) return;
1424   
1425   //TODO: which particles to include? See also the loops below...
1426   Double_t values[AliPID::kSPECIESC]={0};
1427
1428   //nsigmas
1429   for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1430     values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
1431   detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC);
1432   
1433   //probabilities
1434   EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1435   detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1436 }
1437
1438 //______________________________________________________________________________
1439 void AliPIDResponse::FillTrackDetectorPID()
1440 {
1441   //
1442   // create detector PID information and setup the transient pointer in the track
1443   //
1444
1445   if (!fCurrentEvent) return;
1446   
1447   for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1448     AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1449     if (!track) continue;
1450
1451     for (Int_t idet=0; idet<kNdetectors; ++idet){
1452       FillTrackDetectorPID(track, (EDetector)idet);
1453     }
1454   }
1455 }
1456
1457 //______________________________________________________________________________
1458 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1459   //
1460   // Set TOF response function
1461   // Input option for event_time used
1462   //
1463   
1464     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1465     if(t0spread < 10) t0spread = 80;
1466
1467     // T0 from TOF algorithm
1468
1469     Bool_t flagT0TOF=kFALSE;
1470     Bool_t flagT0T0=kFALSE;
1471     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1472     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1473     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1474
1475     // T0-TOF arrays
1476     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1477     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1478     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1479       estimatedT0event[i]=0.0;
1480       estimatedT0resolution[i]=0.0;
1481       startTimeMask[i] = 0;
1482     }
1483
1484     Float_t resT0A=fResT0A;
1485     Float_t resT0C=fResT0C;
1486     Float_t resT0AC=fResT0AC;
1487     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1488         flagT0T0=kTRUE;
1489     }
1490
1491
1492     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1493
1494     if (tofHeader) { // read global info and T0-TOF
1495       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1496       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1497       if(t0spread < 10) t0spread = 80;
1498
1499       flagT0TOF=kTRUE;
1500       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1501         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1502         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1503         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1504       }
1505
1506       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1507       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1508       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1509       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1510         Int_t icurrent = (Int_t)ibin->GetAt(j);
1511         startTime[icurrent]=t0Bin->GetAt(j);
1512         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1513         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1514       }
1515     }
1516
1517     // for cut of 3 sigma on t0 spread
1518     Float_t t0cut = 3 * t0spread;
1519     if(t0cut < 500) t0cut = 500;
1520
1521     if(option == kFILL_T0){ // T0-FILL is used
1522         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1523           estimatedT0event[i]=0.0;
1524           estimatedT0resolution[i]=t0spread;
1525         }
1526         fTOFResponse.SetT0event(estimatedT0event);
1527         fTOFResponse.SetT0resolution(estimatedT0resolution);
1528     }
1529
1530     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1531         if(flagT0TOF){
1532             fTOFResponse.SetT0event(startTime);
1533             fTOFResponse.SetT0resolution(startTimeRes);
1534             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1535               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1536               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1537             }
1538         }
1539         else{
1540             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1541               estimatedT0event[i]=0.0;
1542               estimatedT0resolution[i]=t0spread;
1543               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1544             }
1545             fTOFResponse.SetT0event(estimatedT0event);
1546             fTOFResponse.SetT0resolution(estimatedT0resolution);
1547         }
1548     }
1549     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1550         Float_t t0AC=-10000;
1551         Float_t t0A=-10000;
1552         Float_t t0C=-10000;
1553         if(flagT0T0){
1554             t0A= vevent->GetT0TOF()[1];
1555             t0C= vevent->GetT0TOF()[2];
1556         //      t0AC= vevent->GetT0TOF()[0];
1557         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1558         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1559         t0AC /= resT0AC*resT0AC;
1560         }
1561
1562         Float_t t0t0Best = 0;
1563         Float_t t0t0BestRes = 9999;
1564         Int_t t0used=0;
1565         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1566             t0t0Best = t0AC;
1567             t0t0BestRes = resT0AC;
1568             t0used=6;
1569         }
1570         else if(TMath::Abs(t0C) < t0cut){
1571             t0t0Best = t0C;
1572             t0t0BestRes = resT0C;
1573             t0used=4;
1574         }
1575         else if(TMath::Abs(t0A) < t0cut){
1576             t0t0Best = t0A;
1577             t0t0BestRes = resT0A;
1578             t0used=2;
1579         }
1580
1581         if(flagT0TOF){ // if T0-TOF info is available
1582             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1583                 if(t0t0BestRes < 999){
1584                   if(startTimeRes[i] < t0spread){
1585                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1586                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1587                     estimatedT0event[i]=t0best / wtot;
1588                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1589                     startTimeMask[i] = t0used+1;
1590                   }
1591                   else {
1592                     estimatedT0event[i]=t0t0Best;
1593                     estimatedT0resolution[i]=t0t0BestRes;
1594                     startTimeMask[i] = t0used;
1595                   }
1596                 }
1597                 else{
1598                   estimatedT0event[i]=startTime[i];
1599                   estimatedT0resolution[i]=startTimeRes[i];
1600                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1601                 }
1602                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1603             }
1604             fTOFResponse.SetT0event(estimatedT0event);
1605             fTOFResponse.SetT0resolution(estimatedT0resolution);
1606         }
1607         else{ // if no T0-TOF info is available
1608             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1609               fTOFResponse.SetT0binMask(i,t0used);
1610               if(t0t0BestRes < 999){
1611                 estimatedT0event[i]=t0t0Best;
1612                 estimatedT0resolution[i]=t0t0BestRes;
1613               }
1614               else{
1615                 estimatedT0event[i]=0.0;
1616                 estimatedT0resolution[i]=t0spread;
1617               }
1618             }
1619             fTOFResponse.SetT0event(estimatedT0event);
1620             fTOFResponse.SetT0resolution(estimatedT0resolution);
1621         }
1622     }
1623
1624     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1625         Float_t t0AC=-10000;
1626         Float_t t0A=-10000;
1627         Float_t t0C=-10000;
1628         if(flagT0T0){
1629             t0A= vevent->GetT0TOF()[1];
1630             t0C= vevent->GetT0TOF()[2];
1631         //      t0AC= vevent->GetT0TOF()[0];
1632         t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
1633         resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1634         t0AC /= resT0AC*resT0AC;
1635         }
1636
1637         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1638             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1639               estimatedT0event[i]=t0AC;
1640               estimatedT0resolution[i]=resT0AC;
1641               fTOFResponse.SetT0binMask(i,6);
1642             }
1643         }
1644         else if(TMath::Abs(t0C) < t0cut){
1645             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1646               estimatedT0event[i]=t0C;
1647               estimatedT0resolution[i]=resT0C;
1648               fTOFResponse.SetT0binMask(i,4);
1649             }
1650         }
1651         else if(TMath::Abs(t0A) < t0cut){
1652             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1653               estimatedT0event[i]=t0A;
1654               estimatedT0resolution[i]=resT0A;
1655               fTOFResponse.SetT0binMask(i,2);
1656             }
1657         }
1658         else{
1659             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1660               estimatedT0event[i]=0.0;
1661               estimatedT0resolution[i]=t0spread;
1662               fTOFResponse.SetT0binMask(i,0);
1663             }
1664         }
1665         fTOFResponse.SetT0event(estimatedT0event);
1666         fTOFResponse.SetT0resolution(estimatedT0resolution);
1667     }
1668     delete [] startTime;
1669     delete [] startTimeRes;
1670     delete [] startTimeMask;
1671     delete [] estimatedT0event;
1672     delete [] estimatedT0resolution;
1673 }
1674
1675 //______________________________________________________________________________
1676 // private non cached versions of the PID calculation
1677 //
1678
1679
1680 //______________________________________________________________________________
1681 Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
1682 {
1683   //
1684   // NumberOfSigmas for 'detCode'
1685   //
1686   
1687   switch (detCode){
1688     case kITS:   return GetNumberOfSigmasITS(track, type); break;
1689     case kTPC:   return GetNumberOfSigmasTPC(track, type); break;
1690     case kTOF:   return GetNumberOfSigmasTOF(track, type); break;
1691     case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
1692     default: return -999.;
1693   }
1694   
1695 }
1696
1697
1698
1699 //______________________________________________________________________________
1700 Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
1701 {
1702   //
1703   // Calculate the number of sigmas in the ITS
1704   //
1705   
1706   AliVTrack *track=(AliVTrack*)vtrack;
1707   
1708   Float_t dEdx=track->GetITSsignal();
1709   if (dEdx<=0) return -999.;
1710   
1711   UChar_t clumap=track->GetITSClusterMap();
1712   Int_t nPointsForPid=0;
1713   for(Int_t i=2; i<6; i++){
1714     if(clumap&(1<<i)) ++nPointsForPid;
1715   }
1716   Float_t mom=track->P();
1717   
1718   //check for ITS standalone tracks
1719   Bool_t isSA=kTRUE;
1720   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1721   
1722   //TODO: in case of the electron, use the SA parametrisation,
1723   //      this needs to be changed if ITS provides a parametrisation
1724   //      for electrons also for ITS+TPC tracks
1725   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
1726 }
1727
1728 //______________________________________________________________________________
1729 Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
1730 {
1731   //
1732   // Calculate the number of sigmas in the TPC
1733   //
1734   
1735   AliVTrack *track=(AliVTrack*)vtrack;
1736   
1737   Double_t nSigma = -999.;
1738   
1739   //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
1740   if (fTuneMConData) {
1741     Double_t mom  = track->GetTPCmomentum();
1742     Double_t sig  = this->GetTPCsignalTunedOnData(track);
1743     UInt_t   sigN = track->GetTPCsignalN();
1744   
1745     if (sigN>0)
1746       nSigma = fTPCResponse.GetNumberOfSigmas(mom, sig, sigN, type);
1747   }
1748   else {
1749     nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1750   }
1751   
1752   return nSigma;
1753 }
1754
1755 //______________________________________________________________________________
1756 Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
1757 {
1758   //
1759   // Calculate the number of sigmas in the EMCAL
1760   //
1761   
1762   AliVTrack *track=(AliVTrack*)vtrack;
1763   
1764   AliVCluster *matchedClus = NULL;
1765   
1766   Double_t mom     = -1.;
1767   Double_t pt      = -1.;
1768   Double_t EovP    = -1.;
1769   Double_t fClsE   = -1.;
1770   
1771   Int_t nMatchClus = -1;
1772   Int_t charge     = 0;
1773   
1774   // Track matching
1775   nMatchClus = track->GetEMCALcluster();
1776   if(nMatchClus > -1){
1777     
1778     mom    = track->P();
1779     pt     = track->Pt();
1780     charge = track->Charge();
1781     
1782     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1783     
1784     if(matchedClus){
1785       
1786       // matched cluster is EMCAL
1787       if(matchedClus->IsEMCAL()){
1788         
1789         fClsE       = matchedClus->E();
1790         EovP        = fClsE/mom;
1791         
1792         
1793         // NSigma value really meaningful only for electrons!
1794         return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1795       }
1796     }
1797   }
1798   
1799   return -999;
1800   
1801 }
1802
1803
1804 //______________________________________________________________________________
1805 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1806 {
1807   //
1808   // Compute PID response of 'detCode'
1809   //
1810
1811   switch (detCode){
1812     case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
1813     case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
1814     case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
1815     case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
1816     case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
1817     case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
1818     case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
1819     default: return kDetNoSignal;
1820   }
1821 }
1822
1823 //______________________________________________________________________________
1824 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1825 {
1826   //
1827   // Compute PID response for the ITS
1828   //
1829   
1830   // look for cached value first
1831   // only the non SA tracks are cached
1832   if (track->GetDetectorPID()){
1833     return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
1834   }
1835   
1836   // set flat distribution (no decision)
1837   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1838   
1839   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
1840     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
1841   
1842   //check for ITS standalone tracks
1843   Bool_t isSA=kTRUE;
1844   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
1845
1846   Double_t mom=track->P();
1847   Double_t dedx=track->GetITSsignal();
1848   Double_t momITS=mom;
1849   UChar_t clumap=track->GetITSClusterMap();
1850   Int_t nPointsForPid=0;
1851   for(Int_t i=2; i<6; i++){
1852     if(clumap&(1<<i)) ++nPointsForPid;
1853   }
1854
1855   if(nPointsForPid<3) { // track not to be used for combined PID purposes
1856     //       track->ResetStatus(AliVTrack::kITSpid);
1857     return kDetNoSignal;
1858   }
1859
1860   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1861   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
1862     Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
1863     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1864     Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
1865     //TODO: in case of the electron, use the SA parametrisation,
1866     //      this needs to be changed if ITS provides a parametrisation
1867     //      for electrons also for ITS+TPC tracks
1868     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
1869     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1870       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1871     } else {
1872       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1873       mismatch=kFALSE;
1874     }
1875
1876     // Check for particles heavier than (AliPID::kSPECIES - 1)
1877     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
1878
1879   }
1880
1881   if (mismatch){
1882     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
1883     return kDetNoSignal;
1884   }
1885
1886
1887   return kDetPidOk;
1888 }
1889 //______________________________________________________________________________
1890 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1891 {
1892   //
1893   // Compute PID response for the TPC
1894   //
1895   
1896   // set flat distribution (no decision)
1897   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1898   
1899   // check quality of the track
1900   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
1901   
1902   Double_t dedx=track->GetTPCsignal();
1903   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
1904   
1905   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
1906   
1907   Double_t bethe = 0.;
1908   Double_t sigma = 0.;
1909   
1910   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1911     AliPID::EParticleType type=AliPID::EParticleType(j);
1912     
1913     //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
1914     if (fTuneMConData) {
1915       Double_t mom = track->GetTPCmomentum();
1916       bethe=fTPCResponse.GetExpectedSignal(mom,type);
1917       sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
1918     }
1919     else {
1920       bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1921       sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
1922     }
1923     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
1924       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
1925     } else {
1926       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
1927       mismatch=kFALSE;
1928     }
1929   }
1930   
1931   if (mismatch){
1932     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1933     return kDetNoSignal;
1934   }
1935   
1936   return kDetPidOk;
1937 }
1938 //______________________________________________________________________________
1939 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
1940 {
1941   //
1942   // Compute PID response for the
1943   //
1944   
1945   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
1946   
1947   // set flat distribution (no decision)
1948   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1949   
1950   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
1951   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
1952   
1953   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
1954   for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
1955     AliPID::EParticleType type=AliPID::EParticleType(j);
1956     Double_t nsigmas=GetNumberOfSigmasTOF(track,type) + meanCorrFactor;
1957     
1958     Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
1959     Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
1960     if (TMath::Abs(nsigmas) > (fRange+2)) {
1961       if(nsigmas < fTOFtail)
1962         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
1963       else
1964         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
1965     } else{
1966       if(nsigmas < fTOFtail)
1967         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
1968       else
1969         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
1970     }
1971     
1972     if (TMath::Abs(nsigmas)<5.){
1973       Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
1974       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
1975     }
1976   }
1977   
1978   if (mismatch){
1979     return kDetMismatch;
1980   }
1981   
1982   // TODO: Light nuclei
1983   
1984   return kDetPidOk;
1985 }
1986 //______________________________________________________________________________
1987 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const
1988 {
1989   //
1990   // Compute PID response for the
1991   //
1992   
1993   UInt_t TRDslicesForPID[2];
1994   SetTRDSlices(TRDslicesForPID,PIDmethod);
1995   // set flat distribution (no decision)
1996   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1997   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
1998   
1999   Float_t mom[6]={0.};
2000   Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
2001   Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
2002   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  TRDslicesForPID[0], TRDslicesForPID[1], nslices));
2003   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
2004     mom[ilayer] = track->GetTRDmomentum(ilayer);
2005     for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
2006       dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
2007     }
2008   }
2009   fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
2010   return kDetPidOk;
2011   
2012 }
2013 //______________________________________________________________________________
2014 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2015 {
2016   //
2017   // Compute PID response for the EMCAL
2018   //
2019   
2020   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2021   
2022   AliVCluster *matchedClus = NULL;
2023   
2024   Double_t mom     = -1.;
2025   Double_t pt      = -1.;
2026   Double_t EovP    = -1.;
2027   Double_t fClsE   = -1.;
2028   
2029   Int_t nMatchClus = -1;
2030   Int_t charge     = 0;
2031   
2032   // Track matching
2033   nMatchClus = track->GetEMCALcluster();
2034   
2035   if(nMatchClus > -1){
2036     
2037     mom    = track->P();
2038     pt     = track->Pt();
2039     charge = track->Charge();
2040     
2041     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2042     
2043     if(matchedClus){
2044       
2045       // matched cluster is EMCAL
2046       if(matchedClus->IsEMCAL()){
2047         
2048         fClsE       = matchedClus->E();
2049         EovP        = fClsE/mom;
2050         
2051         
2052         // compute the probabilities
2053         if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
2054           
2055           // in case everything is OK
2056           return kDetPidOk;
2057         }
2058       }
2059     }
2060   }
2061   
2062   // in all other cases set flat distribution (no decision)
2063   for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
2064   return kDetNoSignal;
2065   
2066 }
2067 //______________________________________________________________________________
2068 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
2069 {
2070   //
2071   // Compute PID response for the PHOS
2072   //
2073   
2074   // set flat distribution (no decision)
2075   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2076   return kDetNoSignal;
2077 }
2078 //______________________________________________________________________________
2079 AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2080 {
2081   //
2082   // Compute PID response for the HMPID
2083   //
2084   // set flat distribution (no decision)
2085   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2086   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
2087   
2088   track->GetHMPIDpid(p);
2089   
2090   return kDetPidOk;
2091 }