Moving the classes that belong to the following libraries: STEERBase, ESD, CDB, 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 <AliTRDPIDReference.h>
39
40 #include "AliPIDResponse.h"
41
42 ClassImp(AliPIDResponse);
43
44 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
45 TNamed("PIDResponse","PIDResponse"),
46 fITSResponse(isMC),
47 fTPCResponse(),
48 fTRDResponse(),
49 fTOFResponse(),
50 fRange(5.),
51 fITSPIDmethod(kITSTruncMean),
52 fIsMC(isMC),
53 fOADBPath(),
54 fBeamType("PP"),
55 fLHCperiod(),
56 fMCperiodTPC(),
57 fMCperiodUser(),
58 fCurrentFile(),
59 fRecoPass(0),
60 fRecoPassUser(-1),
61 fRun(0),
62 fOldRun(0),
63 fArrPidResponseMaster(0x0),
64 fResolutionCorrection(0x0),
65 fTRDPIDParams(0x0),
66 fTRDPIDReference(0x0),
67 fTOFTimeZeroType(kBest_T0),
68 fTOFres(100.)
69 {
70   //
71   // default ctor
72   //
73   AliLog::SetClassDebugLevel("AliPIDResponse",10);
74   AliLog::SetClassDebugLevel("AliESDpid",10);
75   AliLog::SetClassDebugLevel("AliAODpidUtil",10);
76
77   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
78 }
79
80 //______________________________________________________________________________
81 AliPIDResponse::~AliPIDResponse()
82 {
83   //
84   // dtor
85   //
86   delete fArrPidResponseMaster;
87   delete fTRDPIDParams;
88   delete fTRDPIDReference;
89 }
90
91 //______________________________________________________________________________
92 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
93 TNamed(other),
94 fITSResponse(other.fITSResponse),
95 fTPCResponse(other.fTPCResponse),
96 fTRDResponse(other.fTRDResponse),
97 fTOFResponse(other.fTOFResponse),
98 fRange(other.fRange),
99 fITSPIDmethod(other.fITSPIDmethod),
100 fIsMC(other.fIsMC),
101 fOADBPath(other.fOADBPath),
102 fBeamType("PP"),
103 fLHCperiod(),
104 fMCperiodTPC(),
105 fMCperiodUser(other.fMCperiodUser),
106 fCurrentFile(),
107 fRecoPass(0),
108 fRecoPassUser(other.fRecoPassUser),
109 fRun(0),
110 fOldRun(0),
111 fArrPidResponseMaster(0x0),
112 fResolutionCorrection(0x0),
113 fTRDPIDParams(0x0),
114 fTRDPIDReference(0x0),
115 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
116 fTOFres(100.)
117 {
118   //
119   // copy ctor
120   //
121   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
122 }
123
124 //______________________________________________________________________________
125 AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
126 {
127   //
128   // copy ctor
129   //
130   if(this!=&other) {
131     delete fArrPidResponseMaster;
132     TNamed::operator=(other);
133     fITSResponse=other.fITSResponse;
134     fTPCResponse=other.fTPCResponse;
135     fTRDResponse=other.fTRDResponse;
136     fTOFResponse=other.fTOFResponse;
137     fRange=other.fRange;
138     fITSPIDmethod=other.fITSPIDmethod;
139     fOADBPath=other.fOADBPath;
140     fIsMC=other.fIsMC;
141     fBeamType="PP";
142     fLHCperiod="";
143     fMCperiodTPC="";
144     fMCperiodUser=other.fMCperiodUser;
145     fCurrentFile="";
146     fRecoPass=0;
147     fRecoPassUser=other.fRecoPassUser;
148     fRun=0;
149     fOldRun=0;
150     fArrPidResponseMaster=0x0;
151     fResolutionCorrection=0x0;
152     fTRDPIDParams=0x0;
153     fTRDPIDReference=0x0;
154     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
155     fTOFTimeZeroType=AliPIDResponse::kBest_T0;
156     fTOFres=100.;
157   }
158   return *this;
159 }
160
161 //______________________________________________________________________________
162 Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
163 {
164   //
165   // NumberOfSigmas for 'detCode'
166   //
167
168   switch (detCode){
169     case kDetITS: return NumberOfSigmasITS(track, type); break;
170     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
171     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
172 //     case kDetTRD: return ComputeTRDProbability(track, type); break;
173 //     case kDetPHOS: return ComputePHOSProbability(track, type); break;
174 //     case kDetEMCAL: return ComputeEMCALProbability(track, type); break;
175 //     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
176     default: return -999.;
177   }
178
179 }
180
181 //______________________________________________________________________________
182 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
183 {
184   //
185   // Compute PID response of 'detCode'
186   //
187
188   switch (detCode){
189     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
190     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
191     case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
192     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
193     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
194     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
195     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
196     default: return kDetNoSignal;
197   }
198 }
199
200 //______________________________________________________________________________
201 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
202 {
203   //
204   // Compute PID response for the ITS
205   //
206
207   // set flat distribution (no decision)
208   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
209
210   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
211     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
212   
213   Double_t mom=track->P();
214   Double_t dedx=track->GetITSsignal();
215   Bool_t isSA=kTRUE;
216   Double_t momITS=mom;
217   ULong_t trStatus=track->GetStatus();
218   if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
219   UChar_t clumap=track->GetITSClusterMap();
220   Int_t nPointsForPid=0;
221   for(Int_t i=2; i<6; i++){
222     if(clumap&(1<<i)) ++nPointsForPid;
223   }
224   
225   if(nPointsForPid<3) { // track not to be used for combined PID purposes
226     //       track->ResetStatus(AliVTrack::kITSpid);
227     return kDetNoSignal;
228   }
229
230   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
231   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
232     Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
233     Double_t bethe=fITSResponse.Bethe(momITS,mass);
234     Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
235     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
236       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
237     } else {
238       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
239       mismatch=kFALSE;
240     }
241
242     // Check for particles heavier than (AliPID::kSPECIES - 1)
243     //       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
244
245   }
246
247   if (mismatch){
248     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
249     return kDetNoSignal;
250   }
251
252     
253   return kDetPidOk;
254 }
255 //______________________________________________________________________________
256 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
257 {
258   //
259   // Compute PID response for the TPC
260   //
261
262   // set flat distribution (no decision)
263   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
264
265   // check quality of the track
266   if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
267
268   Double_t mom = track->GetTPCmomentum();
269
270   Double_t dedx=track->GetTPCsignal();
271   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
272
273   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
274     AliPID::EParticleType type=AliPID::EParticleType(j);
275     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
276     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
277     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
278       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
279     } else {
280       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
281       mismatch=kFALSE;
282     }
283
284     // TODO: Light nuclei, also in TPC pid response
285     
286     // Check for particles heavier than (AliPID::kSPECIES - 1)
287 //     if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
288
289   }
290
291   if (mismatch){
292     for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
293     return kDetNoSignal;
294   }
295
296   return kDetPidOk;
297 }
298 //______________________________________________________________________________
299 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
300 {
301   //
302   // Compute PID response for the
303   //
304
305   // set flat distribution (no decision)
306   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
307   
308   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
309   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
310   
311   Double_t time[AliPID::kSPECIESN];
312   track->GetIntegratedTimes(time);
313   
314   Double_t sigma[AliPID::kSPECIES];
315   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
316     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
317   }
318   
319   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
320   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
321     AliPID::EParticleType type=AliPID::EParticleType(j);
322     Double_t nsigmas=NumberOfSigmasTOF(track,type);
323     
324     Double_t sig = sigma[j];
325     if (TMath::Abs(nsigmas) > (fRange+2)) {
326       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
327     } else
328       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
329
330     if (TMath::Abs(nsigmas)<5.){
331       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
332       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
333     }
334   }
335
336   if (mismatch){
337     return kDetMismatch;    
338   }
339
340     // TODO: Light nuclei
341     
342   return kDetPidOk;
343 }
344 //______________________________________________________________________________
345 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
346 {
347   //
348   // Compute PID response for the
349   //
350
351   // set flat distribution (no decision)
352   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
353   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
354
355   Float_t mom[6];
356   Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
357   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
358   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
359   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
360     mom[ilayer] = track->GetTRDmomentum(ilayer);
361     for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
362       dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
363     }
364   }
365   fTRDResponse.GetResponse(nslices, dedx, mom, p);
366   return kDetPidOk;
367 }
368 //______________________________________________________________________________
369 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability(const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
370 {
371   //
372   // Compute PID response for the EMCAL
373   //
374
375   // set flat distribution (no decision)
376   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
377   return kDetNoSignal;
378 }
379 //______________________________________________________________________________
380 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
381 {
382   //
383   // Compute PID response for the PHOS
384   //
385
386   // set flat distribution (no decision)
387   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
388   return kDetNoSignal;
389 }
390 //______________________________________________________________________________
391 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
392 {
393   //
394   // Compute PID response for the HMPID
395   //
396
397   // set flat distribution (no decision)
398   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
399   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
400
401   track->GetHMPIDpid(p);
402   
403   return kDetPidOk;
404 }
405
406 //______________________________________________________________________________
407 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
408 {
409   //
410   // Apply settings for the current event
411   //
412   fRecoPass=pass;
413
414   if (!event) return;
415   fRun=event->GetRunNumber();
416   
417   if (fRun!=fOldRun){
418     ExecNewRun();
419     fOldRun=fRun;
420   }
421   
422   //TPC resolution parametrisation PbPb
423   if ( fResolutionCorrection ){
424     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
425     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
426   }
427   
428   //TOF resolution
429   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
430   
431 }
432
433 //______________________________________________________________________________
434 void AliPIDResponse::ExecNewRun()
435 {
436   //
437   // Things to Execute upon a new run
438   //
439   SetRecoInfo();
440   
441   SetITSParametrisation();
442   
443   SetTPCPidResponseMaster();
444   SetTPCParametrisation();
445   
446   fTOFResponse.SetTimeResolution(fTOFres);
447 }
448
449 //_____________________________________________________
450 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
451 {
452   //
453   // Get TPC multiplicity in bins of 150
454   //
455   
456   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
457   Double_t tpcMulti=0.;
458   if(vertexTPC){
459     Double_t vertexContribTPC=vertexTPC->GetNContributors();
460     tpcMulti=vertexContribTPC/150.;
461     if (tpcMulti>20.) tpcMulti=20.;
462   }
463   
464   return tpcMulti;
465 }
466
467 //______________________________________________________________________________
468 void AliPIDResponse::SetRecoInfo()
469 {
470   //
471   // Set reconstruction information
472   //
473   
474   //reset information
475   fLHCperiod="";
476   fMCperiodTPC="";
477   
478   fBeamType="";
479     
480   fBeamType="PP";
481   
482   TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z]*)/.*");
483   //find the period by run number (UGLY, but not stored in ESD and AOD... )
484   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
485   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
486   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
487   else if (fRun>=127719&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
488   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
489   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
490   else if (fRun>=136851&&fRun<=139517) {
491     fLHCperiod="LHC10H";
492     fMCperiodTPC="LHC10H8";
493     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
494     fBeamType="PBPB";
495   }
496   else if (fRun>=139699) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
497
498   //exception new pp MC productions from 2011
499   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
500 }
501
502 //______________________________________________________________________________
503 void AliPIDResponse::SetITSParametrisation()
504 {
505   //
506   // Set the ITS parametrisation
507   //
508 }
509
510 //______________________________________________________________________________
511 void AliPIDResponse::SetTPCPidResponseMaster()
512 {
513   //
514   // Load the TPC pid response functions from the OADB
515   //
516   //don't load twice for the moment
517    if (fArrPidResponseMaster) return;
518  
519
520   //reset the PID response functions
521   delete fArrPidResponseMaster;
522   fArrPidResponseMaster=0x0;
523   
524   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
525   
526   TFile *f=TFile::Open(fileName.Data());
527   if (f && f->IsOpen() && !f->IsZombie()){
528     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
529   }
530   delete f;
531   
532   if (!fArrPidResponseMaster){
533     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
534     return;
535   }
536   fArrPidResponseMaster->SetOwner();
537 }
538
539 //______________________________________________________________________________
540 void AliPIDResponse::SetTPCParametrisation()
541 {
542   //
543   // Change BB parametrisation for current run
544   //
545   
546   if (fLHCperiod.IsNull()) {
547     AliFatal("No period set, not changing parametrisation");
548     return;
549   }
550   
551   //
552   // Set default parametrisations for data and MC
553   //
554   
555   //data type
556   TString datatype="DATA";
557   //in case of mc fRecoPass is per default 1
558   if (fIsMC) {
559     datatype="MC";
560     fRecoPass=1;
561   }
562   
563   //
564   //reset old splines
565   //
566   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
567     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
568   }
569   
570   //
571   //set the new PID splines
572   //
573   TString period=fLHCperiod;
574   if (fArrPidResponseMaster){
575     TObject *grAll=0x0;
576     //for MC don't use period information
577 //     if (fIsMC) period="[A-Z0-9]*";
578     //for MC use MC period information
579     if (fIsMC) period=fMCperiodTPC;
580 //pattern for the default entry (valid for all particles)
581     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
582     
583     //loop over entries and filter them
584     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
585       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
586       if (responseFunction==0x0) continue;
587       TString responseName=responseFunction->GetName();
588       
589       if (!reg.MatchB(responseName)) continue;
590       
591       TObjArray *arr=reg.MatchS(responseName);
592       TString particleName=arr->At(1)->GetName();
593       delete arr;
594       if (particleName.IsNull()) continue;
595       if (particleName=="ALL") grAll=responseFunction;
596       else {
597         //find particle id
598         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
599           TString particle=AliPID::ParticleName(ispec);
600           particle.ToUpper();
601           if ( particle == particleName ){
602             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
603             fTPCResponse.SetUseDatabase(kTRUE);
604             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
605             break;
606           }
607         }
608       }
609     }
610     
611     //set default response function to all particles which don't have a specific one
612     if (grAll){
613       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
614         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
615           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
616           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
617         }
618       }
619     }
620   }
621   
622   //
623   // Setup resolution parametrisation
624   //
625   
626   //default
627   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
628   
629   if (fRun>=122195){
630     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
631   }
632   if (fArrPidResponseMaster)
633   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
634   
635   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
636 }
637
638 //______________________________________________________________________________
639 void AliPIDResponse::SetTRDPidResponseMaster()
640 {
641   //
642   // Load the TRD pid params and references from the OADB
643   //
644   if(fTRDPIDParams) return;
645   AliOADBContainer contParams; 
646   contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
647   fTRDPIDParams = (TObjArray *)contParams.GetObject(fRun);
648
649   AliOADBContainer contRefs;
650   contRefs.InitFromFile(Form("%s/COMMON/PID/dReferencesLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
651   fTRDPIDReference = (AliTRDPIDReference *)contRefs.GetObject(fRun);
652 }
653
654 //______________________________________________________________________________
655 void AliPIDResponse::InitializeTRDResponse(){
656   //
657   // Set PID Params and references to the TRD PID response
658   // 
659   fTRDResponse.SetPIDParams(fTRDPIDParams);
660   fTRDResponse.Load(fTRDPIDReference);
661   if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
662     fTRDslicesForPID[0] = 0;
663     fTRDslicesForPID[1] = 7;
664   }
665 }
666
667 //_________________________________________________________________________
668 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
669   //
670   // Check whether track is identified as electron under a given electron efficiency hypothesis
671   //
672   Double_t probs[AliPID::kSPECIES];
673   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
674
675   Int_t ntracklets=0;
676   Double_t p = 0;
677   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
678     if(vtrack->GetTRDmomentum(iPl) > 0.){
679       ntracklets++;
680       p = vtrack->GetTRDmomentum(iPl); 
681     }
682   }
683
684   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
685 }
686