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