AliPhysics  master (3d17d9d)
AliAODPidHF.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * *
5  * Author: The ALICE Off-line Project. *
6  * Contributors are mentioned in the code where appropriate. *
7  *
8  * Permission to use, copy, modify and distribute this software and its *
9  * documentation strictly for non-commercial purposes is hereby granted *
10  * without fee, provided that the above copyright notice appears in all *
11  * copies and that both the copyright notice and this permission notice *
12  * appear in the supporting documentation. The authors make no claims *
13  * about the suitability of this software for any purpose. It is *
14  * provided "as is" without express or implied warranty. *
15  * *************************************************************************/
16 
17 /* $Id$ */
18 
24 #include <TCanvas.h>
25 #include <TString.h>
26 #include <TH1F.h>
27 #include <TF1.h>
28 #include <TFile.h>
29 
30 #include "AliAODPidHF.h"
31 #include "AliAODPid.h"
32 #include "AliPID.h"
33 #include "AliPIDResponse.h"
34 #include "AliAODpidUtil.h"
35 #include "AliESDtrack.h"
36 
38 ClassImp(AliAODPidHF);
40 
41 //------------------------------
43 TObject(),
44 fnNSigma(5),
45 fnSigma(0),
46 fTOFSigma(160.),
47 fCutTOFmismatch(0.01),
48 fMinNClustersTPCPID(0),
49 fnPriors(5),
50 fPriors(0),
51 fnPLimit(2),
52 fPLimit(0),
53 fAsym(kFALSE),
54 fTPC(kFALSE),
55 fTOF(kFALSE),
56 fITS(kFALSE),
57 fTRD(kFALSE),
58 fMatch(0),
59 fForceTOFforKaons(kFALSE),
60 fCompat(kFALSE),
61 fPCompatTOF(1.5),
62 fUseAsymTOF(kFALSE),
63 fLownSigmaTOF(-3.),
64 fUpnSigmaTOF(3.),
65 fLownSigmaCompatTOF(-3.),
66 fUpnSigmaCompatTOF(3.),
67 fnNSigmaCompat(2),
68 fnSigmaCompat(0),
69 fMC(kFALSE),
70 fOnePad(kFALSE),
71 fMCLowEn2011(kFALSE),
72 fppLowEn2011(kFALSE),
73 fPbPb(kFALSE),
74 fTOFdecide(kFALSE),
75 fOldPid(kFALSE),
76 fPtThresholdTPC(999999.),
77 fMaxTrackMomForCombinedPID(999999.),
78 fPidResponse(0),
79 fPidCombined(new AliPIDCombined()),
80 fTPCResponse(new AliTPCPIDResponse()),
81 fPriorsH(),
82 fCombDetectors(kTPCTOF),
83 fUseCombined(kFALSE),
84 fDefaultPriors(kTRUE),
85 fApplyNsigmaTPCDataCorr(kFALSE),
86 fMeanNsigmaTPCPionData{},
96 {
104 
105  for(Int_t i=0;i<fnNSigma;i++){
106  fnSigma[i]=0.;
107  }
108  for(Int_t i=0;i<fnPriors;i++){
109  fPriors[i]=0.;
110  }
111  for(Int_t i=0;i<fnPLimit;i++){
112  fPLimit[i]=0.;
113  }
114  for(Int_t i=0;i<fnNSigmaCompat;i++){
115  fnSigmaCompat[i]=3.;
116  }
117  for(Int_t i=0; i<3; i++){ // pi, K, proton
118  fMaxnSigmaCombined[i]=3.;
119  fMinnSigmaTPC[i]=-3;
120  fMaxnSigmaTPC[i]=3;
121  fMinnSigmaTOF[i]=-3;
122  fMaxnSigmaTOF[i]=3;
123  }
124  for (Int_t s=0;s<AliPID::kSPECIES;s++) {
125  for (Int_t d=0;d<4;d++) {
126  fIdBandMin[s][d] = NULL;
127  fIdBandMax[s][d] = NULL;
128  fCompBandMin[s][d] = NULL;
129  fCompBandMax[s][d] = NULL;
130  }
131  }
132 
133  for(int iP=0; iP<=kMaxPBins; iP++) {
134  fPlimitsNsigmaTPCDataCorr[iP] = 0.;
135  }
136  for(int iEta=0; iEta<=kMaxEtaBins; iEta++) {
138  }
139 }
140 //----------------------
142 {
144  if(fPLimit) delete [] fPLimit;
145  if(fnSigma) delete [] fnSigma;
146  if(fPriors) delete [] fPriors;
147  if(fnSigmaCompat) delete [] fnSigmaCompat;
148  delete fPidCombined;
149 
150  delete fTPCResponse;
151  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
152  delete fPriorsH[ispecies];
153  }
154 
155  for (Int_t s=0;s<AliPID::kSPECIES;s++) {
156  for (Int_t d=0;d<4;d++) {
157  delete fIdBandMin[s][d];
158  delete fIdBandMax[s][d];
159  delete fCompBandMin[s][d];
160  delete fCompBandMax[s][d];
161  }
162  }
163 }
164 //------------------------
166 TObject(),
167 fnNSigma(pid.fnNSigma),
168 fnSigma(0),
169 fTOFSigma(pid.fTOFSigma),
172 fnPriors(pid.fnPriors),
173 fPriors(0),
174 fnPLimit(pid.fnPLimit),
175 fPLimit(0),
176 fAsym(pid.fAsym),
177 fTPC(pid.fTPC),
178 fTOF(pid.fTOF),
179 fITS(pid.fITS),
180 fTRD(pid.fTRD),
181 fMatch(pid.fMatch),
183 fCompat(pid.fCompat),
191 fnSigmaCompat(0x0),
192 fMC(pid.fMC),
193 fOnePad(pid.fOnePad),
196 fPbPb(pid.fPbPb),
198 fOldPid(pid.fOldPid),
201 fPidResponse(0x0),
202 fPidCombined(0x0),
203 fTPCResponse(0x0),
210 {
211 
213  for(Int_t i=0;i<fnNSigmaCompat;i++){
214  fnSigmaCompat[i]=pid.fnSigmaCompat[i];
215  }
216  fnSigma = new Double_t[fnNSigma];
217  for(Int_t i=0;i<fnNSigma;i++){
218  fnSigma[i]=pid.fnSigma[i];
219  }
220  fPriors = new Double_t[fnPriors];
221  for(Int_t i=0;i<fnPriors;i++){
222  fPriors[i]=pid.fPriors[i];
223  }
224  fPLimit = new Double_t[fnPLimit];
225  for(Int_t i=0;i<fnPLimit;i++){
226  fPLimit[i]=pid.fPLimit[i];
227  }
228  for(Int_t i=0;i<AliPID::kSPECIES;i++){
229  fPriorsH[i] = pid.fPriorsH[i] ? new TH1F(*pid.fPriorsH[i]) : NULL;
230  }
231  for(Int_t i=0; i<3; i++){ // pi, K, proton
233  fMinnSigmaTPC[i]=pid.fMinnSigmaTPC[i];
234  fMaxnSigmaTPC[i]=pid.fMaxnSigmaTPC[i];
235  fMinnSigmaTOF[i]=pid.fMinnSigmaTOF[i];
236  fMaxnSigmaTOF[i]=pid.fMaxnSigmaTOF[i];
237  }
238 
239  // if(pid.fTPCResponse) fTPCResponse = new AliTPCPIDResponse(*(pid.fTPCResponse));
240  fTPCResponse = new AliTPCPIDResponse();
241  SetBetheBloch();
242  fPidCombined = new AliPIDCombined();
243  //fPidResponse = new AliPIDResponse(*(pid.fPidResponse));
244  //fPidCombined = new AliPIDCombined(*(pid.fPidCombined));
245 
246  //Copy bands
247  for (Int_t s=0;s<AliPID::kSPECIES;s++) {
248  for (Int_t d=0;d<4;d++) {
249  fIdBandMin[s][d] = pid.fIdBandMin[s][d] ? new TF1(*pid.fIdBandMin[s][d]) : NULL;
250  fIdBandMax[s][d] = pid.fIdBandMax[s][d] ? new TF1(*pid.fIdBandMax[s][d]) : NULL;
251  fCompBandMin[s][d] = pid.fCompBandMin[s][d] ? new TF1(*pid.fCompBandMin[s][d]) : NULL;
252  fCompBandMax[s][d] = pid.fCompBandMax[s][d] ? new TF1(*pid.fCompBandMax[s][d]) : NULL;
253  }
254  }
255 
256  for(Int_t iBin=0; iBin<fNPbinsNsigmaTPCDataCorr; iBin++) {
257  for(Int_t iEta=0; iEta<fNEtabinsNsigmaTPCDataCorr; iEta++) {
258  fMeanNsigmaTPCPionData[iEta][iBin] = pid.fMeanNsigmaTPCPionData[iEta][iBin];
259  fMeanNsigmaTPCKaonData[iEta][iBin] = pid.fMeanNsigmaTPCKaonData[iEta][iBin];
260  fMeanNsigmaTPCProtonData[iEta][iBin] = pid.fMeanNsigmaTPCProtonData[iEta][iBin];
261  fSigmaNsigmaTPCPionData[iEta][iBin] = pid.fSigmaNsigmaTPCPionData[iEta][iBin];
262  fSigmaNsigmaTPCKaonData[iEta][iBin] = pid.fSigmaNsigmaTPCKaonData[iEta][iBin];
263  fSigmaNsigmaTPCProtonData[iEta][iBin] = pid.fSigmaNsigmaTPCProtonData[iEta][iBin];
264  }
266  }
268  for(int iEta=0; iEta<=fNEtabinsNsigmaTPCDataCorr; iEta++) {
270  }
271 }
272 //----------------------
273 Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
275  Int_t specie=-1;
276  if(detector.Contains("ITS")) return ApplyPidITSRaw(track,specie);
277  if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,specie);
278  if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,specie);
279 
280  return specie;
281 
282 }
283 //---------------------------
284 Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
286  Int_t specie=0;
287 
288  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
289  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
290  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);
291 
292  if(specie==3) return kTRUE;
293  return kFALSE;
294 }
295 //---------------------------
296 Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
298 
299  Int_t specie=0;
300 
301  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
302  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
303  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);
304 
305  if(specie==2) return kTRUE;
306  return kFALSE;
307 }
308 //---------------------------
309 Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
311 
312  Int_t specie=0;
313  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
314  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4);
315  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);
316 
317  if(specie==4) return kTRUE;
318 
319  return kFALSE;
320 }
321 //--------------------------
322 Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
324 
325  Int_t specie=-1;
326  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,0);
327  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,0);
328  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,0);
329 
330  if(specie==0) return kTRUE;
331 
332  return kFALSE;
333 }
334 //--------------------------
335 Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
337 
338  Double_t nsigma=-999.;
339  Int_t pid=-1;
340 
341  if(specie<0){ // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
342  Double_t nsigmaMin=999.;
343  for(Int_t ipart=0;ipart<5;ipart++){
344  if(GetnSigmaTPC(track,ipart,nsigma)==1){
345  nsigma=TMath::Abs(nsigma);
346  if((nsigma<nsigmaMin) && (nsigma<fnSigma[0])) {
347  pid=ipart;
348  nsigmaMin=nsigma;
349  }
350  }
351  }
352  }else{ // asks only for one particle specie
353  if(GetnSigmaTPC(track,specie,nsigma)==1){
354  nsigma=TMath::Abs(nsigma);
355  if (nsigma>fnSigma[0]) pid=-1;
356  else pid=specie;
357  }
358  }
359 
360  return pid;
361 }
362 //----------------------------
363 Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
365 
366  Double_t nsigma=-999.;
367  Int_t pid=-1;
368 
369  if(specie<0){ // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
370  Double_t nsigmaMin=999.;
371  for(Int_t ipart=0;ipart<5;ipart++){
372  if(GetnSigmaITS(track,ipart,nsigma)==1){
373  nsigma=TMath::Abs(nsigma);
374  if((nsigma<nsigmaMin) && (nsigma<fnSigma[4])) {
375  pid=ipart;
376  nsigmaMin=nsigma;
377  }
378  }
379  }
380  }else{ // asks only for one particle specie
381  if(GetnSigmaITS(track,specie,nsigma)==1){
382  nsigma=TMath::Abs(nsigma);
383  if (nsigma>fnSigma[4]) pid=-1;
384  else pid=specie;
385  }
386  }
387 
388  return pid;
389 }
390 //----------------------------
391 Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
393 
394  Double_t nsigma=-999.;
395  Int_t pid=-1;
396 
397  if(specie<0){
398  Double_t nsigmaMin=999.;
399  for(Int_t ipart=0;ipart<5;ipart++){
400  if(GetnSigmaTOF(track,ipart,nsigma)==1){
401  nsigma=TMath::Abs(nsigma);
402  if((nsigma<nsigmaMin)&& (nsigma<fnSigma[3])){
403  pid=ipart;
404  nsigmaMin=nsigma;
405  }
406  }
407  }
408  }else{ // asks only for one particle specie
409  Double_t nSigmaMin,nSigmaMax;
410  if(fUseAsymTOF){
411  nSigmaMin=fLownSigmaTOF;
412  nSigmaMax=fUpnSigmaTOF;
413  }else{
414  nSigmaMin=-fnSigma[3];
415  nSigmaMax=fnSigma[3];
416  }
417  if(GetnSigmaTOF(track,specie,nsigma)==1){
418  if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
419  else pid=specie;
420  }
421  }
422  return pid;
423 }
424 //----------------------------
425 Int_t AliAODPidHF::ApplyTOFCompatibilityBand(AliAODTrack *track,Int_t specie) const{
427 
428  if(specie<0) return -1;
429  Double_t nsigma=-999.;
430  Int_t pid=-1;
431 
432  Double_t nSigmaMin,nSigmaMax;
433  if(fUseAsymTOF){
434  nSigmaMin=fLownSigmaCompatTOF;
435  nSigmaMax=fUpnSigmaCompatTOF;
436  }else{
437  nSigmaMin=-fnSigmaCompat[1];
438  nSigmaMax=fnSigmaCompat[1];
439  }
440  if(GetnSigmaTOF(track,specie,nsigma)==1){
441  if(nsigma<nSigmaMin || nsigma>nSigmaMax) pid=-1;
442  else pid=specie;
443  }
444  return pid;
445 }
446 //------------------------------
447 void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
449 
450  const Double_t *pid=track->PID();
451  Float_t max=0.;
452  Int_t k=-1;
453  for (Int_t i=0; i<10; i++) {
454  if (pid[i]>max) {k=i; max=pid[i];}
455  }
456 
457  if(k==2) type[0]=kTRUE;
458  if(k==3) type[1]=kTRUE;
459  if(k==4) type[2]=kTRUE;
460 
461  return;
462 }
463 //--------------------------------
464 Bool_t AliAODPidHF::CheckITSPIDStatus(AliAODTrack *track) const{
466  AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kITS,track);
467  if (status != AliPIDResponse::kDetPidOk) return kFALSE;
468  return kTRUE;
469 }
470 //--------------------------------
471 Bool_t AliAODPidHF::CheckTPCPIDStatus(AliAODTrack *track) const{
473  AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTPC,track);
474  if (status != AliPIDResponse::kDetPidOk) return kFALSE;
475  UInt_t nclsTPCPID = track->GetTPCsignalN();
476  if(nclsTPCPID<fMinNClustersTPCPID) return kFALSE;
477  return kTRUE;
478 }
479 //--------------------------------
480 Bool_t AliAODPidHF::CheckTOFPIDStatus(AliAODTrack *track) const{
482  AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTOF,track);
483  if (status != AliPIDResponse::kDetPidOk) return kFALSE;
484  Float_t probMis = fPidResponse->GetTOFMismatchProbability(track);
485  if (probMis > fCutTOFmismatch) return kFALSE;
486  return kTRUE;
487 }
488 //--------------------------------
489 Bool_t AliAODPidHF::CheckTRDPIDStatus(AliAODTrack *track) const{
491  AliPIDResponse::EDetPidStatus status = fPidResponse->CheckPIDStatus(AliPIDResponse::kTRD,track);
492  if (status != AliPIDResponse::kDetPidOk) return kFALSE;
493  return kTRUE;
494 }
495 //--------------------------------
496 Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
498  if(detectors.Contains("ITS")) return CheckITSPIDStatus(track);
499  else if(detectors.Contains("TPC")) return CheckTPCPIDStatus(track);
500  else if(detectors.Contains("TOF")) return CheckTOFPIDStatus(track);
501  else if(detectors.Contains("TRD")) return CheckTRDPIDStatus(track);
502  else{
503  AliError("Wrong detector name");
504  return kFALSE;
505  }
506 }
507 //--------------------------------------------
508 Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
510 
511  AliAODPid *pidObj = track->GetDetPid();
512  Double_t mom = pidObj->GetTPCmomentum();
513  if(mom>fPtThresholdTPC) return kTRUE;
514 
516  if(GetnSigmaTPC(track,specie,nsigma)!=1) return kFALSE;
517  nsigma=TMath::Abs(nsigma);
518 
519 
520  if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
521  if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
522  if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;
523 
524  return kFALSE;
525 }
526 //------------------
527 Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track, Int_t specie){
529 
530  Double_t ptrack=track->P();
531  if(ptrack>fMaxTrackMomForCombinedPID) return 1;
532 
533  Bool_t okTPC=CheckTPCPIDStatus(track);
534  if(ptrack>fPtThresholdTPC) okTPC=kFALSE;
535  Bool_t okTOF=CheckTOFPIDStatus(track);
536 
537  if(fMatch==1){
538  //TOF || TPC (a la' Andrea R.)
539  // convention:
540  // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
541  // the method returns the sum of the response of the 2 detectors
542 
543  if(fTPC && fTOF) {
544  if(!okTPC && !okTOF) return 0;
545  }
546 
547  Int_t tTPCinfo=0;
548  if(fTPC && okTPC){
549  tTPCinfo=-1;
550  if(fAsym) {
551  if(TPCRawAsym(track,specie)) tTPCinfo=1;
552  }else{
553  if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
554  }
555  if(fCompat && tTPCinfo<0){
556  Double_t sig0tmp=fnSigma[0];
557  SetSigma(0,fnSigmaCompat[0]);
558  if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=0;
559  SetSigma(0,sig0tmp);
560  }
561  }
562 
563  Int_t tTOFinfo=0;
564  if(fTOF){
565  if(!okTOF && fTPC) return tTPCinfo;
566  tTOFinfo=-1;
567  if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
568  if(fCompat && tTOFinfo>0){
569  if(ptrack>fPCompatTOF) {
570  if(ApplyTOFCompatibilityBand(track,specie)==specie) tTOFinfo=0;
571  }
572  }
573  }
574 
575 
576  if(tTPCinfo+tTOFinfo==0 && fTOFdecide){
577  if(!okTOF) return tTPCinfo;
578  return tTOFinfo;
579  }
580 
581  if(tTPCinfo+tTOFinfo==0 && fITS){
582  if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
583  Int_t tITSinfo = -1;
584  if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
585  return tITSinfo;
586  }
587  return tTPCinfo+tTOFinfo;
588  }
589 
590  if(fMatch==2){
591  //TPC & TOF (a la' Yifei)
592  // convention: -1 = kFALSE, 1 = kTRUE, 0 = not identified
593  Int_t tTPCinfo=0;
594 
595  if(fTPC && okTPC) {
596  tTPCinfo=1;
597  if(fAsym){
598  if(!TPCRawAsym(track,specie)) tTPCinfo=-1;
599  }else{
600  if(ApplyPidTPCRaw(track,specie)!=specie) tTPCinfo=-1;
601  }
602  }
603 
604  Int_t tTOFinfo=1;
605  if(fTOF){
606  if(fTPC && !okTOF) return tTPCinfo;
607  if(ApplyPidTPCRaw(track,specie)!=specie) tTOFinfo=-1;
608  }
609 
610  if(tTOFinfo==1 && tTPCinfo==1) return 1;
611 
612  if(tTPCinfo+tTOFinfo==0 && fITS){
613  if(!CheckITSPIDStatus(track)) return tTPCinfo+tTOFinfo;
614  Int_t tITSinfo = -1;
615  if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
616  return tITSinfo;
617  }
618  return -1;
619  }
620 
621  if(fMatch==3){
622  //TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
623  // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
624  if(fTPC && fTOF) if(!okTPC && !okTOF) return 0;
625 
626 
627  Int_t tTPCinfo=-1;
628  if(ptrack>=fPLimit[0] && ptrack<fPLimit[1] && fTPC) {
629  if(!okTPC) return 0;
630  if(fAsym) {
631  if(TPCRawAsym(track,specie)) tTPCinfo=1;
632  }else{
633  if(ApplyPidTPCRaw(track,specie)==specie) tTPCinfo=1;
634  }
635  return tTPCinfo;
636  }
637 
638  Int_t tTOFinfo=-1;
639  if(ptrack>=fPLimit[1] && fTOF){
640  if(!okTOF) return 0;
641  if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
642  return tTOFinfo;
643  }
644 
645  Int_t tITSinfo=-1;
646  if(ptrack<fPLimit[0] && fITS){
647  if(!CheckITSPIDStatus(track)) return 0;
648  if(ApplyPidITSRaw(track,specie)==specie) tITSinfo=1;
649  return tITSinfo;
650  }
651  }
652 
653  if(fMatch==4 || fMatch==5){
654 
655  // fMatch == 4 ---> "circular cut" in nSigmaTPC, nSimgaTOF plane
656  // ---> nsigmaTPC^2+nsigmaTOF^2 < cut^2
657  // fMatch == 5 ---> "rectangular cut" in nSigmaTPC, nsigmaTOF plane
658  // ---> ns1<nSigmaTPC<NS1 && ns2<nSigmaTOF<NS2
659 
660  Double_t nSigmaTPC=0.;
661  if(okTPC) {
662  nSigmaTPC = fPidResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)specie);
663  if(fApplyNsigmaTPCDataCorr && nSigmaTPC>-990.) {
664  Float_t mean=0., sigma=1.;
665  GetNsigmaTPCMeanSigmaData(mean, sigma, (AliPID::EParticleType)specie, track->GetTPCmomentum(),track->Eta());
666  nSigmaTPC = (nSigmaTPC-mean)/sigma;
667  }
668  if(nSigmaTPC<-990.) nSigmaTPC=0.;
669  }
670  Double_t nSigmaTOF=0.;
671  if(okTOF) {
672  nSigmaTOF=fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)specie);
673  }
674  Int_t iPart=specie-2; //species is 2 for pions,3 for kaons and 4 for protons
675  if(iPart<0 || iPart>2) return -1;
676  if(fMatch==4){
677  Double_t nSigma2=nSigmaTPC*nSigmaTPC+nSigmaTOF*nSigmaTOF;
678  if(nSigma2<fMaxnSigmaCombined[iPart]*fMaxnSigmaCombined[iPart]) return 1;
679  else return -1;
680  }
681  else if(fMatch==5){
682  if(fForceTOFforKaons && iPart==1 && !okTOF) return -1;
683  if((nSigmaTPC>fMinnSigmaTPC[iPart] && nSigmaTPC<fMaxnSigmaTPC[iPart]) &&
684  (nSigmaTOF>fMinnSigmaTOF[iPart] && nSigmaTOF<fMaxnSigmaTOF[iPart])) return 1;
685  else return -1;
686  }
687  }
688 
689  //Asymmetric cuts using user defined bands
690  if (fMatch == 10) {
691  if (fTPC && fTOF && !okTPC && !okTOF) {
692  return 0;
693  }
694 
695  Int_t tTPCinfo = 0;
696  if (fTPC && okTPC) {
697  tTPCinfo = CheckBands((AliPID::EParticleType) specie, AliPIDResponse::kTPC, track);
698  }
699 
700  Int_t tTOFinfo = 0;
701  if (fTOF) {
702  if (!okTOF && fTPC) {
703  return tTPCinfo;
704  }
705  tTOFinfo = CheckBands((AliPID::EParticleType) specie, AliPIDResponse::kTOF, track);
706  }
707 
708 
709  if (tTPCinfo+tTOFinfo == 0 && fTOFdecide) {
710  if (!okTOF) {
711  return tTPCinfo;
712  }
713  return tTOFinfo;
714  }
715 
716  if (tTPCinfo+tTOFinfo == 0 && fITS) {
717  if (!CheckITSPIDStatus(track)) {
718  return tTPCinfo+tTOFinfo;
719  }
720  Int_t tITSinfo = CheckBands((AliPID::EParticleType) specie, AliPIDResponse::kITS, track);
721  return tITSinfo;
722  }
723  return tTPCinfo+tTOFinfo;
724  }
725 
726  return -1;
727 
728 }
729 
730 //--------------------------------------------------------------
731 Int_t AliAODPidHF::MatchTPCTOFMin(AliAODTrack *track, Int_t specie){
733 
734  Bool_t okTPC=CheckTPCPIDStatus(track);
735  Bool_t okTOF=CheckTOFPIDStatus(track);
736 
737  if(fTPC && fTOF){
738  if(!okTPC && !okTOF) return 0;
739  }
740 
741  Int_t pid=-1;
742  Double_t nsigmaTPC[5]={999.,999.,999.,999.,999.};
743  Double_t nsigmaTOF[5]={999.,999.,999.,999.,999.};
744  Double_t nsigmaMin=999.;
745  Double_t nsigma[5]={999.,999.,999.,999.,999.};
746 
747  if(okTPC) {
748  for(Int_t ipart=0;ipart<5;ipart++){
749  if(GetnSigmaTPC(track,ipart,nsigmaTPC[ipart])<1) nsigmaTPC[ipart]=0.;
750  }
751  }else{
752  for(Int_t ipart=0;ipart<5;ipart++){nsigmaTPC[ipart]=0.;}
753  }
754 
755  if(okTOF){
756  for(Int_t ipart=0;ipart<5;ipart++){
757  if(GetnSigmaTOF(track,ipart,nsigmaTOF[ipart])<1) nsigmaTOF[ipart]=0.;
758  }
759  }else{
760  for(Int_t ipart=0;ipart<5;ipart++){nsigmaTOF[ipart]=0.;}
761  }
762 
763  for(Int_t ipart=0;ipart<5;ipart++){
764  nsigma[ipart]=TMath::Sqrt(nsigmaTPC[ipart]*nsigmaTPC[ipart]+nsigmaTOF[ipart]*nsigmaTOF[ipart]);
765  if(nsigma[ipart]<nsigmaMin) {nsigmaMin=nsigma[ipart];pid=ipart;}
766  }
767 
768  if(pid==specie) return 1;
769 
770  return 0;
771 }
772 
773 
774 //----------------------------------
775 Int_t AliAODPidHF::MakeRawPid(AliAODTrack *track, Int_t specie){
777  if(fMatch>0){
778  return MatchTPCTOF(track,specie);
779  }else{
780  if(fTPC && !fTOF && !fITS) {
781  Int_t tTPCres=0;
782  if(!fAsym){
783  tTPCres=ApplyPidTPCRaw(track,specie);
784  if(tTPCres==specie) return 1;
785  else return tTPCres;
786  }else{
787  if(TPCRawAsym(track,specie)) tTPCres=1;
788  else tTPCres=-1;
789  }
790  return tTPCres;
791  }else if(fTOF && !fTPC && !fITS) {
792  Int_t tTOFres=ApplyPidTOFRaw(track,specie);
793  if(tTOFres==specie) return 1;
794  else return tTOFres;
795  }else if(fITS && !fTPC && !fTOF) {
796  Int_t tITSres=ApplyPidITSRaw(track,specie);
797  if(tITSres==specie) return 1;
798  else return tITSres;
799  }else{
800  AliError("You should enable just one detector if you don't want to match");
801  return 0;
802  }
803  }
804 }
805 //--------------------------------------------
806 void AliAODPidHF::GetTPCBetheBlochParams(Double_t alephParameters[5]) const {
808  if(fMC) { // MC
809 
810  if(fPbPb) { // PbPb MC
811 
812  alephParameters[0] = 1.44405/50.;
813  alephParameters[1] = 2.35409e+01;
814  alephParameters[2] = TMath::Exp(-2.90330e+01);
815  alephParameters[3] = 2.10681e+00;
816  alephParameters[4] = 4.62254e+00;
817 
818  } else { // pp MC
819  if(fMCLowEn2011){
820  alephParameters[0]=0.0207667;
821  alephParameters[1]=29.9936;
822  alephParameters[2]=3.87866e-11;
823  alephParameters[3]=2.17291;
824  alephParameters[4]=7.1623;
825  }else if(fOnePad){
826  alephParameters[0]=0.029021;
827  alephParameters[1]=25.4181;
828  alephParameters[2]=4.66596e-08;
829  alephParameters[3]=1.90008;
830  alephParameters[4]=4.63783;
831  }else{
832  alephParameters[0] = 2.15898/50.;
833  alephParameters[1] = 1.75295e+01;
834  alephParameters[2] = 3.40030e-09;
835  alephParameters[3] = 1.96178e+00;
836  alephParameters[4] = 3.91720e+00;
837  }
838  }
839 
840  } else { // Real Data
841 
842  if(fOnePad) { // pp 1-pad (since LHC10d)
843 
844  alephParameters[0] =1.34490e+00/50.;
845  alephParameters[1] = 2.69455e+01;
846  alephParameters[2] = TMath::Exp(-2.97552e+01);
847  alephParameters[3] = 2.35339e+00;
848  alephParameters[4] = 5.98079e+00;
849 
850  } else if(fPbPb) { // PbPb
851 
852  // alephParameters[0] = 1.25202/50.;
853  // alephParameters[1] = 2.74992e+01;
854  // alephParameters[2] = TMath::Exp(-3.31517e+01);
855  // alephParameters[3] = 2.46246;
856  // alephParameters[4] = 6.78938;
857 
858  alephParameters[0] = 5.10207e+00/50.;
859  alephParameters[1] = 7.94982e+00;
860  alephParameters[2] = TMath::Exp(-9.07942e+00);
861  alephParameters[3] = 2.38808e+00;
862  alephParameters[4] = 1.68165e+00;
863 
864  } else if(fppLowEn2011){ // pp low energy
865 
866  alephParameters[0]=0.031642;
867  alephParameters[1]=22.353;
868  alephParameters[2]=4.16239e-12;
869  alephParameters[3]=2.61952;
870  alephParameters[4]=5.76086;
871 
872  } else { // pp no 1-pad (LHC10bc)
873 
874  alephParameters[0] = 0.0283086/0.97;
875  alephParameters[1] = 2.63394e+01;
876  alephParameters[2] = 5.04114e-11;
877  alephParameters[3] = 2.12543e+00;
878  alephParameters[4] = 4.88663e+00;
879 
880  }
881 
882  }
883 
884 }
885 
886 //-----------------------
889 
890  Double_t alephParameters[5];
891  GetTPCBetheBlochParams(alephParameters);
892  fTPCResponse->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
893 
894  return;
895 }
896 
897 
898 //--------------------------------------------------------------------------
899 Int_t AliAODPidHF::GetnSigmaITS(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
901 
902 
903  if (!CheckITSPIDStatus(track)) return -1;
904 
905  Double_t nsigmaITS=-999;
906 
907  if (fOldPid) {
908  Double_t mom=track->P();
909  AliAODPid *pidObj = track->GetDetPid();
910  Double_t dedx=pidObj->GetITSsignal();
911 
912  AliITSPIDResponse itsResponse;
913  AliPID::EParticleType type=AliPID::EParticleType(species);
914  nsigmaITS = itsResponse.GetNumberOfSigmas(mom,dedx,type);
915 
916  } // old pid
917  else { // new pid
918 
919  AliPID::EParticleType type=AliPID::EParticleType(species);
920  nsigmaITS = fPidResponse->NumberOfSigmasITS(track,type);
921 
922  } //new pid
923 
924  nsigma = nsigmaITS;
925 
926  return 1;
927 
928 }
929 //--------------------------------------------------------------------------
930 Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsigma) const{
932 
933  if(!CheckTPCPIDStatus(track)) return -1;
934 
935  Double_t nsigmaTPC=-999;
936 
937  if(fOldPid){
938  AliAODPid *pidObj = track->GetDetPid();
939  Double_t dedx=pidObj->GetTPCsignal();
940  Double_t mom = pidObj->GetTPCmomentum();
941  if(mom>fPtThresholdTPC) return -2;
942  UShort_t nTPCClus=pidObj->GetTPCsignalN();
943  if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
944  AliPID::EParticleType type=AliPID::EParticleType(species);
945  nsigmaTPC = fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type);
946  nsigma=nsigmaTPC;
947  } else{
948  if(!fPidResponse) return -1;
949  AliPID::EParticleType type=AliPID::EParticleType(species);
950  nsigmaTPC = fPidResponse->NumberOfSigmasTPC(track,type);
951  if(fApplyNsigmaTPCDataCorr && nsigmaTPC>-990.) {
952  Float_t mean=0., sigma=1.;
953  GetNsigmaTPCMeanSigmaData(mean, sigma, type, track->GetTPCmomentum(), track->Eta());
954  nsigmaTPC = (nsigmaTPC-mean)/sigma;
955  }
956  nsigma=nsigmaTPC;
957  }
958  return 1;
959 }
960 
961 //-----------------------------
962 
963 Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
965 
966  if(!CheckTOFPIDStatus(track)) return -1;
967 
968  if(fPidResponse){
969  nsigma = fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)species);
970  return 1;
971  }else{
972  AliFatal("To use TOF PID you need to attach AliPIDResponseTask");
973  nsigma=-999.;
974  return -1;
975  }
976 }
977 
978 //-----------------------
979 Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut, TString detectors) {
981 
982  if (detectors.Contains("ITS")) {
983 
984  AliInfo("Nothing to be done");
985  /*
986  Double_t nsigma=0.;
987  if (GetnSigmaITS(track,labelTrack,nsigma)==1){
988  if(nsigma>nsigmaCut) return kTRUE;
989  }
990  */
991  return kFALSE;
992 
993  } else if (detectors.Contains("TPC")) {
994 
995  Double_t nsigma=0.;
996  if (GetnSigmaTPC(track,labelTrack,nsigma)==1){
997  if(nsigma>nsigmaCut) return kTRUE;
998  }
999  return kFALSE;
1000 
1001  } else if (detectors.Contains("TOF")) {
1002 
1003  Double_t nsigma=0.;
1004  if (GetnSigmaTOF(track,labelTrack,nsigma)==1){
1005  if(nsigma>nsigmaCut) return kTRUE;
1006  }
1007  return kFALSE;
1008 
1009  }
1010  return kFALSE;
1011 
1012 }
1013 //-----------------------
1014 Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
1016 
1017  if(!CheckTOFPIDStatus(track)) return 0;
1018 
1019  Double_t nsigma;
1020  if(GetnSigmaTOF(track,3,nsigma)==1){
1021  if(nsigma>nsigmaK) return kTRUE;
1022  }
1023  return kFALSE;
1024  /* Double_t time[AliPID::kSPECIESN];
1025  Double_t sigmaTOFPid[AliPID::kSPECIES];
1026  AliAODPid *pidObj = track->GetDetPid();
1027  pidObj->GetIntegratedTimes(time);
1028  Double_t sigTOF=pidObj->GetTOFsignal();
1029 
1030  AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
1031  if (event) {
1032  AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
1033  if (tofH && fPidResponse) {
1034  AliTOFPIDResponse TOFres = (AliTOFPIDResponse)fPidResponse->GetTOFResponse();
1035  sigTOF -= TOFres.GetStartTime(track->P());
1036  sigmaTOFPid[3]=TOFres.GetExpectedSigma(track->P(),time[3],AliPID::ParticleMass(3));
1037  }
1038  else pidObj->GetTOFpidResolution(sigmaTOFPid);
1039  } else pidObj->GetTOFpidResolution(sigmaTOFPid);
1040  Double_t sigmaTOFtrack;
1041  if (sigmaTOFPid[3]>0) sigmaTOFtrack=sigmaTOFPid[3];
1042  else sigmaTOFtrack=fTOFSigma; // backward compatibility for old AODs
1043 
1044  if((sigTOF-time[3])>nsigmaK*sigmaTOFtrack)return kTRUE;// K, Pi excluded (->LIKELY A PROTON)
1045 
1046  return kFALSE;
1047  */
1048 }
1049 
1050 //--------------------------------------------------------------------------
1051 void AliAODPidHF::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior){
1052 
1057 
1058  GetPidCombined()->SetPriorDistribution(type,prior);
1059 }
1060 //--------------------------------------------------------------------------
1061 void AliAODPidHF::DrawPrior(AliPID::EParticleType type){
1062 
1065 
1066  new TCanvas();
1067  GetPidCombined()->GetPriorDistribution(type)->Draw();
1068 }
1069 
1070 //-----------------------------
1071 void AliAODPidHF::SetPriors(Double_t *priors, Int_t npriors){
1073  if (fnPriors != npriors) {
1074  if (fPriors) delete fPriors;
1075  fPriors = new Double_t[npriors];
1076  fnPriors = npriors;
1077  }
1078  for(Int_t i = 0; i < fnPriors; i++) fPriors[i] = priors[i];
1079 }
1080 
1081 //-----------------------------
1084  if (fnPLimit != npLim) {
1085  if (fPLimit) delete fPLimit;
1086  fPLimit = new Double_t[npLim];
1087  fnPLimit = npLim;
1088  }
1089  for(Int_t i = 0; i < fnPLimit; i++) fPLimit[i] = plim[i];
1090 }
1091 //-----------------------------
1094 
1095  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
1096  if(fPriorsH[ispecies]) delete fPriorsH[ispecies];
1097  TString nt ="name";
1098  nt+="_prior_";
1099  nt+=AliPID::ParticleName(ispecies);
1100  }
1101  TDirectory *current = gDirectory;
1102  TFile *priorFile=TFile::Open(priorFileName);
1103  if (priorFile) {
1104  TH1F* h3=static_cast<TH1F*>(priorFile->Get("priors3step9"));
1105  TH1F* h2=static_cast<TH1F*>(priorFile->Get("priors2step9"));
1106  TH1F* h1=static_cast<TH1F*>(priorFile->Get("priors1step9"));
1107  current->cd();
1108  fPriorsH[AliPID::kProton] = new TH1F(*h3);
1109  fPriorsH[AliPID::kKaon ] = new TH1F(*h2);
1110  fPriorsH[AliPID::kPion ] = new TH1F(*h1);
1111  priorFile->Close();
1112  delete priorFile;
1113  TF1 *salt=new TF1("salt","1.e-10",0,10);
1114  fPriorsH[AliPID::kProton]->Add(salt);
1115  fPriorsH[AliPID::kKaon ]->Add(salt);
1116  fPriorsH[AliPID::kPion ]->Add(salt);
1117  delete salt;
1118  }
1119 }
1120 //----------------------------------
1123 
1124  fPidCombined->SetSelectedSpecies(AliPID::kSPECIES);
1125  if(!fDefaultPriors){
1126  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
1127  fPidCombined->SetPriorDistribution(static_cast<AliPID::EParticleType>(ispecies),fPriorsH[ispecies]);
1128  }
1129  }else{
1130  fPidCombined->SetDefaultTPCPriors();
1131  }
1132  switch (fCombDetectors){
1133  case kTPCTOF:
1134  fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
1135  break;
1136  case kTPCITS:
1137  fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetITS);
1138  break;
1139  case kTPC:
1140  fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1141  break;
1142  case kTOF:
1143  fPidCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1144  break;
1145  }
1146 }
1147 
1148 
1149 //-----------------------------
1152  printf("Detectors used for PID: ");
1153  if(fITS) printf("ITS ");
1154  if(fTPC) printf("TPC ");
1155  if(fTRD) printf("TRD ");
1156  if(fTOF) printf("TOF ");
1157  printf("\n");
1158  printf("Minimum TPC PID clusters = %d\n",fMinNClustersTPCPID);
1159  printf("Maximum momentum for using TPC PID = %f\n",fPtThresholdTPC);
1160  printf("TOF Mismatch probablility cut = %f\n",fCutTOFmismatch);
1161  printf("Maximum momentum for combined PID TPC PID = %f\n",fMaxTrackMomForCombinedPID);
1162  if(fOldPid){
1163  printf("Use OLD PID");
1164  printf(" fMC = %d\n",fMC);
1165  printf(" fPbPb = %d\n",fPbPb);
1166  printf(" fOnePad = %d\n",fOnePad);
1167  printf(" fMCLowEn2011 = %d\n",fMCLowEn2011);
1168  printf(" fppLowEn2011 = %d\n",fppLowEn2011);
1169  }
1170  printf("--- Matching algorithm = %d ---\n",fMatch);
1171  if(fMatch==1){
1172  if(fITS) printf("nSigmaITS = %.2f\n",fnSigma[4]);
1173  if(fTOF){
1174  printf("nSigmaTOF = %.2f\n",fnSigma[3]);
1175  if(fCompat) printf("Compatibility band at nSigmaTOF=%.2f for p>%.2f\n",fnSigmaCompat[1],fPCompatTOF);
1176  }
1177  if(fTPC){
1178  if(fAsym){
1179  printf("nSigmaTPC:\n");
1180  printf(" pt<%.2f \t nsigmaTPC= %.2f\n",fPLimit[0],fnSigma[0]);
1181  printf(" %.2f<pt<%.2f \t nsigmaTPC= %.2f\n",fPLimit[0],fPLimit[1],fnSigma[1]);
1182  printf(" pt>%.2f \t nsigmaTPC= %.2f\n",fPLimit[1],fnSigma[2]);
1183  }else{
1184  printf("nSigmaTPC = %.2f\n",fnSigma[0]);
1185  }
1186  if(fCompat) printf("Compatibility band at nSigmaTPC=%.2f\n",fnSigmaCompat[0]);
1187  }
1188  }else if(fMatch==4){
1189  printf("Cuts on sqrt(nSigmaTPC^2+nSigmaTOF^2):\n");
1190  printf(" Pions: nSigma = %.2f\n",fMaxnSigmaCombined[0]);
1191  printf(" Kaons: nSigma = %.2f\n",fMaxnSigmaCombined[1]);
1192  printf(" Protons: nSigma = %.2f\n",fMaxnSigmaCombined[2]);
1193  }else if(fMatch==5){
1194  printf("nSigma ranges:\n");
1195  printf(" Pions: %.2f<nSigmaTPC<%.2f %.2f<nSigmaTOF<%.2f\n",
1197  printf(" Kaons: %.2f<nSigmaTPC<%.2f %.2f<nSigmaTOF<%.2f\n",
1199  printf(" Protons: %.2f<nSigmaTPC<%.2f %.2f<nSigmaTOF<%.2f\n",
1201  } else if (fMatch == 10) {
1202  printf("Asymmetric PID using identification/compatibility bands as a function of track momentum p\n");
1203  printf("The following bands are set:\n");
1204  TString species[] = {"electron", "muon", "pion", "kaon", "proton"};
1205  TString detectors[] = {"ITS", "TPC", "TRD", "TOF"};
1206  for (Int_t s=0;s<AliPID::kSPECIES;s++) {
1207  for (Int_t d=0;d<4;d++) {
1208  if (fIdBandMin[s][d] && fIdBandMax[s][d]) {
1209  printf(" Identification band %s %s\n", species[s].Data(), detectors[d].Data());
1210  }
1211  if (fCompBandMin[s][d] && fCompBandMax[s][d]) {
1212  printf(" Compatibility band %s %s\n", species[s].Data(), detectors[d].Data());
1213  }
1214  }
1215  }
1216  }
1217 }
1218 
1219 //------------------
1220 void AliAODPidHF::SetIdBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TH1F *min, TH1F *max) {
1221  Int_t spe = (Int_t) specie;
1222  Int_t det = (Int_t) detector;
1223 
1224  if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1225  AliError("Identification band not set");
1226  return;
1227  }
1228 
1229  TAxis *axis;
1230  HistFunc *histFunc;
1231 
1232  axis = min->GetXaxis();
1233  histFunc = new HistFunc(min);
1234  TF1 *minFunc = new TF1(Form("IdMin_%d_%d", spe, det), *histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1235 
1236  axis = max->GetXaxis();
1237  histFunc = new HistFunc(max);
1238  TF1 *maxFunc = new TF1(Form("IdMax_%d_%d", spe, det), *histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1239 
1240  SetIdBand(specie, detector, minFunc, maxFunc);
1241 }
1242 
1243 //------------------
1244 void AliAODPidHF::SetIdBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TF1 *min, TF1 *max) {
1245  Int_t spe = (Int_t) specie;
1246  Int_t det = (Int_t) detector;
1247 
1248  if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1249  AliError("Identification band not set");
1250  return;
1251  }
1252 
1253  if (fIdBandMin[spe][det]) {
1254  delete fIdBandMin[spe][det];
1255  }
1256  fIdBandMin[spe][det] = new TF1(*min);
1257 
1258  if (fIdBandMax[spe][det]) {
1259  delete fIdBandMax[spe][det];
1260  }
1261  fIdBandMax[spe][det] = new TF1(*max);
1262 }
1263 
1264 //------------------
1265 void AliAODPidHF::SetCompBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TH1F *min, TH1F *max) {
1266  Int_t spe = (Int_t) specie;
1267  Int_t det = (Int_t) detector;
1268 
1269  if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1270  AliError("Compatibility band not set");
1271  return;
1272  }
1273 
1274  TAxis *axis;
1275  HistFunc *histFunc;
1276 
1277  axis = min->GetXaxis();
1278  histFunc = new HistFunc(min);
1279  TF1 *minFunc = new TF1(Form("CompMin_%d_%d", spe, det), *histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1280 
1281  axis = max->GetXaxis();
1282  histFunc = new HistFunc(max);
1283  TF1 *maxFunc = new TF1(Form("CompMax_%d_%d", spe, det), *histFunc, axis->GetBinLowEdge(axis->GetFirst()), axis->GetBinUpEdge(axis->GetLast()), 0, "HistFunc");
1284 
1285  SetCompBand(specie, detector, minFunc, maxFunc);
1286 }
1287 
1288 //------------------
1289 void AliAODPidHF::SetCompBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TF1 *min, TF1 *max) {
1290  Int_t spe = (Int_t) specie;
1291  Int_t det = (Int_t) detector;
1292 
1293  if (spe >= AliPID::kSPECIES || det > 3 || !min || !max) {
1294  AliError("Compatibility band not set");
1295  return;
1296  }
1297 
1298  if (fCompBandMin[spe][det]) {
1299  delete fCompBandMin[spe][det];
1300  }
1301  fCompBandMin[spe][det] = new TF1(*min);
1302 
1303  if (fCompBandMax[spe][det]) {
1304  delete fCompBandMax[spe][det];
1305  }
1306  fCompBandMax[spe][det] = new TF1(*max);
1307 }
1308 
1309 //------------------
1310 Bool_t AliAODPidHF::CheckDetectorPIDStatus(AliPIDResponse::EDetector detector, AliAODTrack* track) {
1311  switch (detector) {
1312  case AliPIDResponse::kITS:
1313  return CheckITSPIDStatus(track);
1314  break;
1315  case AliPIDResponse::kTPC:
1316  return CheckTPCPIDStatus(track);
1317  break;
1318  case AliPIDResponse::kTRD:
1319  return CheckTRDPIDStatus(track);
1320  break;
1321  case AliPIDResponse::kTOF:
1322  return CheckTOFPIDStatus(track);
1323  break;
1324  default:
1325  return kFALSE;
1326  break;
1327  }
1328 }
1329 
1330 //------------------
1331 Float_t AliAODPidHF::NumberOfSigmas(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track) {
1332  switch (detector) {
1333  case AliPIDResponse::kITS:
1334  {
1335  return fPidResponse->NumberOfSigmasITS(track, specie);
1336  break;
1337  }
1338  case AliPIDResponse::kTPC:
1339  {
1340  Double_t nsigmaTPC = fPidResponse->NumberOfSigmasTPC(track, specie);
1341  if(fApplyNsigmaTPCDataCorr && nsigmaTPC>-990.) {
1342  Float_t mean=0., sigma=1.;
1343  GetNsigmaTPCMeanSigmaData(mean, sigma, specie, track->GetTPCmomentum(), track->Eta());
1344  nsigmaTPC = (nsigmaTPC-mean)/sigma;
1345  }
1346  return nsigmaTPC;
1347  break;
1348  }
1349  case AliPIDResponse::kTOF:
1350  {
1351  return fPidResponse->NumberOfSigmasTOF(track, specie);
1352  break;
1353  }
1354  default:
1355  {
1356  return -999.;
1357  break;
1358  }
1359  }
1360 }
1361 
1362 //------------------
1363 Int_t AliAODPidHF::CheckBands(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track) {
1365 
1366  Int_t spe = (Int_t) specie;
1367  Int_t det = (Int_t) detector;
1368 
1369  if (!fPidResponse || spe >= AliPID::kSPECIES) {
1370  return -1;
1371  }
1372 
1373  if (!CheckDetectorPIDStatus(detector, track)) {
1374  return 0;
1375  }
1376 
1377  Double_t P = track->P();
1378 
1379  Float_t nSigma = NumberOfSigmas(specie, detector, track);
1380  Float_t minContent, maxContent;
1381  Bool_t hasAnyBand = kFALSE;
1382 
1383  //Check if within identification band, return 1
1384  TF1 *IdBandMin = fIdBandMin[spe][det];
1385  TF1 *IdBandMax = fIdBandMax[spe][det];
1386 
1387  if (IdBandMin && IdBandMax) {
1388  minContent = IdBandMin->IsInside(&P) ? IdBandMin->Eval(P) : 0;
1389  maxContent = IdBandMax->IsInside(&P) ? IdBandMax->Eval(P) : 0;
1390  if (minContent != 0 || maxContent != 0) {
1391  //At least one identification band is set at this momentum
1392  hasAnyBand = kTRUE;
1393  if ((minContent == 0 || nSigma >= minContent) && (maxContent == 0 || nSigma <= maxContent)) {
1394  return 1;
1395  }
1396  }
1397  }
1398 
1399  //Check if within compatibility band, return 0
1400  TF1 *CompBandMin = fCompBandMin[spe][det];
1401  TF1 *CompBandMax = fCompBandMax[spe][det];
1402 
1403  if (CompBandMin && CompBandMax) {
1404  minContent = CompBandMin->IsInside(&P) ? CompBandMin->Eval(P) : 0;
1405  maxContent = CompBandMax->IsInside(&P) ? CompBandMax->Eval(P) : 0;
1406  if (minContent != 0 || maxContent != 0) {
1407  //At least one compatibility band is set at this momentum
1408  hasAnyBand = kTRUE;
1409  if ((minContent == 0 || nSigma >= minContent) && (maxContent == 0 || nSigma <= maxContent)) {
1410  return 0;
1411  }
1412  }
1413  }
1414 
1415  if (!hasAnyBand) {
1416  //No bands
1417  return 0;
1418  }
1419 
1420  //Bands were set and checked, but no match
1421  return -1;
1422 }
1423 
1424 //------------------
1426  SetMatch(10);
1427  SetTPC(kTRUE);
1428  SetTOF(kTRUE);
1429 
1430  //TPC K: shift by -0.2
1431  TF1 *TPCCompBandMinK = new TF1("TPCCompBandMinK", "[0]", 0, 24); TPCCompBandMinK->SetParameter(0, -3.2);
1432  TF1 *TPCCompBandMaxK = new TF1("TPCCompBandMaxK", "[0]", 0, 24); TPCCompBandMaxK->SetParameter(0, 2.8);
1433  SetCompBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCCompBandMinK, TPCCompBandMaxK);
1434 
1435  TF1 *TPCIdBandMinK = new TF1("TPCIdBandMinK", "[0]", 0, 24); TPCIdBandMinK->SetParameter(0, -2.2);
1436  TF1 *TPCIdBandMaxK = new TF1("TPCIdBandMaxK", "[0]", 0, 24); TPCIdBandMaxK->SetParameter(0, 1.8);
1437  SetIdBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCIdBandMinK, TPCIdBandMaxK);
1438 
1439  //TPC pi: shift by -0.14
1440  TF1 *TPCCompBandMinPi = new TF1("TPCCompBandMinPi", "[0]", 0, 24); TPCCompBandMinPi->SetParameter(0, -3.14);
1441  TF1 *TPCCompBandMaxPi = new TF1("TPCCompBandMaxPi", "[0]", 0, 24); TPCCompBandMaxPi->SetParameter(0, 2.86);
1442  SetCompBand(AliPID::kPion, AliPIDResponse::kTPC, TPCCompBandMinPi, TPCCompBandMaxPi);
1443 
1444  TF1 *TPCIdBandMinPi = new TF1("TPCIdBandMinPi", "[0]", 0, 24); TPCIdBandMinPi->SetParameter(0, -2.14);
1445  TF1 *TPCIdBandMaxPi = new TF1("TPCIdBandMaxPi", "[0]", 0, 24); TPCIdBandMaxPi->SetParameter(0, 1.86);
1446  SetIdBand(AliPID::kPion, AliPIDResponse::kTPC, TPCIdBandMinPi, TPCIdBandMaxPi);
1447 
1448  //TOF K: shift by -0.1
1449  TF1 *TOFCompBandMinK = new TF1("TOFCompBandMinK", "[0]", 2, 24); TOFCompBandMinK->SetParameter(0, -3.1);
1450  TF1 *TOFCompBandMaxK = new TF1("TOFCompBandMaxK", "[0]", 2, 24); TOFCompBandMaxK->SetParameter(0, 2.9);
1451  SetCompBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFCompBandMinK, TOFCompBandMaxK);
1452 
1453  TF1 *TOFIdBandMinK = new TF1("TOFIdBandMinK", "[0]", 0, 2); TOFIdBandMinK->SetParameter(0, -3.1);
1454  TF1 *TOFIdBandMaxK = new TF1("TOFIdBandMaxK", "[0]", 0, 2); TOFIdBandMaxK->SetParameter(0, 2.9);
1455  SetIdBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFIdBandMinK, TOFIdBandMaxK);
1456 
1457  //TOF pi: shift by -0.15
1458  TF1 *TOFCompBandMinPi = new TF1("TOFCompBandMinPi", "[0]", 2, 24); TOFCompBandMinPi->SetParameter(0, -3.15);
1459  TF1 *TOFCompBandMaxPi = new TF1("TOFCompBandMaxPi", "[0]", 2, 24); TOFCompBandMaxPi->SetParameter(0, 2.85);
1460  SetCompBand(AliPID::kPion, AliPIDResponse::kTOF, TOFCompBandMinPi, TOFCompBandMaxPi);
1461 
1462  TF1 *TOFIdBandMinPi = new TF1("TOFIdBandMinPi", "[0]", 0, 2); TOFIdBandMinPi->SetParameter(0, -3.15);
1463  TF1 *TOFIdBandMaxPi = new TF1("TOFIdBandMaxPi", "[0]", 0, 2); TOFIdBandMaxPi->SetParameter(0, 2.85);
1464  SetIdBand(AliPID::kPion, AliPIDResponse::kTOF, TOFIdBandMinPi, TOFIdBandMaxPi);
1465 }
1466 
1467 //------------------
1470 
1471  SetMatch(10);
1472  SetTPC(kTRUE);
1473  SetTOF(kTRUE);
1474 
1475  //TPC K
1476  Double_t TPCIdBandMinKBins[] = {0, 0.4, 0.5, 0.6, 0.9, 24};
1477  TH1F *TPCIdBandMinK = new TH1F("TPCIdBandMinK", "TPC Id Band Min K", 5, TPCIdBandMinKBins);
1478  TPCIdBandMinK->SetBinContent(1, -3); //0 - 0.4
1479  TPCIdBandMinK->SetBinContent(2, -2); //0.4 - 0.5
1480  TPCIdBandMinK->SetBinContent(3, -3); //0.5 - 0.6
1481  TPCIdBandMinK->SetBinContent(4, -2); //0.6 - 0.9
1482  TPCIdBandMinK->SetBinContent(5, -3); //0.9 - 24
1483 
1484  Double_t TPCIdBandMaxKBins[] = {0, 0.6, 0.7, 24};
1485  TH1F *TPCIdBandMaxK = new TH1F("TPCIdBandMaxK", "TPC Id Band Max K", 3, TPCIdBandMaxKBins);
1486  TPCIdBandMaxK->SetBinContent(1, 3); //0 - 0.6
1487  TPCIdBandMaxK->SetBinContent(2, 2); //0.6 - 0.7
1488  TPCIdBandMaxK->SetBinContent(3, 3); //0.7 - 24
1489 
1490  SetIdBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCIdBandMinK, TPCIdBandMaxK);
1493 
1494  //TPC pi
1495  Double_t TPCIdBandMinpiBins[] = {0, 24};
1496  TH1F *TPCIdBandMinpi = new TH1F("TPCIdBandMinpi", "TPC Id Band Min pi", 1, TPCIdBandMinpiBins);
1497  TPCIdBandMinpi->SetBinContent(1, -3); //0 - 24
1498 
1499  Double_t TPCIdBandMaxpiBins[] = {0, 0.7, 0.9, 1.3, 1.4, 24};
1500  TH1F *TPCIdBandMaxpi = new TH1F("TPCIdBandMaxpi", "TPC Id Band Max pi", 5, TPCIdBandMaxpiBins);
1501  TPCIdBandMaxpi->SetBinContent(1, 3); //0 - 0.7
1502  TPCIdBandMaxpi->SetBinContent(2, 2); //0.7 - 0.9
1503  TPCIdBandMaxpi->SetBinContent(3, 3); //0.9 - 1.3
1504  TPCIdBandMaxpi->SetBinContent(4, 2); //1.3 - 1.4
1505  TPCIdBandMaxpi->SetBinContent(5, 3); //1.4 - 24
1506 
1507  SetIdBand(AliPID::kPion, AliPIDResponse::kTPC, TPCIdBandMinpi, TPCIdBandMaxpi);
1510 
1511  //TOF K
1512  TF1 *TOFIdBandMinK = new TF1("TOFIdBandMinK", "[0]", 0, 24); TOFIdBandMinK->SetParameter(0, -3);
1513  TF1 *TOFIdBandMaxK = new TF1("TOFIdBandMaxK", "[0]", 0, 24); TOFIdBandMaxK->SetParameter(0, 3);
1514 
1515  SetIdBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFIdBandMinK, TOFIdBandMaxK);
1516 
1517  //TOF pi
1518  TF1 *TOFIdBandMinPi = new TF1("TOFIdBandMinPi", "[0]", 0, 24); TOFIdBandMinPi->SetParameter(0, -3);
1519  TF1 *TOFIdBandMaxPi = new TF1("TOFIdBandMaxPi", "[0]", 0, 24); TOFIdBandMaxPi->SetParameter(0, 3);
1520 
1521  SetIdBand(AliPID::kPion, AliPIDResponse::kTOF, TOFIdBandMinPi, TOFIdBandMaxPi);
1522 }
1523 
1524 //------------------
1527 
1528  SetMatch(10);
1529  SetTPC(kTRUE);
1530  SetTOF(kTRUE);
1531 
1532  //TPC K
1533  TF1 *TPCCompBandMinK = new TF1("TPCCompBandMinK", "[0]", 0, 24); TPCCompBandMinK->SetParameter(0, -3);
1534  TF1 *TPCCompBandMaxK = new TF1("TPCCompBandMaxK", "[0]", 0, 24); TPCCompBandMaxK->SetParameter(0, 3);
1535 
1536  SetCompBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCCompBandMinK, TPCCompBandMaxK);
1537 
1538  Double_t TPCIdBandMinKBins[6] = {0, 0.45, 0.55, 0.7, 1.1, 24};
1539  TH1F *TPCIdBandMinK = new TH1F("TPCIdBandMinK", "TPC Id Band Min K", 5, TPCIdBandMinKBins);
1540  TPCIdBandMinK->SetBinContent(1, -2); //0-0.45
1541  TPCIdBandMinK->SetBinContent(2, -1); //0.45-0.55
1542  TPCIdBandMinK->SetBinContent(3, -2); //0.55-0.7
1543  TPCIdBandMinK->SetBinContent(4, -1); //0.7-1.1
1544  TPCIdBandMinK->SetBinContent(5, -2); //1.1-24
1545 
1546  Double_t TPCIdBandMaxKBins[4] = {0, 0.5, 0.7, 24};
1547  TH1F *TPCIdBandMaxK = new TH1F("TPCIdBandMaxK", "TPC Id Band Max K", 3, TPCIdBandMaxKBins);
1548  TPCIdBandMaxK->SetBinContent(1, 2); //0-0.5
1549  TPCIdBandMaxK->SetBinContent(2, 1); //0.5-0.7
1550  TPCIdBandMaxK->SetBinContent(3, 2); //0.7-24
1551 
1552  SetIdBand(AliPID::kKaon, AliPIDResponse::kTPC, TPCIdBandMinK, TPCIdBandMaxK);
1555 
1556  //TPC pi
1557  TF1 *TPCCompBandMinpi = new TF1("TPCCompBandMinpi", "[0]", 0, 24); TPCCompBandMinpi->SetParameter(0, -3);
1558  TF1 *TPCCompBandMaxpi = new TF1("TPCCompBandMaxpi", "[0]", 0, 24); TPCCompBandMaxpi->SetParameter(0, 3);
1559 
1560  SetCompBand(AliPID::kPion, AliPIDResponse::kTPC, TPCCompBandMinpi, TPCCompBandMaxpi);
1561 
1562  Double_t TPCIdBandMinpiBins[2] = {0, 24};
1563  TH1F *TPCIdBandMinpi = new TH1F("TPCIdBandMinpi", "TPC Id Band Min pi", 1, TPCIdBandMinpiBins);
1564  TPCIdBandMinpi->SetBinContent(1, -2); //0-24
1565 
1566  Double_t TPCIdBandMaxpiBins[4] = {0, 0.7, 1.7, 24};
1567  TH1F *TPCIdBandMaxpi = new TH1F("TPCIdBandMaxpi", "TPC Id Band Max pi", 3, TPCIdBandMaxpiBins);
1568  TPCIdBandMaxpi->SetBinContent(1, 2); //0-0.7
1569  TPCIdBandMaxpi->SetBinContent(2, 1); //0.7-1.7
1570  TPCIdBandMaxpi->SetBinContent(3, 2); //1.7-24
1571 
1572  SetIdBand(AliPID::kPion, AliPIDResponse::kTPC, TPCIdBandMinpi, TPCIdBandMaxpi);
1575 
1576  //TOF K
1577  TF1 *TOFCompBandMinK = new TF1("TOFCompBandMinK", "[0]", 2, 24); TOFCompBandMinK->SetParameter(0, -3);
1578  TF1 *TOFCompBandMaxK = new TF1("TOFCompBandMaxK", "[0]", 2, 24); TOFCompBandMaxK->SetParameter(0, 3);
1579 
1580  SetCompBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFCompBandMinK, TOFCompBandMaxK);
1581 
1582  TF1 *TOFIdBandMinK = new TF1("TOFIdBandMinK", "[0]", 0, 2); TOFIdBandMinK->SetParameter(0, -3);
1583  TF1 *TOFIdBandMaxK = new TF1("TOFIdBandMaxK", "[0]", 0, 2); TOFIdBandMaxK->SetParameter(0, 3);
1584 
1585  SetIdBand(AliPID::kKaon, AliPIDResponse::kTOF, TOFIdBandMinK, TOFIdBandMaxK);
1586 
1587  //TOF pi
1588  TF1 *TOFCompBandMinpi = new TF1("TOFCompBandMinpi", "[0]", 2, 24); TOFCompBandMinpi->SetParameter(0, -3);
1589  TF1 *TOFCompBandMaxpi = new TF1("TOFCompBandMaxpi", "[0]", 2, 24); TOFCompBandMaxpi->SetParameter(0, 3);
1590 
1591  SetCompBand(AliPID::kPion, AliPIDResponse::kTOF, TOFCompBandMinpi, TOFCompBandMaxpi);
1592 
1593  TF1 *TOFIdBandMinpi = new TF1("TOFIdBandMinpi", "[0]", 0, 2); TOFIdBandMinpi->SetParameter(0, -3);
1594  TF1 *TOFIdBandMaxpi = new TF1("TOFIdBandMaxpi", "[0]", 0, 2); TOFIdBandMaxpi->SetParameter(0, 3);
1595 
1596  SetIdBand(AliPID::kPion, AliPIDResponse::kTOF, TOFIdBandMinpi, TOFIdBandMaxpi);
1597 }
1598 
1599 //------------------
1601 
1602  fApplyNsigmaTPCDataCorr = kTRUE;
1604 }
1605 
1606 //------------------
1607 void AliAODPidHF::GetNsigmaTPCMeanSigmaData(Float_t &mean, Float_t &sigma, AliPID::EParticleType species, Float_t pTPC, Float_t eta) const {
1608 
1609  Int_t bin = TMath::BinarySearch(fNPbinsNsigmaTPCDataCorr,fPlimitsNsigmaTPCDataCorr,pTPC);
1610  if(bin<0) bin=0; //underflow --> equal to min value
1611  else if(bin>fNPbinsNsigmaTPCDataCorr-1) bin=fNPbinsNsigmaTPCDataCorr-1; //overflow --> equal to max value
1612 
1613  Int_t etabin = TMath::BinarySearch(fNEtabinsNsigmaTPCDataCorr,fEtalimitsNsigmaTPCDataCorr,TMath::Abs(eta));
1614  if(etabin<0) etabin=0; //underflow --> equal to min value
1615  else if(etabin>fNEtabinsNsigmaTPCDataCorr-1) etabin=fNEtabinsNsigmaTPCDataCorr-1; //overflow --> equal to max value
1616 
1617  switch(species) {
1618  case AliPID::kPion:
1619  {
1620  mean = fMeanNsigmaTPCPionData[etabin][bin];
1621  sigma = fSigmaNsigmaTPCPionData[etabin][bin];
1622  break;
1623  }
1624  case AliPID::kKaon:
1625  {
1626  mean = fMeanNsigmaTPCKaonData[etabin][bin];
1627  sigma = fSigmaNsigmaTPCKaonData[etabin][bin];
1628  break;
1629  }
1630  case AliPID::kProton:
1631  {
1632  mean = fMeanNsigmaTPCProtonData[etabin][bin];
1633  sigma = fSigmaNsigmaTPCProtonData[etabin][bin];
1634  break;
1635  }
1636  default:
1637  {
1638  mean = 0.;
1639  sigma = 1.;
1640  break;
1641  }
1642  }
1643 }
1644 
1645 //___________________________________________________________________________________//
1646 void AliAODPidHF::SetNsigmaTPCDataDrivenCorrection(Int_t run, Int_t system, Int_t &nPbins, Float_t Plims[kMaxPBins+1], Int_t &nEtabins, Float_t absEtalims[kMaxEtaBins+1],
1647  vector<vector<Float_t> > &meanNsigmaTPCpion, vector<vector<Float_t> > &meanNsigmaTPCkaon, vector<vector<Float_t> > &meanNsigmaTPCproton,
1648  vector<vector<Float_t> > &sigmaNsigmaTPCpion, vector<vector<Float_t> > &sigmaNsigmaTPCkaon, vector<vector<Float_t> > &sigmaNsigmaTPCproton)
1649  {
1650 
1651  meanNsigmaTPCpion.resize(kMaxEtaBins,vector<Float_t>(kMaxPBins,0.));
1652  meanNsigmaTPCkaon.resize(kMaxEtaBins,vector<Float_t>(kMaxPBins,0.));
1653  meanNsigmaTPCproton.resize(kMaxEtaBins,vector<Float_t>(kMaxPBins,0.));
1654  sigmaNsigmaTPCpion.resize(kMaxEtaBins,vector<Float_t>(kMaxPBins,1.));
1655  sigmaNsigmaTPCkaon.resize(kMaxEtaBins,vector<Float_t>(kMaxPBins,1.));
1656  sigmaNsigmaTPCproton.resize(kMaxEtaBins,vector<Float_t>(kMaxPBins,1.));
1657 
1658  if(run>=295585 && run<=296623 && system==kPbPb010) { //LHC18q 0-10%
1659  nPbins = 8;
1660  vector<Float_t> pTPClims = {0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
1661  nEtabins = 5;
1662  vector<Float_t> absetalims = {0.,0.1,0.2,0.4,0.6,0.8};
1663 
1664  vector<vector<Float_t> > meanPion = { {-0.656082, -0.604754, -0.63195, -0.669819, -0.708323, -0.746162, -0.800557, -0.893548},
1665  {-0.711512, -0.686848, -0.711177, -0.752377, -0.781678, -0.838341, -0.80682, -0.917012},
1666  {-0.650884, -0.706274, -0.752449, -0.798673, -0.835543, -0.859084, -0.823375, -0.854207},
1667  {-0.479705, -0.700093, -0.815494, -0.895986, -0.958185, -0.980633, -0.967819, -0.996364},
1668  {-0.264033, -0.637537, -0.828537, -0.952223, -1.06091, -1.10456, -1.09789, -1.15579} };
1669 
1670  vector<vector<Float_t> > meanKaon = { {-0.48114, -0.672897, -0.625657, -0.776678, -0.786824, -0.708909, -0.822472, -0.491422},
1671  {-0.432004, -0.71966, -0.708949, -0.870034, -0.856239, -0.825942, -0.871391, -1.17962},
1672  {-0.336167, -0.71892, -0.735327, -0.93073, -0.864011, -0.891611, -0.924604, -0.735026},
1673  {-0.391054, -0.716122, -0.796868, -1.0208, -0.984637, -0.998813, -1.01377, -1.06832},
1674  {-0.551251, -0.696801, -0.815058, -1.05691, -1.06688, -1.06648, -1.07023, -1.05183} };
1675 
1676  vector<vector<Float_t> > meanProton = { {-0.200581, -0.16751, -0.451043, -0.952266, -0.852847, -0.760682, -0.676723, -0.603716},
1677  {-0.123522, -0.128086, -0.444873, -0.846087, -0.988114, -1.05446, -0.761678, -0.785548},
1678  {-0.100534, -0.1431, -0.448783, -0.8385, -0.843197, -1.05925, -0.878891, -0.69573},
1679  {-0.233023, -0.317509, -0.598837, -0.945108, -1.01043, -1.21354, -0.99634, -0.915479},
1680  {-0.391233, -0.516667, -0.768414, -1.00696, -1.03589, -1.26272, -1.02806, -0.994112} };
1681 
1682  vector<vector<Float_t> > sigmaPion = { {0.986692, 1.01991, 1.00333, 0.986744, 0.981785, 1.04139, 1.0638, 1.09162},
1683  {0.968236, 0.999018, 0.984673, 0.963658, 0.963348, 0.991749, 1.00931, 1.07714},
1684  {0.948544, 0.971808, 0.957514, 0.938766, 0.936994, 0.987149, 0.957657, 0.994133},
1685  {0.92104, 0.931385, 0.916291, 0.896921, 0.890271, 0.926601, 0.891902, 0.905638},
1686  {0.909424, 0.89589, 0.881729, 0.860767, 0.842961, 0.873783, 0.83704, 0.758586} };
1687 
1688  vector<vector<Float_t> > sigmaKaon = { {0.891168, 1.00326, 1.04053, 0.952367, 0.919632, 0.951424, 0.902434, 1.07759},
1689  {0.853805, 0.968212, 1.01839, 0.930692, 0.918841, 0.92428, 0.913715, 0.929855},
1690  {0.830364, 0.911894, 0.987625, 0.912264, 0.939659, 0.906548, 0.893721, 1.00634},
1691  {0.718803, 0.882484, 0.959253, 0.870302, 0.907273, 0.881795, 0.867757, 0.906373},
1692  {0.688955, 0.855596, 0.932222, 0.839553, 0.867639, 0.855212, 0.845177, 0.90448} };
1693 
1694  vector<vector<Float_t> > sigmaProton = { {0.771648, 0.841043, 0.917283, 1.12449, 1.0023, 0.952976, 0.963016, 1.01111},
1695  {0.752951, 0.825488, 0.883897, 1.02998, 1.07061, 0.866346, 0.952794, 0.916068},
1696  {0.72623, 0.799905, 0.860896, 0.996524, 1.02278, 0.866737, 0.834891, 0.988606},
1697  {0.70374, 0.759504, 0.811721, 0.944681, 0.96305, 0.81128, 0.85648, 0.900749},
1698  {0.702538, 0.723393, 0.781419, 0.83867, 0.940137, 0.785817, 0.841202, 0.859564} };
1699 
1700  std::copy(pTPClims.begin(),pTPClims.end(),Plims);
1701  std::copy(absetalims.begin(),absetalims.end(),absEtalims);
1702 
1703  for(int iEta=0; iEta<nEtabins; iEta++) {
1704  meanNsigmaTPCpion[iEta] = meanPion[iEta];
1705  meanNsigmaTPCkaon[iEta] = meanKaon[iEta];
1706  meanNsigmaTPCproton[iEta] = meanProton[iEta];
1707  sigmaNsigmaTPCpion[iEta] = sigmaPion[iEta];
1708  sigmaNsigmaTPCkaon[iEta] = sigmaKaon[iEta];
1709  sigmaNsigmaTPCproton[iEta] = sigmaProton[iEta];
1710  }
1711  }
1712  else if(run>=295585 && run<=296623 && system==kPbPb3050) { //LHC18q 30-50%
1713  nPbins = 8;
1714  vector<Float_t> pTPClims = {0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
1715  nEtabins = 5;
1716  vector<Float_t> absetalims = {0.,0.1,0.2,0.4,0.6,0.8};
1717 
1718  vector<vector<Float_t> > meanPion = { {-0.537046, -0.427744, -0.411915, -0.42242, -0.445157, -0.423209, -0.403354, -0.39011},
1719  {-0.451747, -0.358803, -0.347827, -0.360332, -0.379419, -0.373067, -0.309076, -0.201842},
1720  {-0.37701, -0.342516, -0.352131, -0.375476, -0.397523, -0.380363, -0.348125, -0.334354},
1721  {-0.340843, -0.438716, -0.504516, -0.554322, -0.614602, -0.624125, -0.612949, -0.616095},
1722  {-0.273705, -0.510522, -0.643307, -0.739041, -0.830935, -0.860098, -0.885069, -0.956967} };
1723 
1724  vector<vector<Float_t> > meanKaon = { {-0.226216, -0.427422, -0.414774, -0.562412, -0.521969, -0.467143, -0.414713, -0.330372},
1725  {-0.16309, -0.36387, -0.436445, -0.395585, -0.452207, -0.408419, -0.30696, -0.323571},
1726  {-0.106973, -0.382832, -0.361131, -0.402223, -0.452784, -0.400874, -0.342613, -0.185365},
1727  {-0.252347, -0.480454, -0.500724, -0.573321, -0.642967, -0.616602, -0.546648, -0.482116},
1728  {-0.5916, -0.593699, -0.6563, -0.683526, -0.807421, -0.815922, -0.785543, -0.842433} };
1729 
1730  vector<vector<Float_t> > meanProton = { {0.00677222, 0.0347718, -0.211127, -0.466866, -0.323172, -0.53392, -0.504211, -0.334974},
1731  {0.0935506, 0.0970568, -0.138627, -0.392521, -0.399267, -0.43474, -0.200821, -0.23501},
1732  {0.075394, 0.0609517, -0.170246, -0.409987, -0.420188, -0.448851, -0.267424, -0.313302},
1733  {-0.133011, -0.210294, -0.42431, -0.562203, -0.459603, -0.673718, -0.649959, -0.520375},
1734  {-0.324865, -0.495658, -0.69697, -0.814164, -0.710279, -0.778491, -0.80033, -0.76221} };
1735 
1736  vector<vector<Float_t> > sigmaPion = { {0.915632, 0.93365, 0.932587, 0.931425, 0.922551, 0.92571, 0.881836, 0.796746},
1737  {0.906096, 0.925435, 0.924713, 0.919556, 0.906064, 0.911215, 0.866192, 0.773724},
1738  {0.881485, 0.902513, 0.901312, 0.893701, 0.89239, 0.875522, 0.825261, 0.764695},
1739  {0.835562, 0.850935, 0.853431, 0.846889, 0.834667, 0.819543, 0.798447, 0.774576},
1740  {0.809872, 0.807042, 0.80727, 0.799673, 0.785673, 0.757604, 0.74583, 0.726052} };
1741 
1742  vector<vector<Float_t> > sigmaKaon = { {0.814089, 0.877054, 0.944152, 0.90886, 0.92537, 0.931201, 0.926482, 0.722974},
1743  {0.798753, 0.841944, 0.952407, 0.952185, 0.929863, 0.921988, 0.932327, 0.828574},
1744  {0.733393, 0.821445, 0.903447, 0.928461, 0.917579, 0.921567, 0.923788, 0.727546},
1745  {0.694873, 0.772357, 0.864553, 0.886463, 0.866449, 0.86353, 0.872874, 0.990303},
1746  {0.687564, 0.747428, 0.792902, 0.87729, 0.824179, 0.814195, 0.813793, 0.901867} };
1747 
1748  vector<vector<Float_t> > sigmaProton = { {0.758072, 0.796573, 0.838565, 0.942299, 0.990804, 0.914681, 0.900297, 0.872938},
1749  {0.733353, 0.776217, 0.823004, 0.922272, 1.00522, 0.91665, 0.986521, 1.05648},
1750  {0.715624, 0.75819, 0.79931, 0.897033, 0.972885, 0.920213, 0.942288, 0.978577},
1751  {0.691934, 0.724153, 0.7582, 0.793117, 0.879233, 0.832657, 0.8289, 0.932063},
1752  {0.695097, 0.694242, 0.719957, 0.790272, 0.855222, 0.820038, 0.789541, 0.834146} };
1753 
1754  std::copy(pTPClims.begin(),pTPClims.end(),Plims);
1755  std::copy(absetalims.begin(),absetalims.end(),absEtalims);
1756 
1757  for(int iEta=0; iEta<nEtabins; iEta++) {
1758  meanNsigmaTPCpion[iEta] = meanPion[iEta];
1759  meanNsigmaTPCkaon[iEta] = meanKaon[iEta];
1760  meanNsigmaTPCproton[iEta] = meanProton[iEta];
1761  sigmaNsigmaTPCpion[iEta] = sigmaPion[iEta];
1762  sigmaNsigmaTPCkaon[iEta] = sigmaKaon[iEta];
1763  sigmaNsigmaTPCproton[iEta] = sigmaProton[iEta];
1764  }
1765  }
1766  else if(run>=295585 && run<=296623 && system==kPbPb6080) { //LHC18q 60-80%
1767  nPbins = 8;
1768  vector<Float_t> pTPClims = {0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
1769  nEtabins = 5;
1770  vector<Float_t> absetalims = {0.,0.1,0.2,0.4,0.6,0.8};
1771 
1772  vector<vector<Float_t> > meanPion = { {-0.290461, -0.157362, -0.132967, -0.150766, -0.177387, -0.154474, -0.159506, -0.0819701},
1773  {-0.133763, -0.00498091, 0.0067985, 0.000219383, -0.0387796, -0.0120673, -0.0189252, -0.104004},
1774  {-0.0546859, 0.0146941, 0.012695, -0.0140586, -0.0365883, -0.0455404, -0.0531616, -0.037389},
1775  {-0.0808276, -0.145407, -0.191546, -0.251847, -0.302195, -0.351974, -0.354307, -0.354091},
1776  {-0.093217, -0.288268, -0.403743, -0.500826, -0.589854, -0.641595, -0.693617, -0.700852} };
1777 
1778  vector<vector<Float_t> > meanKaon = { {-0.105712, -0.226045, -0.197616, -0.208206, -0.165878, -0.117968, -0.110323, -0.0885635},
1779  {0.0171305, -0.106295, -0.0665182, 0.00432982, -0.0121926, 0.068422, 0.0699157, 0.0536083},
1780  {0.0668358, -0.115267, -0.049729, 0.00432982, -0.0170053, 0.0328713, 0.0474682, 0.0357199},
1781  {-0.120693, -0.271152, -0.2492, -0.299562, -0.278587, -0.239348, -0.292299, -0.254613},
1782  {-0.485594, -0.44179, -0.427148, -0.531635, -0.522951, -0.552046, -0.679839, -0.666041} };
1783 
1784  vector<vector<Float_t> > meanProton = { {0.0665573, 0.139883, -0.081605, -0.203813, -0.211302, -0.171301, -0.123612, -0.22991},
1785  {0.178373, 0.233109, 0.0432268, -0.059794, -0.0699837, -0.0425294, 0.00809618, 0.176972},
1786  {0.173038, 0.194914, 0.017688, -0.0766737, -0.081089, -0.0359465, 0.0191725, -0.022404},
1787  {-0.0596051, -0.105099, -0.263882, -0.342295, -0.335415, -0.312026, -0.271335, -0.17117},
1788  {-0.281277, -0.420024, -0.588705, -0.614508, -0.618193, -0.559688, -0.506958, -0.587841} };
1789 
1790  vector<vector<Float_t> > sigmaPion = { {0.892435, 0.910563, 0.910262, 0.909437, 0.907557, 0.893655, 0.852062, 0.841},
1791  {0.887739, 0.90651, 0.912569, 0.896137, 0.897385, 0.857293, 0.799653, 0.747155},
1792  {0.861411, 0.890826, 0.889809, 0.891439, 0.87201, 0.847899, 0.830059, 0.813778},
1793  {0.816587, 0.836887, 0.841028, 0.842631, 0.828161, 0.816797, 0.779706, 0.726596},
1794  {0.78294, 0.787287, 0.797276, 0.791521, 0.762976, 0.750954, 0.706446, 0.627991} };
1795 
1796  vector<vector<Float_t> > sigmaKaon = { {0.820386, 0.862568, 0.913401, 0.941527, 0.945742, 0.949472, 0.891267, 0.85767},
1797  {0.777885, 0.833124, 0.907605, 0.99764, 0.965709, 0.976485, 0.946178, 0.860913},
1798  {0.730192, 0.80115, 0.883422, 0.99764, 0.946244, 0.955063, 0.911931, 0.815714},
1799  {0.689248, 0.772107, 0.846482, 0.87357, 0.888549, 0.89338, 0.837796, 0.810673},
1800  {0.638027, 0.725104, 0.79249, 0.819129, 0.831166, 0.809165, 0.74976, 0.811842} };
1801 
1802  vector<vector<Float_t> > sigmaProton = { {0.761836, 0.798676, 0.831566, 0.866324, 0.951735, 0.896439, 0.949077, 0.983275},
1803  {0.728604, 0.764111, 0.811052, 0.848357, 0.872833, 0.936189, 0.937075, 0.949396},
1804  {0.729501, 0.747678, 0.796469, 0.839726, 0.916352, 0.925674, 0.894974, 0.85559},
1805  {0.709612, 0.721265, 0.763333, 0.788841, 0.82415, 0.870806, 0.857862, 0.823096},
1806  {0.697076, 0.692054, 0.711974, 0.749047, 0.78658, 0.79145, 0.776506, 0.57546} };
1807 
1808  std::copy(pTPClims.begin(),pTPClims.end(),Plims);
1809  std::copy(absetalims.begin(),absetalims.end(),absEtalims);
1810 
1811  for(int iEta=0; iEta<nEtabins; iEta++) {
1812  meanNsigmaTPCpion[iEta] = meanPion[iEta];
1813  meanNsigmaTPCkaon[iEta] = meanKaon[iEta];
1814  meanNsigmaTPCproton[iEta] = meanProton[iEta];
1815  sigmaNsigmaTPCpion[iEta] = sigmaPion[iEta];
1816  sigmaNsigmaTPCkaon[iEta] = sigmaKaon[iEta];
1817  sigmaNsigmaTPCproton[iEta] = sigmaProton[iEta];
1818  }
1819  }
1820  else if(run>=296690 && run<=297595 && system==kPbPb010) { //LHC18r 0-10%
1821  nPbins = 8;
1822  vector<Float_t> pTPClims = {0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
1823  nEtabins = 5;
1824  vector<Float_t> absetalims = {0.,0.1,0.2,0.4,0.6,0.8};
1825 
1826  vector<vector<Float_t> > meanPion = { {-0.744242, -0.715713, -0.69801, -0.696472, -0.717912, -0.767909, -0.822175, -0.883157},
1827  {-0.976407, -0.964248, -0.976588, -0.969265, -1.00251, -1.04185, -1.08507, -1.02488},
1828  {-0.938573, -1.02253, -1.0532, -1.06874, -1.09608, -1.11066, -1.07855, -1.06274},
1829  {-0.462091, -0.766549, -0.875959, -0.918783, -0.979887, -0.984493, -0.945828, -0.954307},
1830  {0.154123, -0.361271, -0.568491, -0.667592, -0.782836, -0.751772, -0.732903, -0.749} };
1831 
1832  vector<vector<Float_t> > meanKaon = { {-0.468947, -0.636701, -0.601858, -0.806051, -0.94714, -0.842379, -0.955165, -0.898824},
1833  {-0.588647, -0.883708, -0.894757, -1.09769, -1.11786, -1.08056, -1.15336, -1.44054},
1834  {-0.462369, -0.900829, -0.959231, -1.19209, -1.17182, -1.21788, -1.26831, -1.46315},
1835  {-0.28288, -0.585668, -0.942842, -1.13897, -1.18188, -1.1556, -1.20724, -1.06756},
1836  {-0.0830475, -0.129884, -0.40388, -0.905485, -1.03586, -0.963208, -0.95807, -0.591766} };
1837 
1838  vector<vector<Float_t> > meanProton = { {-0.43448, -0.41261, -0.468653, -0.766399, -0.906529, -0.87423, -0.925983, -0.834281},
1839  {-0.439377, -0.510631, -0.648086, -0.99403, -1.09146, -1.14373, -1.19123, -0.993241},
1840  {-0.400341, -0.54514, -0.649413, -0.979681, -1.22253, -1.27323, -1.27736, -1.12044},
1841  {-0.295184, -0.441872, -0.470702, -0.910766, -1.04581, -1.17824, -1.17277, -0.978326},
1842  {-0.169806, -0.33929, -0.309714, -0.680191, -0.862101, -0.972894, -0.951602, -0.676351} };
1843 
1844  vector<vector<Float_t> > sigmaPion = { {1.19971, 1.20244, 1.19018, 1.18674, 1.19735, 1.27442, 1.31886, 1.35234},
1845  {1.17433, 1.18271, 1.16982, 1.17218, 1.17712, 1.26327, 1.33523, 1.30084},
1846  {1.16732, 1.17134, 1.16173, 1.15583, 1.16336, 1.24219, 1.23569, 1.24103},
1847  {1.1773, 1.16117, 1.14767, 1.14647, 1.19338, 1.22358, 1.21504, 1.20365},
1848  {1.21022, 1.17586, 1.16182, 1.15354, 1.15374, 1.22138, 1.19871, 1.20838} };
1849 
1850  vector<vector<Float_t> > sigmaKaon = { {1.2537, 1.29365, 1.27439, 1.1386, 1.06835, 1.10702, 1.03569, 1.12542},
1851  {1.20352, 1.25335, 1.24745, 1.12229, 1.12846, 1.11836, 1.0746, 1.33407},
1852  {1.17982, 1.21396, 1.23344, 1.13316, 1.15827, 1.10607, 1.06816, 1.14628},
1853  {1.10113, 1.23381, 1.30511, 1.11268, 1.11029, 1.1049, 1.02285, 1.26782},
1854  {1.09653, 1.23731, 1.28769, 1.10684, 1.09593, 1.11015, 1.04836, 1.12603} };
1855 
1856  vector<vector<Float_t> > sigmaProton = { {1.13405, 1.18163, 1.24085, 1.27282, 1.12543, 1.0912, 1.04366, 1.10697},
1857  {1.09801, 1.16854, 1.22617, 1.26278, 1.25383, 1.07191, 1.04581, 1.11811},
1858  {1.11623, 1.17665, 1.22384, 1.24119, 1.23884, 1.04855, 1.05348, 1.11713},
1859  {1.08417, 1.16621, 1.22578, 1.34172, 1.24733, 1.04892, 1.05745, 1.14111},
1860  {1.07421, 1.15563, 1.23173, 1.33796, 1.18478, 1.07793, 1.05178, 1.18215} };
1861 
1862  std::copy(pTPClims.begin(),pTPClims.end(),Plims);
1863  std::copy(absetalims.begin(),absetalims.end(),absEtalims);
1864 
1865  for(int iEta=0; iEta<nEtabins; iEta++) {
1866  meanNsigmaTPCpion[iEta] = meanPion[iEta];
1867  meanNsigmaTPCkaon[iEta] = meanKaon[iEta];
1868  meanNsigmaTPCproton[iEta] = meanProton[iEta];
1869  sigmaNsigmaTPCpion[iEta] = sigmaPion[iEta];
1870  sigmaNsigmaTPCkaon[iEta] = sigmaKaon[iEta];
1871  sigmaNsigmaTPCproton[iEta] = sigmaProton[iEta];
1872  }
1873  }
1874  else if(run>=296690 && run<=297595 && system==kPbPb3050) { //LHC18r 30-50%
1875  nPbins = 8;
1876  vector<Float_t> pTPClims = {0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
1877  nEtabins = 5;
1878  vector<Float_t> absetalims = {0.,0.1,0.2,0.4,0.6,0.8};
1879 
1880  vector<vector<Float_t> > meanPion = { {-0.643992, -0.547379, -0.477054, -0.440121, -0.426498, -0.406638, -0.361916, -0.311107},
1881  {-0.714799, -0.624286, -0.588865, -0.554541, -0.539189, -0.50596, -0.442875, -0.504426},
1882  {-0.641763, -0.623418, -0.607051, -0.580924, -0.573627, -0.534728, -0.520039, -0.457556},
1883  {-0.292997, -0.446284, -0.484493, -0.499929, -0.516507, -0.485561, -0.428434, -0.387974},
1884  {0.151661, -0.179137, -0.303383, -0.371113, -0.418849, -0.413089, -0.380399, -0.442395} };
1885 
1886  vector<vector<Float_t> > meanKaon = { {-0.180292, -0.327542, -0.451762, -0.598542, -0.655667, -0.629728, -0.562158, -0.399934},
1887  {-0.247604, -0.448958, -0.522777, -0.71607, -0.737935, -0.683963, -0.577149, -0.59199},
1888  {-0.189726, -0.48024, -0.575707, -0.742612, -0.765291, -0.709371, -0.607647, -0.28355},
1889  {0.00208301, -0.258596, -0.444814, -0.667385, -0.713196, -0.654455, -0.540685, -0.339524},
1890  {-0.000405731, 0.0385342, -0.162143, -0.523791, -0.612858, -0.553858, -0.45908, -0.401148} };
1891 
1892  vector<vector<Float_t> > meanProton = { {-0.149563, -0.156293, -0.17301, -0.377551, -0.531683, -0.617193, -0.576007, -0.621884},
1893  {-0.14944, -0.204765, -0.255216, -0.438923, -0.675954, -0.637264, -0.640382, -0.585865},
1894  {-0.149585, -0.263429, -0.274904, -0.428442, -0.50586, -0.734513, -0.659482, -0.473013},
1895  {-0.13168, -0.278353, -0.222243, -0.405964, -0.603998, -0.612839, -0.631597, -0.651136},
1896  {-0.108823, -0.280568, -0.193149, -0.250911, -0.53358, -0.55575, -0.545469, -0.421828} };
1897 
1898  vector<vector<Float_t> > sigmaPion = { {1.08474, 1.09912, 1.10605, 1.11466, 1.12312, 1.13388, 1.11915, 0.98612},
1899  {1.0785, 1.09545, 1.10056, 1.10451, 1.11234, 1.12323, 1.10217, 1.14173},
1900  {1.0727, 1.09036, 1.09512, 1.10093, 1.10376, 1.10297, 1.1321, 1.06425},
1901  {1.07348, 1.08058, 1.08902, 1.08921, 1.08795, 1.08588, 1.05973, 1.01606},
1902  {1.08087, 1.06969, 1.08016, 1.08215, 1.08229, 1.07935, 1.04663, 1.02234} };
1903 
1904  vector<vector<Float_t> > sigmaKaon = { {1.10082, 1.14338, 1.17166, 1.08815, 1.07553, 1.05804, 0.98676, 1.09437},
1905  {1.06003, 1.09105, 1.13176, 1.08161, 1.08362, 1.0915, 1.0352, 1.11064},
1906  {1.04461, 1.08936, 1.14854, 1.08121, 1.09336, 1.09623, 1.04665, 0.882711},
1907  {1.01912, 1.09634, 1.15565, 1.08002, 1.09909, 1.09498, 1.05249, 0.999984},
1908  {1.0683, 1.09023, 1.13342, 1.07979, 1.0902, 1.09837, 1.0282, 1.07592} };
1909 
1910  vector<vector<Float_t> > sigmaProton = { {1.06529, 1.10828, 1.14221, 1.18478, 1.10225, 1.07407, 1.08529, 1.14558},
1911  {1.05081, 1.10053, 1.1237, 1.15429, 1.04817, 1.09967, 1.10276, 1.21743},
1912  {1.06485, 1.11335, 1.13392, 1.12825, 1.12658, 1.08174, 1.10395, 1.2337},
1913  {1.05451, 1.11328, 1.13998, 1.17462, 1.19486, 1.09464, 1.07388, 1.07995},
1914  {1.05442, 1.09606, 1.14899, 1.14143, 1.04405, 1.08258, 1.05422, 1.13655} };
1915 
1916  std::copy(pTPClims.begin(),pTPClims.end(),Plims);
1917  std::copy(absetalims.begin(),absetalims.end(),absEtalims);
1918 
1919  for(int iEta=0; iEta<nEtabins; iEta++) {
1920  meanNsigmaTPCpion[iEta] = meanPion[iEta];
1921  meanNsigmaTPCkaon[iEta] = meanKaon[iEta];
1922  meanNsigmaTPCproton[iEta] = meanProton[iEta];
1923  sigmaNsigmaTPCpion[iEta] = sigmaPion[iEta];
1924  sigmaNsigmaTPCkaon[iEta] = sigmaKaon[iEta];
1925  sigmaNsigmaTPCproton[iEta] = sigmaProton[iEta];
1926  }
1927  }
1928  else if(run>=296690 && run<=297595 && system==kPbPb6080) { //LHC18r 60-80%
1929  nPbins = 8;
1930  vector<Float_t> pTPClims = {0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
1931  nEtabins = 5;
1932  vector<Float_t> absetalims = {0.,0.1,0.2,0.4,0.6,0.8};
1933 
1934  vector<vector<Float_t> > meanPion = { {-0.388593, -0.259655, -0.188052, -0.148457, -0.137755, -0.103713, -0.0103896, 0.0183914},
1935  {-0.360916, -0.239814, -0.206334, -0.151416, -0.152129, -0.117832, -0.0512912, -0.0440523},
1936  {-0.275863, -0.225395, -0.194411, -0.162487, -0.166327, -0.137437, -0.0628822, 0.00493976},
1937  {0.0136961, -0.0937519, -0.123901, -0.12913, -0.147936, -0.131946, -0.0580969, -0.140029},
1938  {0.383871, 0.108436, 0.00656927, -0.0673751, -0.108841, -0.136234, -0.0603665, -0.118457} };
1939 
1940  vector<vector<Float_t> > meanKaon = { {-0.0164804, -0.0972051, -0.117233, -0.268106, -0.286532, -0.192977, -0.13506, -0.131241},
1941  {-0.0437434, -0.144134, -0.1716, -0.282154, -0.252251, -0.229581, -0.118788, -0.05944},
1942  {0.0694593, -0.139885, -0.151188, -0.290645, -0.290291, -0.271112, -0.10308, -0.0566368},
1943  {0.183337, 0.0309278, -0.0603532, -0.246455, -0.269791, -0.28694, -0.073309, 0.0979597},
1944  {0.21312, 0.264178, 0.120222, -0.167352, -0.219633, -0.253195, -0.091609, 0.083244} };
1945 
1946  vector<vector<Float_t> > meanProton = { {-0.00831157, -0.0246502, -0.00931054, -0.059935, -0.186943, -0.197498, -0.307994, -0.0702371},
1947  {0.0332978, -0.0489401, -0.03383, -0.0776816, -0.178209, -0.222742, -0.271989, -0.520199},
1948  {-0.0323603, -0.0684328, -0.0158219, -0.073291, -0.180924, -0.258898, -0.308916, -0.418997},
1949  {-0.0288743, -0.135427, 0.00271318, -0.0413093, -0.181906, -0.220064, -0.259092, -0.174753},
1950  {-0.00748328, -0.151846, -0.0157611, 0.0119542, -0.0992955, -0.196113, -0.216777, -0.149337} };
1951 
1952  vector<vector<Float_t> > sigmaPion = { {1.05503, 1.0735, 1.08331, 1.09315, 1.09289, 1.09399, 1.07299, 1.19263},
1953  {1.05972, 1.08279, 1.09055, 1.08779, 1.095, 1.09622, -1.02437, 1.01151},
1954  {1.05442, 1.08032, 1.08152, 1.08671, 1.08594, 1.0849, 1.09464, 1.00198},
1955  {1.05033, 1.06876, 1.07743, 1.08634, 1.08502, 1.09015, 1.03311, 1.02371},
1956  {1.05402, 1.05733, 1.07034, 1.06867, 1.07451, 1.07605, 1.02273, 1.0096} };
1957 
1958  vector<vector<Float_t> > sigmaKaon = { {1.08457, 1.1115, 1.10762, 1.11613, 1.10712, 1.11228, 1.07115, 1.00753},
1959  {1.0578, 1.08034, 1.09357, 1.11391, 1.12987, 1.11409, 1.08436, 1.10904},
1960  {1.03893, 1.06877, 1.10364, 1.11311, 1.13444, 1.10691, 1.09994, 1.07071},
1961  {1.01472, 1.07273, 1.10807, 1.10768, 1.13076, 1.09433, 1.08405, 0.96925},
1962  {1.09291, 1.07712, 1.1035, 1.09887, 1.12159, 1.07188, 1.05105, 0.956145} };
1963 
1964  vector<vector<Float_t> > sigmaProton = { {1.04921, 1.10968, 1.12128, 1.10103, 1.08091, 1.09558, 1.09125, 0.861456},
1965  {1.08841, 1.0943, 1.12142, 1.10862, 1.08891, 1.09479, 1.05673, 1.04939},
1966  {1.11483, 1.10683, 1.11962, 1.09238, 1.07412, 1.10004, 1.08878, 1.18803},
1967  {1.0769, 1.11629, 1.13705, 1.10901, 1.09598, 1.09247, 1.04659, 1.25114},
1968  {1.06329, 1.10681, 1.15275, 1.11122, 1.0892, 1.06841, 1.06915, 1.09393} };
1969 
1970  std::copy(pTPClims.begin(),pTPClims.end(),Plims);
1971  std::copy(absetalims.begin(),absetalims.end(),absEtalims);
1972 
1973  for(int iEta=0; iEta<nEtabins; iEta++) {
1974  meanNsigmaTPCpion[iEta] = meanPion[iEta];
1975  meanNsigmaTPCkaon[iEta] = meanKaon[iEta];
1976  meanNsigmaTPCproton[iEta] = meanProton[iEta];
1977  sigmaNsigmaTPCpion[iEta] = sigmaPion[iEta];
1978  sigmaNsigmaTPCkaon[iEta] = sigmaKaon[iEta];
1979  sigmaNsigmaTPCproton[iEta] = sigmaProton[iEta];
1980  }
1981  }
1982  else { //default: no correction applied
1983  nPbins = 1;
1984  vector<Float_t> pTPClims = {0.,1000.};
1985  nEtabins = 1;
1986  vector<Float_t> absetalims = {0.,2.};
1987 
1988  vector<Float_t> meanPion = {0.};
1989  vector<Float_t> meanKaon = {0.};
1990  vector<Float_t> meanProton = {0.};
1991  vector<Float_t> sigmaPion = {1.};
1992  vector<Float_t> sigmaKaon = {1.};
1993  vector<Float_t> sigmaProton = {1.};
1994  std::copy(pTPClims.begin(),pTPClims.end(),Plims);
1995  std::copy(absetalims.begin(),absetalims.end(),absEtalims);
1996 
1997  meanNsigmaTPCpion[0] = meanPion;
1998  meanNsigmaTPCkaon[0] = meanKaon;
1999  meanNsigmaTPCproton[0] = meanProton;
2000  sigmaNsigmaTPCpion[0] = sigmaPion;
2001  sigmaNsigmaTPCkaon[0] = sigmaKaon;
2002  sigmaNsigmaTPCproton[0] = sigmaProton;
2003  }
2004 }
vector< vector< Float_t > > fSigmaNsigmaTPCPionData
array of NsigmaTPC proton mean in data for different eta bins
Definition: AliAODPidHF.h:320
void SetIdAsymmetricPID()
void SetPriorsHistos(TString priorFileName)
Int_t fnNSigmaCompat
upper nsigma TOF (for fUseAsymTOF)
Definition: AliAODPidHF.h:283
Double_t fLownSigmaTOF
flag for using asymmetrig nSigmaCut in TOF for fMatch==1
Definition: AliAODPidHF.h:279
#define P(T, U, S)
Bool_t fppLowEn2011
MC for low energy MC.
Definition: AliAODPidHF.h:294
double Double_t
Definition: External.C:58
Bool_t fUseCombined
detectors to be involved for combined PID
Definition: AliAODPidHF.h:307
vector< vector< Float_t > > fMeanNsigmaTPCProtonData
array of NsigmaTPC kaon mean in data for different eta bins
Definition: AliAODPidHF.h:319
Bool_t CheckITSPIDStatus(AliAODTrack *track) const
Bool_t IsTOFPiKexcluded(AliAODTrack *track, Double_t nsigmaK)
general method to perform PID using raw signals
Double_t fPtThresholdTPC
old PID method implemented
Definition: AliAODPidHF.h:298
Int_t GetnSigmaTOF(AliAODTrack *track, Int_t species, Double_t &sigma) const
Double_t * fPLimit
limit of p intervals for asimmetric PID: fPLimit<p[0], fPLimit[0]<p<fPLimit[1], p>fPLimit[1] ...
Definition: AliAODPidHF.h:268
Int_t RawSignalPID(AliAODTrack *track, TString detector) const
AliTPCPIDResponse * fTPCResponse
! TPC response
Definition: AliAODPidHF.h:303
Bool_t IsKaonRaw(AliAODTrack *track, TString detector) const
Double_t fMaxnSigmaCombined[3]
Definition: AliAODPidHF.h:286
Bool_t fMCLowEn2011
real data with one pad clusters
Definition: AliAODPidHF.h:293
Float_t NumberOfSigmas(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track)
TF1 * GetIdBandMin(AliPID::EParticleType specie, AliPIDResponse::EDetector detector)
Definition: AliAODPidHF.h:230
Int_t GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const
Bool_t CheckTOFPIDStatus(AliAODTrack *track) const
void GetNsigmaTPCMeanSigmaData(Float_t &mean, Float_t &sigma, AliPID::EParticleType species, Float_t pTPC, Float_t eta) const
Bool_t CheckTPCPIDStatus(AliAODTrack *track) const
void SetPriors(Double_t *priors, Int_t npriors)
static void SetNsigmaTPCDataDrivenCorrection(Int_t run, Int_t system, Int_t &nPbins, Float_t Plims[kMaxPBins+1], Int_t &nEtabins, Float_t absEtalims[kMaxEtaBins+1], vector< vector< Float_t > > &meanNsigmaTPCpion, vector< vector< Float_t > > &meanNsigmaTPCkaon, vector< vector< Float_t > > &meanNsigmaTPCproton, vector< vector< Float_t > > &sigmaNsigmaTPCpion, vector< vector< Float_t > > &sigmaNsigmaTPCkaon, vector< vector< Float_t > > &sigmaNsigmaTPCproton)
TF1 * fIdBandMax[AliPID::kSPECIES][4]
Definition: AliAODPidHF.h:312
Bool_t IsElectronRaw(AliAODTrack *track, TString detector) const
TF1 * fIdBandMin[AliPID::kSPECIES][4]
use default priors for combined PID
Definition: AliAODPidHF.h:311
Bool_t fMC
max. of nSigma range for pi,K,p in TOF (match==5)
Definition: AliAODPidHF.h:291
Bool_t fCompat
force TOF for kaons in mode fMatch==5
Definition: AliAODPidHF.h:276
Int_t fnNSigma
Definition: AliAODPidHF.h:257
Bool_t fITS
switch to include or exclude TOF
Definition: AliAODPidHF.h:272
Int_t ApplyPidTPCRaw(AliAODTrack *track, Int_t specie) const
Double_t fMaxnSigmaTOF[3]
min. of nSigma range for pi,K,p in TOF (match==5)
Definition: AliAODPidHF.h:290
Int_t ApplyPidITSRaw(AliAODTrack *track, Int_t specie) const
void SetCompBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TH1F *min, TH1F *max)
TH1F * fPriorsH[AliPID::kSPECIES]
Definition: AliAODPidHF.h:305
Int_t MatchTPCTOFMin(AliAODTrack *track, Int_t specie)
PID nSigma strategy closer to the Bayesian approach with Max. prob.
void PrintAll() const
Double_t fMaxTrackMomForCombinedPID
pT threshold to use TPC PID
Definition: AliAODPidHF.h:299
Double_t fMinnSigmaTOF[3]
max. of nSigma range for pi,K,p in TPC (match==5)
Definition: AliAODPidHF.h:289
Double_t fMinnSigmaTPC[3]
nSigma cut for pi,K,p (TPC^2+TOF^2)
Definition: AliAODPidHF.h:287
TF1 * fCompBandMin[AliPID::kSPECIES][4]
Definition: AliAODPidHF.h:313
void SetTOF(Bool_t tof)
Definition: AliAODPidHF.h:108
TF1 * GetIdBandMax(AliPID::EParticleType specie, AliPIDResponse::EDetector detector)
Definition: AliAODPidHF.h:231
Int_t fNEtabinsNsigmaTPCDataCorr
array of eta limits for data-driven NsigmaTPC correction
Definition: AliAODPidHF.h:326
Int_t fnPLimit
Definition: AliAODPidHF.h:266
Bool_t fOnePad
MC(kTRUE) or real data (kFALSE, default option)
Definition: AliAODPidHF.h:292
Double_t * sigma
void SetIdBand(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, TH1F *min, TH1F *max)
Assymetric PID using histograms.
Double_t fTOFSigma
Definition: AliAODPidHF.h:260
Int_t CheckBands(AliPID::EParticleType specie, AliPIDResponse::EDetector detector, AliAODTrack *track)
AliPIDResponse * fPidResponse
momentum threshold to use PID
Definition: AliAODPidHF.h:300
void DrawPrior(AliPID::EParticleType type)
Double_t * fPriors
set of priors
Definition: AliAODPidHF.h:265
int Int_t
Definition: External.C:63
virtual ~AliAODPidHF()
unsigned int UInt_t
Definition: External.C:33
AliPIDCombined * GetPidCombined() const
Definition: AliAODPidHF.h:174
float Float_t
Definition: External.C:68
Int_t fNPbinsNsigmaTPCDataCorr
array of p limits for data-driven NsigmaTPC correction
Definition: AliAODPidHF.h:324
void GetTPCBetheBlochParams(Double_t alephParameters[5]) const
Bool_t fOldPid
real data PbPb
Definition: AliAODPidHF.h:297
void SetShiftedAsymmetricPID()
Some suggested asymmetric PID.
Bool_t fDefaultPriors
detectors to be involved for combined PID
Definition: AliAODPidHF.h:308
static const int kMaxPBins
Definition: AliAODPidHF.h:45
Bool_t CheckStatus(AliAODTrack *track, TString detectors) const
Double_t fPCompatTOF
compatibility region : useful only if fMatch=1
Definition: AliAODPidHF.h:277
Bool_t IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut, TString detector)
Bool_t fForceTOFforKaons
switch to combine the info from more detectors: 1 = || , 2 = &, 3 = p region
Definition: AliAODPidHF.h:275
Double_t nsigma
Bool_t fTOF
switch to include or exclude TPC
Definition: AliAODPidHF.h:271
void SetUpCombinedPID()
Int_t MakeRawPid(AliAODTrack *track, Int_t specie)
Int_t fMatch
switch to include or exclude TRD
Definition: AliAODPidHF.h:274
Int_t MatchTPCTOF(AliAODTrack *track, Int_t specie)
void SetIdCompAsymmetricPID()
void SetSigma(Double_t *sigma)
Definition: AliAODPidHF.h:52
Int_t fnPriors
Minimum TPC PID clusters cut.
Definition: AliAODPidHF.h:263
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Bool_t fTOFdecide
real data PbPb
Definition: AliAODPidHF.h:296
Int_t ApplyPidTOFRaw(AliAODTrack *track, Int_t specie) const
void SetPriorDistribution(AliPID::EParticleType type, TH1F *prior)
Double_t fMaxnSigmaTPC[3]
min. of nSigma range for pi,K,p in TPC (match==5)
Definition: AliAODPidHF.h:288
void SetBetheBloch()
UInt_t fMinNClustersTPCPID
Cut of TOF mismatch probability.
Definition: AliAODPidHF.h:262
Double_t fLownSigmaCompatTOF
upper nsigma TOF (for fUseAsymTOF)
Definition: AliAODPidHF.h:281
vector< vector< Float_t > > fMeanNsigmaTPCKaonData
array of NsigmaTPC pion mean in data for different eta bins
Definition: AliAODPidHF.h:318
Int_t GetnSigmaITS(AliAODTrack *track, Int_t species, Double_t &sigma) const
AliPIDCombined * fPidCombined
! combined PID object
Definition: AliAODPidHF.h:301
Bool_t IsPionRaw(AliAODTrack *track, TString detector) const
Bool_t CheckDetectorPIDStatus(AliPIDResponse::EDetector detector, AliAODTrack *track)
Int_t ApplyTOFCompatibilityBand(AliAODTrack *track, Int_t specie) const
Double_t * fnSigma
sigma for the raw signal PID: 0-2 for TPC, 3 for TOF, 4 for ITS
Definition: AliAODPidHF.h:259
Bool_t fTPC
asimmetric PID required
Definition: AliAODPidHF.h:270
Bool_t fUseAsymTOF
compatibility p limit for TOF
Definition: AliAODPidHF.h:278
Double_t fUpnSigmaCompatTOF
lower nsigma TOF (for fUseAsymTOF)
Definition: AliAODPidHF.h:282
Bool_t fTRD
switch to include or exclude ITS
Definition: AliAODPidHF.h:273
vector< vector< Float_t > > fMeanNsigmaTPCPionData
flag to enable data-driven NsigmaTPC correction
Definition: AliAODPidHF.h:317
static const int kMaxEtaBins
Definition: AliAODPidHF.h:44
Float_t fPlimitsNsigmaTPCDataCorr[kMaxPBins+1]
array of NsigmaTPC proton mean in data for different eta bins
Definition: AliAODPidHF.h:323
Float_t fEtalimitsNsigmaTPCDataCorr[kMaxEtaBins+1]
number of p bins for data-driven NsigmaTPC correction
Definition: AliAODPidHF.h:325
unsigned short UShort_t
Definition: External.C:28
void EnableNsigmaTPCDataCorr(Int_t run, Int_t system)
Set Nsigma data-driven correction.
vector< vector< Float_t > > fSigmaNsigmaTPCKaonData
array of NsigmaTPC pion mean in data for different eta bins
Definition: AliAODPidHF.h:321
void SetMatch(Int_t match)
Definition: AliAODPidHF.h:111
void SetPLimit(Double_t *plim, Int_t npLim)
Bool_t TPCRawAsym(AliAODTrack *track, Int_t specie) const
void SetTPC(Bool_t tpc)
Definition: AliAODPidHF.h:107
Bool_t fApplyNsigmaTPCDataCorr
Definition: AliAODPidHF.h:316
bool Bool_t
Definition: External.C:53
Bool_t fAsym
Definition: AliAODPidHF.h:269
Bool_t fPbPb
Data for low energy pp 2011.
Definition: AliAODPidHF.h:295
Double_t fUpnSigmaTOF
lower nsigma TOF (for fUseAsymTOF)
Definition: AliAODPidHF.h:280
ECombDetectors fCombDetectors
priors histos
Definition: AliAODPidHF.h:306
void CombinedProbability(AliAODTrack *track, Bool_t *type) const
Bool_t IsProtonRaw(AliAODTrack *track, TString detector) const
Double_t * fnSigmaCompat
0: n sigma for TPC compatibility band, 1: for TOF
Definition: AliAODPidHF.h:285
vector< vector< Float_t > > fSigmaNsigmaTPCProtonData
array of NsigmaTPC kaon mean in data for different eta bins
Definition: AliAODPidHF.h:322
Double_t fCutTOFmismatch
TOF precision.
Definition: AliAODPidHF.h:261
Bool_t CheckTRDPIDStatus(AliAODTrack *track) const
TF1 * fCompBandMax[AliPID::kSPECIES][4]
Definition: AliAODPidHF.h:314