#92331: preparing to move AliTOFHeader to AOD
[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 <TSpline.h>
31 #include <TFile.h>
32
33 #include <AliVEvent.h>
34 #include <AliVTrack.h>
35 #include <AliLog.h>
36 #include <AliPID.h>
37 #include <AliOADBContainer.h>
38 #include <AliTRDPIDParams.h>
39 #include <AliTRDPIDReference.h>
40
41 #include "AliPIDResponse.h"
42
43 #include "AliCentrality.h"
44
45 ClassImp(AliPIDResponse);
46
47 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
48 TNamed("PIDResponse","PIDResponse"),
49 fITSResponse(isMC),
50 fTPCResponse(),
51 fTRDResponse(),
52 fTOFResponse(),
53 fEMCALResponse(),
54 fRange(5.),
55 fITSPIDmethod(kITSTruncMean),
56 fIsMC(isMC),
57 fOADBPath(),
58 fBeamType("PP"),
59 fLHCperiod(),
60 fMCperiodTPC(),
61 fMCperiodUser(),
62 fCurrentFile(),
63 fRecoPass(0),
64 fRecoPassUser(-1),
65 fRun(0),
66 fOldRun(0),
67 fArrPidResponseMaster(0x0),
68 fResolutionCorrection(0x0),
69 fTRDPIDParams(0x0),
70 fTRDPIDReference(0x0),
71 fTOFTimeZeroType(kBest_T0),
72 fTOFres(100.),
73 fTOFtail(1.1),
74 fEMCALPIDParams(0x0),
75 fCurrentEvent(0x0),
76 fCurrCentrality(0.0)
77 {
78   //
79   // default ctor
80   //
81   AliLog::SetClassDebugLevel("AliPIDResponse",10);
82   AliLog::SetClassDebugLevel("AliESDpid",10);
83   AliLog::SetClassDebugLevel("AliAODpidUtil",10);
84
85   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
86 }
87
88 //______________________________________________________________________________
89 AliPIDResponse::~AliPIDResponse()
90 {
91   //
92   // dtor
93   //
94   delete fArrPidResponseMaster;
95   delete fTRDPIDParams;
96   delete fTRDPIDReference;
97 }
98
99 //______________________________________________________________________________
100 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
101 TNamed(other),
102 fITSResponse(other.fITSResponse),
103 fTPCResponse(other.fTPCResponse),
104 fTRDResponse(other.fTRDResponse),
105 fTOFResponse(other.fTOFResponse),
106 fEMCALResponse(other.fEMCALResponse),
107 fRange(other.fRange),
108 fITSPIDmethod(other.fITSPIDmethod),
109 fIsMC(other.fIsMC),
110 fOADBPath(other.fOADBPath),
111 fBeamType("PP"),
112 fLHCperiod(),
113 fMCperiodTPC(),
114 fMCperiodUser(other.fMCperiodUser),
115 fCurrentFile(),
116 fRecoPass(0),
117 fRecoPassUser(other.fRecoPassUser),
118 fRun(0),
119 fOldRun(0),
120 fArrPidResponseMaster(0x0),
121 fResolutionCorrection(0x0),
122 fTRDPIDParams(0x0),
123 fTRDPIDReference(0x0),
124 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
125 fTOFres(100.),
126 fTOFtail(1.1),
127 fEMCALPIDParams(0x0),
128 fCurrentEvent(0x0),
129 fCurrCentrality(0.0)
130 {
131   //
132   // copy ctor
133   //
134   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
135 }
136
137 //______________________________________________________________________________
138 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
139 {
140   //
141   // copy ctor
142   //
143   if(this!=&other) {
144     delete fArrPidResponseMaster;
145     TNamed::operator=(other);
146     fITSResponse=other.fITSResponse;
147     fTPCResponse=other.fTPCResponse;
148     fTRDResponse=other.fTRDResponse;
149     fTOFResponse=other.fTOFResponse;
150     fEMCALResponse=other.fEMCALResponse;
151     fRange=other.fRange;
152     fITSPIDmethod=other.fITSPIDmethod;
153     fOADBPath=other.fOADBPath;
154     fIsMC=other.fIsMC;
155     fBeamType="PP";
156     fLHCperiod="";
157     fMCperiodTPC="";
158     fMCperiodUser=other.fMCperiodUser;
159     fCurrentFile="";
160     fRecoPass=0;
161     fRecoPassUser=other.fRecoPassUser;
162     fRun=0;
163     fOldRun=0;
164     fArrPidResponseMaster=0x0;
165     fResolutionCorrection=0x0;
166     fTRDPIDParams=0x0;
167     fTRDPIDReference=0x0;
168     fEMCALPIDParams=0x0;
169     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
170     fTOFTimeZeroType=AliPIDResponse::kBest_T0;
171     fTOFres=100.;
172     fTOFtail=1.1;
173     fCurrentEvent=other.fCurrentEvent;
174   }
175   return *this;
176 }
177
178 //______________________________________________________________________________
179 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
180 {
181   //
182   // NumberOfSigmas for 'detCode'
183   //
184
185   switch (detCode){
186     case kDetITS: return NumberOfSigmasITS(track, type); break;
187     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
188     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
189 //     case kDetTRD: return ComputeTRDProbability(track, type); break;
190 //     case kDetPHOS: return ComputePHOSProbability(track, type); break;
191 //     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
192 //     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
193     default: return -999.;
194   }
195
196 }
197
198 //______________________________________________________________________________
199 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
200
201   AliVCluster *matchedClus = NULL;
202
203   Double_t mom     = -1.; 
204   Double_t pt      = -1.; 
205   Double_t EovP    = -1.;
206   Double_t fClsE   = -1.;
207   
208   Int_t nMatchClus = -1;
209   Int_t charge     = 0;
210   
211   // Track matching
212   nMatchClus = track->GetEMCALcluster();
213   if(nMatchClus > -1){
214
215     mom    = track->P();
216     pt     = track->Pt();
217     charge = track->Charge();
218     
219     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
220     
221     if(matchedClus){
222       
223     // matched cluster is EMCAL
224     if(matchedClus->IsEMCAL()){
225       
226       fClsE       = matchedClus->E();
227       EovP        = fClsE/mom;
228       
229       
230       // NSigma value really meaningful only for electrons!
231       return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); 
232     }
233   }
234   }
235
236   return -999;
237   
238 }
239
240 //______________________________________________________________________________
241 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
242 {
243   //
244   // Compute PID response of 'detCode'
245   //
246
247   switch (detCode){
248     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
249     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
250     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
251     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
252     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
253     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
254     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
255     default: return kDetNoSignal;
256   }
257 }
258
259 //______________________________________________________________________________
260 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
261 {
262   //
263   // Compute PID response for the ITS
264   //
265
266   // set flat distribution (no decision)
267   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
268
269   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
270     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
271   
272   Double_t mom=track->P();
273   Double_t dedx=track->GetITSsignal();
274   Bool_t isSA=kTRUE;
275   Double_t momITS=mom;
276   ULong_t trStatus=track->GetStatus();
277   if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
278   UChar_t clumap=track->GetITSClusterMap();
279   Int_t nPointsForPid=0;
280   for(Int_t i=2; i<6; i++){
281     if(clumap&(1<<i)) ++nPointsForPid;
282   }
283   
284   if(nPointsForPid<3) { // track not to be used for combined PID purposes
285     //       track->ResetStatus(AliVTrack::kITSpid);
286     return kDetNoSignal;
287   }
288
289   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
290   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
291     Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
292     Double_t bethe=fITSResponse.Bethe(momITS,mass);
293     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
294     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
295       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
296     } else {
297       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
298       mismatch=kFALSE;
299     }
300
301     // Check for particles heavier than (AliPID::kSPECIES - 1)
302     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
303
304   }
305
306   if (mismatch){
307     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
308     return kDetNoSignal;
309   }
310
311     
312   return kDetPidOk;
313 }
314 //______________________________________________________________________________
315 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
316 {
317   //
318   // Compute PID response for the TPC
319   //
320
321   // set flat distribution (no decision)
322   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
323
324   // check quality of the track
325   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
326
327   Double_t mom = track->GetTPCmomentum();
328
329   Double_t dedx=track->GetTPCsignal();
330   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
331
332   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
333     AliPID::EParticleType type=AliPID::EParticleType(j);
334     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
335     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
336     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
337       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
338     } else {
339       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
340       mismatch=kFALSE;
341     }
342
343     // TODO: Light nuclei, also in TPC pid response
344     
345     // Check for particles heavier than (AliPID::kSPECIES - 1)
346 //     if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
347
348   }
349
350   if (mismatch){
351     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
352     return kDetNoSignal;
353   }
354
355   return kDetPidOk;
356 }
357 //______________________________________________________________________________
358 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
359 {
360   //
361   // Compute PID response for the
362   //
363
364   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
365
366   // set flat distribution (no decision)
367   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
368   
369   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
370   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
371   
372   Double_t time[AliPID::kSPECIESN];
373   track->GetIntegratedTimes(time);
374   
375   Double_t sigma[AliPID::kSPECIES];
376   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
377     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
378   }
379   
380   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
381   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
382     AliPID::EParticleType type=AliPID::EParticleType(j);
383     Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
384
385     Double_t sig = sigma[j];
386     if (TMath::Abs(nsigmas) > (fRange+2)) {
387       if(nsigmas < fTOFtail)
388         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
389       else
390         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
391     } else{
392       if(nsigmas < fTOFtail)
393         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
394       else
395         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
396     }
397
398     /* OLD Gaussian shape
399     if (TMath::Abs(nsigmas) > (fRange+2)) {
400       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
401     } else
402       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
403     */
404
405     if (TMath::Abs(nsigmas)<5.){
406       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
407       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
408     }
409   }
410
411   if (mismatch){
412     return kDetMismatch;    
413   }
414
415     // TODO: Light nuclei
416     
417   return kDetPidOk;
418 }
419 //______________________________________________________________________________
420 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
421 {
422   //
423   // Compute PID response for the
424   //
425
426   // set flat distribution (no decision)
427   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
428   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
429
430   Float_t mom[6];
431   Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
432   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
433   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
434   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
435     mom[ilayer] = track->GetTRDmomentum(ilayer);
436     for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
437       dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
438     }
439   }
440   fTRDResponse.GetResponse(nslices, dedx, mom, p);
441   return kDetPidOk;
442 }
443 //______________________________________________________________________________
444 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
445 {
446   //
447   // Compute PID response for the EMCAL
448   //
449
450   AliVCluster *matchedClus = NULL;
451
452   Double_t mom     = -1.; 
453   Double_t pt      = -1.; 
454   Double_t EovP    = -1.;
455   Double_t fClsE   = -1.;
456   
457   Int_t nMatchClus = -1;
458   Int_t charge     = 0;
459   
460   // Track matching
461   nMatchClus = track->GetEMCALcluster();
462
463   if(nMatchClus > -1){
464
465     mom    = track->P();
466     pt     = track->Pt();
467     charge = track->Charge();
468     
469     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
470     
471     if(matchedClus){    
472
473     // matched cluster is EMCAL
474     if(matchedClus->IsEMCAL()){
475
476       fClsE       = matchedClus->E();
477       EovP        = fClsE/mom;
478       
479       
480       // compute the probabilities 
481       if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){     
482
483         // in case everything is OK
484         return kDetPidOk;
485         
486       }
487     }
488   }
489   }
490   
491   // in all other cases set flat distribution (no decision)
492   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
493   return kDetNoSignal;
494   
495 }
496 //______________________________________________________________________________
497 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
498 {
499   //
500   // Compute PID response for the PHOS
501   //
502
503   // set flat distribution (no decision)
504   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
505   return kDetNoSignal;
506 }
507 //______________________________________________________________________________
508 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
509 {
510   //
511   // Compute PID response for the HMPID
512   //
513
514   // set flat distribution (no decision)
515   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
516   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
517
518   track->GetHMPIDpid(p);
519   
520   return kDetPidOk;
521 }
522
523 //______________________________________________________________________________
524 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
525 {
526   //
527   // Apply settings for the current event
528   //
529   fRecoPass=pass;
530   
531   fCurrentEvent=0x0;
532   if (!event) return;
533   fCurrentEvent=event;
534   fRun=event->GetRunNumber();
535   
536   if (fRun!=fOldRun){
537     ExecNewRun();
538     fOldRun=fRun;
539   }
540   
541   //TPC resolution parametrisation PbPb
542   if ( fResolutionCorrection ){
543     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
544     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
545   }
546   
547   //TOF resolution
548   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
549
550   // Get and set centrality
551   AliCentrality *centrality = event->GetCentrality();
552   if(centrality){
553     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
554   }
555   else{
556     fCurrCentrality = -1;
557   }
558 }
559
560 //______________________________________________________________________________
561 void AliPIDResponse::ExecNewRun()
562 {
563   //
564   // Things to Execute upon a new run
565   //
566   SetRecoInfo();
567   
568   SetITSParametrisation();
569   
570   SetTPCPidResponseMaster();
571   SetTPCParametrisation();
572
573   SetTRDPidResponseMaster(); 
574   InitializeTRDResponse();
575
576   SetEMCALPidResponseMaster(); 
577   InitializeEMCALResponse();
578   
579   fTOFResponse.SetTimeResolution(fTOFres);
580 }
581
582 //_____________________________________________________
583 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
584 {
585   //
586   // Get TPC multiplicity in bins of 150
587   //
588   
589   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
590   Double_t tpcMulti=0.;
591   if(vertexTPC){
592     Double_t vertexContribTPC=vertexTPC->GetNContributors();
593     tpcMulti=vertexContribTPC/150.;
594     if (tpcMulti>20.) tpcMulti=20.;
595   }
596   
597   return tpcMulti;
598 }
599
600 //______________________________________________________________________________
601 void AliPIDResponse::SetRecoInfo()
602 {
603   //
604   // Set reconstruction information
605   //
606   
607   //reset information
608   fLHCperiod="";
609   fMCperiodTPC="";
610   
611   fBeamType="";
612     
613   fBeamType="PP";
614   
615   TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
616   //find the period by run number (UGLY, but not stored in ESD and AOD... )
617   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
618   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
619   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
620   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
621   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
622   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
623   else if (fRun>=136851&&fRun<=139517) {
624     fLHCperiod="LHC10H";
625     fMCperiodTPC="LHC10H8";
626     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
627     fBeamType="PBPB";
628   }
629   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
630   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
631   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
632   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
633   else if (fRun>=166529) {
634     fLHCperiod="LHC11H";
635     fMCperiodTPC="LHC11A10";
636     fBeamType="PBPB";
637   }
638
639
640   //exception new pp MC productions from 2011
641   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
642 }
643
644 //______________________________________________________________________________
645 void AliPIDResponse::SetITSParametrisation()
646 {
647   //
648   // Set the ITS parametrisation
649   //
650 }
651
652 //______________________________________________________________________________
653 void AliPIDResponse::SetTPCPidResponseMaster()
654 {
655   //
656   // Load the TPC pid response functions from the OADB
657   //
658   //don't load twice for the moment
659    if (fArrPidResponseMaster) return;
660  
661
662   //reset the PID response functions
663   delete fArrPidResponseMaster;
664   fArrPidResponseMaster=0x0;
665   
666   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
667   
668   TFile *f=TFile::Open(fileName.Data());
669   if (f && f->IsOpen() && !f->IsZombie()){
670     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
671   }
672   delete f;
673   
674   if (!fArrPidResponseMaster){
675     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
676     return;
677   }
678   fArrPidResponseMaster->SetOwner();
679 }
680
681 //______________________________________________________________________________
682 void AliPIDResponse::SetTPCParametrisation()
683 {
684   //
685   // Change BB parametrisation for current run
686   //
687   
688   if (fLHCperiod.IsNull()) {
689     AliFatal("No period set, not changing parametrisation");
690     return;
691   }
692   
693   //
694   // Set default parametrisations for data and MC
695   //
696   
697   //data type
698   TString datatype="DATA";
699   //in case of mc fRecoPass is per default 1
700   if (fIsMC) {
701     datatype="MC";
702     fRecoPass=1;
703   }
704   
705   //
706   //reset old splines
707   //
708   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
709     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
710   }
711   
712   //
713   //set the new PID splines
714   //
715   TString period=fLHCperiod;
716   if (fArrPidResponseMaster){
717     TObject *grAll=0x0;
718     //for MC don't use period information
719 //     if (fIsMC) period="[A-Z0-9]*";
720     //for MC use MC period information
721     if (fIsMC) period=fMCperiodTPC;
722 //pattern for the default entry (valid for all particles)
723     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
724     
725     //loop over entries and filter them
726     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
727       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
728       if (responseFunction==0x0) continue;
729       TString responseName=responseFunction->GetName();
730       
731       if (!reg.MatchB(responseName)) continue;
732       
733       TObjArray *arr=reg.MatchS(responseName);
734       TString particleName=arr->At(1)->GetName();
735       delete arr;
736       if (particleName.IsNull()) continue;
737       if (particleName=="ALL") grAll=responseFunction;
738       else {
739         //find particle id
740         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
741           TString particle=AliPID::ParticleName(ispec);
742           particle.ToUpper();
743           if ( particle == particleName ){
744             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
745             fTPCResponse.SetUseDatabase(kTRUE);
746             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
747             break;
748           }
749         }
750       }
751     }
752     
753     //set default response function to all particles which don't have a specific one
754     if (grAll){
755       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
756         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
757           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
758           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
759         }
760       }
761     }
762   }
763   
764   //
765   // Setup resolution parametrisation
766   //
767   
768   //default
769   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
770   
771   if (fRun>=122195){
772     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
773   }
774   if (fArrPidResponseMaster)
775   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
776   
777   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
778 }
779
780 //______________________________________________________________________________
781 void AliPIDResponse::SetTRDPidResponseMaster()
782 {
783   //
784   // Load the TRD pid params and references from the OADB
785   //
786   if(fTRDPIDParams) return;
787   AliOADBContainer contParams("contParams"); 
788
789   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
790   if(statusPars){
791     AliError("Failed initializing PID Params from OADB");
792   } else {
793     AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
794     fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
795     if(!fTRDPIDParams){
796       AliError(Form("TRD Params not found in run %d", fRun));
797     }
798   }
799
800   AliOADBContainer contRefs("contRefs");
801   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
802   if(statusRefs){
803     AliInfo("Failed Loading References for TRD");
804   } else {
805     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
806     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
807     if(!fTRDPIDReference){
808       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
809     }
810   }
811 }
812
813 //______________________________________________________________________________
814 void AliPIDResponse::InitializeTRDResponse(){
815   //
816   // Set PID Params and references to the TRD PID response
817   // 
818   fTRDResponse.SetPIDParams(fTRDPIDParams);
819   fTRDResponse.Load(fTRDPIDReference);
820   if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
821     fTRDslicesForPID[0] = 0;
822     fTRDslicesForPID[1] = 7;
823   }
824 }
825
826 //_________________________________________________________________________
827 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
828   //
829   // Check whether track is identified as electron under a given electron efficiency hypothesis
830   //
831   Double_t probs[AliPID::kSPECIES];
832   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
833
834   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
835   // Take mean of the TRD momenta in the given tracklets
836   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
837   Int_t nmomenta = 0;
838   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
839     if(vtrack->GetTRDmomentum(iPl) > 0.){
840       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); 
841     }
842   }
843   p = TMath::Mean(nmomenta, trdmomenta);
844
845   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
846 }
847
848 //______________________________________________________________________________
849 void AliPIDResponse::SetEMCALPidResponseMaster()
850 {
851   //
852   // Load the EMCAL pid response functions from the OADB
853   //
854   TObjArray* fEMCALPIDParamsRun      = NULL;
855   TObjArray* fEMCALPIDParamsPass     = NULL;
856
857   if(fEMCALPIDParams) return;
858   AliOADBContainer contParams("contParams"); 
859
860   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
861   if(statusPars){
862     AliError("Failed initializing PID Params from OADB");
863   } 
864   else {
865     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
866
867     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
868     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
869     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
870
871     if(!fEMCALPIDParams){
872       AliError(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
873       AliInfo("Will take the standard LHC11d instead ...");
874
875       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
876       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
877       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
878
879       if(!fEMCALPIDParams){
880         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));     
881       }
882     }
883   }
884 }
885
886 //______________________________________________________________________________
887 void AliPIDResponse::InitializeEMCALResponse(){
888   //
889   // Set PID Params to the EMCAL PID response
890   // 
891   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
892
893 }
894 //_________________________________________________________________________
895 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
896   //
897   // Set TOF response function
898   // Input option for event_time used
899   //
900   
901     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
902     if(t0spread < 10) t0spread = 80;
903
904     // T0 from TOF algorithm
905
906     Bool_t flagT0TOF=kFALSE;
907     Bool_t flagT0T0=kFALSE;
908     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
909     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
910     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
911
912     // T0-TOF arrays
913     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
914     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
915     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
916       estimatedT0event[i]=0.0;
917       estimatedT0resolution[i]=0.0;
918       startTimeMask[i] = 0;
919     }
920
921     Float_t resT0A=75,resT0C=65,resT0AC=55;
922     if(vevent->GetT0TOF()){ // check if T0 detector information is available
923         flagT0T0=kTRUE;
924     }
925
926
927     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
928
929     if (tofHeader) { // read global info and T0-TOF
930       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
931       t0spread = tofHeader->GetT0spread(); // read t0 sprad
932       if(t0spread < 10) t0spread = 80;
933
934       flagT0TOF=kTRUE;
935       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
936         startTime[i]=tofHeader->GetDefaultEventTimeVal();
937         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
938         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
939       }
940
941       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
942       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
943       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
944       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
945         Int_t icurrent = (Int_t)ibin->GetAt(j);
946         startTime[icurrent]=t0Bin->GetAt(j);
947         startTimeRes[icurrent]=t0ResBin->GetAt(j);
948         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
949       }
950     }
951
952     // for cut of 3 sigma on t0 spread
953     Float_t t0cut = 3 * t0spread;
954     if(t0cut < 500) t0cut = 500;
955
956     if(option == kFILL_T0){ // T0-FILL is used
957         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
958           estimatedT0event[i]=0.0;
959           estimatedT0resolution[i]=t0spread;
960         }
961         fTOFResponse.SetT0event(estimatedT0event);
962         fTOFResponse.SetT0resolution(estimatedT0resolution);
963     }
964
965     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
966         if(flagT0TOF){
967             fTOFResponse.SetT0event(startTime);
968             fTOFResponse.SetT0resolution(startTimeRes);
969             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
970               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
971               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
972             }
973         }
974         else{
975             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
976               estimatedT0event[i]=0.0;
977               estimatedT0resolution[i]=t0spread;
978               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
979             }
980             fTOFResponse.SetT0event(estimatedT0event);
981             fTOFResponse.SetT0resolution(estimatedT0resolution);
982         }
983     }
984     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
985         Float_t t0AC=-10000;
986         Float_t t0A=-10000;
987         Float_t t0C=-10000;
988         if(flagT0T0){
989             t0AC= vevent->GetT0TOF()[0];
990             t0A= vevent->GetT0TOF()[1];
991             t0C= vevent->GetT0TOF()[2];
992         }
993
994         Float_t t0t0Best = 0;
995         Float_t t0t0BestRes = 9999;
996         Int_t t0used=0;
997         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
998             t0t0Best = t0AC;
999             t0t0BestRes = resT0AC;
1000             t0used=6;
1001         }
1002         else if(TMath::Abs(t0C) < t0cut){
1003             t0t0Best = t0C;
1004             t0t0BestRes = resT0C;
1005             t0used=4;
1006         }
1007         else if(TMath::Abs(t0A) < t0cut){
1008             t0t0Best = t0A;
1009             t0t0BestRes = resT0A;
1010             t0used=2;
1011         }
1012
1013         if(flagT0TOF){ // if T0-TOF info is available
1014             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1015                 if(t0t0BestRes < 999){
1016                   if(startTimeRes[i] < t0spread){
1017                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1018                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1019                     estimatedT0event[i]=t0best / wtot;
1020                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1021                     startTimeMask[i] = t0used+1;
1022                   }
1023                   else {
1024                     estimatedT0event[i]=t0t0Best;
1025                     estimatedT0resolution[i]=t0t0BestRes;
1026                     startTimeMask[i] = t0used;
1027                   }
1028                 }
1029                 else{
1030                   estimatedT0event[i]=startTime[i];
1031                   estimatedT0resolution[i]=startTimeRes[i];
1032                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1033                 }
1034                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1035             }
1036             fTOFResponse.SetT0event(estimatedT0event);
1037             fTOFResponse.SetT0resolution(estimatedT0resolution);
1038         }
1039         else{ // if no T0-TOF info is available
1040             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1041               fTOFResponse.SetT0binMask(i,t0used);
1042               if(t0t0BestRes < 999){
1043                 estimatedT0event[i]=t0t0Best;
1044                 estimatedT0resolution[i]=t0t0BestRes;
1045               }
1046               else{
1047                 estimatedT0event[i]=0.0;
1048                 estimatedT0resolution[i]=t0spread;
1049               }
1050             }
1051             fTOFResponse.SetT0event(estimatedT0event);
1052             fTOFResponse.SetT0resolution(estimatedT0resolution);
1053         }
1054     }
1055
1056     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1057         Float_t t0AC=-10000;
1058         Float_t t0A=-10000;
1059         Float_t t0C=-10000;
1060         if(flagT0T0){
1061             t0AC= vevent->GetT0TOF()[0];
1062             t0A= vevent->GetT0TOF()[1];
1063             t0C= vevent->GetT0TOF()[2];
1064         }
1065
1066         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1067             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1068               estimatedT0event[i]=t0AC;
1069               estimatedT0resolution[i]=resT0AC;
1070               fTOFResponse.SetT0binMask(i,6);
1071             }
1072         }
1073         else if(TMath::Abs(t0C) < t0cut){
1074             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1075               estimatedT0event[i]=t0C;
1076               estimatedT0resolution[i]=resT0C;
1077               fTOFResponse.SetT0binMask(i,4);
1078             }
1079         }
1080         else if(TMath::Abs(t0A) < t0cut){
1081             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1082               estimatedT0event[i]=t0A;
1083               estimatedT0resolution[i]=resT0A;
1084               fTOFResponse.SetT0binMask(i,2);
1085             }
1086         }
1087         else{
1088             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1089               estimatedT0event[i]=0.0;
1090               estimatedT0resolution[i]=t0spread;
1091               fTOFResponse.SetT0binMask(i,0);
1092             }
1093         }
1094         fTOFResponse.SetT0event(estimatedT0event);
1095         fTOFResponse.SetT0resolution(estimatedT0resolution);
1096     }
1097     delete [] startTime;
1098     delete [] startTimeRes;
1099     delete [] startTimeMask;
1100     delete [] estimatedT0event;
1101     delete [] estimatedT0resolution;
1102 }