#90105
[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 }