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