AliPhysics  ebc57ae (ebc57ae)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskPIDconfig.cxx
Go to the documentation of this file.
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: AliAnalysisTaskPIDconfig.cxx 43811 2014-10-11 Naghmeh Mohammadi $ */
17 
18 #include "TChain.h"
19 #include "TTree.h"
20 #include "TList.h"
21 #include "TMath.h"
22 #include "TObjArray.h"
23 #include "TCanvas.h"
24 #include "TGraphErrors.h"
25 #include "TString.h"
26 #include "TFile.h"
27 #include "TH1F.h"
28 #include "TH2F.h"
29 #include "TH3F.h"
30 #include "TH2D.h"
31 #include "TH3D.h"
32 #include "TArrayF.h"
33 #include "TF1.h"
34 #include "TROOT.h"
35 #include "stdio.h"
36 #include "TCutG.h"
37 #include "TProfile.h"
38 #include "TF2.h"
39 #include "THnSparse.h"
40 
41 
42 
43 #include "AliTHn.h"
44 #include "AliLog.h"
45 #include "AliAnalysisManager.h"
46 #include "AliESDEvent.h"
47 #include "AliAODInputHandler.h"
48 #include "AliAODEvent.h"
49 #include "AliAODTrack.h"
50 #include "AliAODInputHandler.h"
51 #include "AliCollisionGeometry.h"
52 #include "AliGenEventHeader.h"
53 #include "AliAnalysisUtils.h"
54 #include "AliPIDCombined.h"
55 #include "AliAnalysisTask.h"
56 #include "AliAODHandler.h"
57 #include <AliInputEventHandler.h>
58 #include <AliVEventHandler.h>
59 #include <AliVParticle.h>
60 #include <AliVTrack.h>
61 #include <AliTPCPIDResponse.h>
62 #include <AliTOFPIDResponse.h>
64 #include "AliAnalysisTaskSE.h"
65 #include "AliAODPid.h"
66 #include "AliPhysicsSelection.h"
67 #include "AliCentralitySelectionTask.h"
68 #include "AliCentrality.h"
69 #include "AliKFParticle.h"
70 #include "AliKFVertex.h"
71 #include "AliPID.h"
72 #include "AliPIDResponse.h"
73 #include "AliCFContainer.h"
74 #include "AliCFManager.h"
75 #include "AliVEvent.h"
76 #include "AliAODVZERO.h"
77 
78 using std::cout;
79 using std::endl;
80 
81 
83 //ClassImp()
84 //___________________________________________________________________
87 fVevent(0),
88 fESD(0),
89 fAOD(0),
90 fPIDResponse(0),
91 fTriggerSelection(0),
92 fCentralityPercentileMin(0.),
93 fCentralityPercentileMax(5.),
94 fPurityLevel(0.8),
95 fFilterBit(1),
96 fDCAxyCut(-1),
97 fDCAzCut(-1),
98 fData2011(kFALSE),
99 fUseCentrality(kTRUE),
100 fCutTPCmultiplicityOutliersAOD(kFALSE),
101 fPIDcuts(kFALSE),
102 fCentralityEstimator("V0M"),
103 fPurityFunctionsFile(0),
104 fPurityFunctionsList(0),
105 fListQA(0x0),
106 fListQAtpctof(0x0),
107 fListQAInfo(0x0),
108 fhistCentralityPassBefore(0),
109 fhistCentralityPassAfter(0),
110 fNoEvents(0),
111 fpVtxZ(0),
112 fhistDCABefore(0),
113 fhistDCAAfter(0),
114 fhistPhiDistBefore(0),
115 fhistPhiDistAfter(0),
116 fhistEtaDistBefore(0),
117 fhistEtaDistAfter(0),
118 fTPCvsGlobalMultBeforeOutliers(0),
119 fTPCvsGlobalMultAfterOutliers(0),
120 fTPCvsGlobalMultAfter(0),
121 fHistBetavsPTOFbeforePID(0),
122 fHistdEdxvsPTPCbeforePID(0),
123 fhistNsigmaP(0),
124 fhistTPCnSigmavsP(0),
125 fhistTOFnSigmavsP(0),
126 fHistBetavsPTOFafterPID(0),
127 fHistdEdxvsPTPCafterPID(0),
128 fHistBetavsPTOFafterPID_2(0),
129 fHistdEdxvsPTPCafterPID_2(0),
130 fHistBetavsPTOFafterPIDTPCTOF(0),
131 fHistdEdxvsPTPCafterPIDTPCTOF(0),
132 fHistBetavsPTOFafterPIDTPConly(0),
133 fHistdEdxvsPTPCafterPIDTPConly(0),
134 fHistPion_BetavsPTOFafterPIDTPCTOF(0),
135 fHistPion_dEdxvsPTPCafterPIDTPCTOF(0),
136 fHistKaon_BetavsPTOFafterPIDTPCTOF(0),
137 fHistKaon_dEdxvsPTPCafterPIDTPCTOF(0),
138 fHistProton_BetavsPTOFafterPIDTPCTOF(0),
139 fHistProton_dEdxvsPTPCafterPIDTPCTOF(0),
140 fhistPionEtaDistAfter(0),
141 fhistKaonEtaDistAfter(0),
142 fhistProtonEtaDistAfter(0),
143 fSparseAll(0)
144 //fSparseSpecies(0),
145 //fvalueSpecies(0)
146 {
147  for(int i=0;i<4;i++){fvalueAll[i]=0;}
148  for(int i=0;i<180;i++){
149  fPurityFunction[i]=NULL;
150  }
151  //Low momentum nsigma cuts based on Purity>0.7 with TPC info only.
152 
153  for(int i=0;i<6;i++){
154 
155  fLowPtPIDTPCnsigLow_Pion[i] = 0;
156  fLowPtPIDTPCnsigLow_Kaon[i] = 0;
157  fLowPtPIDTPCnsigHigh_Pion[i] =0;
158  fLowPtPIDTPCnsigHigh_Kaon[i] =0;
159  }
160 
161 }
162 
163 
164 //___________________________________________________________________
165 
167 AliAnalysisTaskSE(name),
168 fVevent(0),
169 fESD(0),
170 fAOD(0),
171 fPIDResponse(0),
172 fTriggerSelection(0),
173 fCentralityPercentileMin(0.),
174 fCentralityPercentileMax(5.),
175 fPurityLevel(0.8),
176 fFilterBit(1),
177 fDCAxyCut(-1),
178 fDCAzCut(-1),
179 fData2011(kFALSE),
180 fUseCentrality(kTRUE),
181 fCutTPCmultiplicityOutliersAOD(kFALSE),
182 fPIDcuts(kFALSE),
183 fCentralityEstimator("V0M"),
184 fPurityFunctionsFile(0),
185 fPurityFunctionsList(0),
186 fListQA(0x0),
187 fListQAtpctof(0x0),
188 fListQAInfo(0x0),
189 fhistCentralityPassBefore(0),
190 fhistCentralityPassAfter(0),
191 fNoEvents(0),
192 fpVtxZ(0),
193 fhistDCABefore(0),
194 fhistDCAAfter(0),
195 fhistPhiDistBefore(0),
196 fhistPhiDistAfter(0),
197 fhistEtaDistBefore(0),
198 fhistEtaDistAfter(0),
199 fTPCvsGlobalMultBeforeOutliers(0),
200 fTPCvsGlobalMultAfterOutliers(0),
201 fTPCvsGlobalMultAfter(0),
202 fHistBetavsPTOFbeforePID(0),
203 fHistdEdxvsPTPCbeforePID(0),
204 fhistNsigmaP(0),
205 fhistTPCnSigmavsP(0),
206 fhistTOFnSigmavsP(0),
207 fHistBetavsPTOFafterPID(0),
208 fHistdEdxvsPTPCafterPID(0),
209 fHistBetavsPTOFafterPID_2(0),
210 fHistdEdxvsPTPCafterPID_2(0),
211 fHistBetavsPTOFafterPIDTPCTOF(0),
212 fHistdEdxvsPTPCafterPIDTPCTOF(0),
213 fHistBetavsPTOFafterPIDTPConly(0),
214 fHistdEdxvsPTPCafterPIDTPConly(0),
215 fHistPion_BetavsPTOFafterPIDTPCTOF(0),
216 fHistPion_dEdxvsPTPCafterPIDTPCTOF(0),
217 fHistKaon_BetavsPTOFafterPIDTPCTOF(0),
218 fHistKaon_dEdxvsPTPCafterPIDTPCTOF(0),
219 fHistProton_BetavsPTOFafterPIDTPCTOF(0),
220 fHistProton_dEdxvsPTPCafterPIDTPCTOF(0),
221 fhistPionEtaDistAfter(0),
222 fhistKaonEtaDistAfter(0),
223 fhistProtonEtaDistAfter(0),
224 fSparseAll(0)
225 //fSparseSpecies(0),
226 //fvalueSpecies(0)
227 {
228  for(int i=0;i<4;i++){fvalueAll[i]=0;}
229  //Default Constructor
230  for(int i=0;i<180;i++){
231  fPurityFunction[i]=NULL;
232  }
233 
234  DefineInput(0,TChain::Class());
235  DefineOutput(1,TList::Class());
236 }
237 
238 //_____________________________________________________________________
240 {
241  //Destructor
242 
243  fPurityFunctionsFile->Close();
244  for(int i=0;i<180;i++){
245  delete fPurityFunction[i];
246  }
247 
248  delete fSparseAll;
249 
250 }
251 //______________________________________________________________________
253 {
254  //
255  // Create the output QA objects
256  //
257 
258  AliLog::SetClassDebugLevel("AliAnalysisTaskPIDconfig",10);
259 
260  //input hander
261  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
262  AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
263  if (!inputHandler) {
264  AliFatal("Input handler needed");
265  return; // to shut up coverity
266  }
267 
268  //pid response object
269  fPIDResponse=inputHandler->GetPIDResponse();
270  if (!fPIDResponse) AliError("PIDResponse object was not created");
271 
272  fListQA=new TList;
273  fListQA->SetOwner();
274 
275  fListQAtpctof=new TList;
276  fListQAtpctof->SetOwner();
277  fListQAtpctof->SetName("PID_TPC_TOF");
278 
279  fListQAInfo=new TList;
280  fListQAInfo->SetOwner();
281  fListQAInfo->SetName("Event_Track_Info");
282 
283  fListQA->Add(fListQAtpctof);
284  fListQA->Add(fListQAInfo);
285 
286  SetupTPCTOFqa();
287  SetupEventInfo();
288 
289 
290  PostData(1,fListQA);
291 
292 }
293 //______________________________________________________________________
295  //Main loop
296  //Called for each event
297 
298  // create pointer to event
299  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
300  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
301 
302 
303 
304  if(!(fESD || fAOD)){
305  printf("ERROR: fESD & fAOD not available\n");
306  return;
307  }
308 
309  Int_t ntracks=fAOD->GetNumberOfTracks();
310 
311  fVevent = dynamic_cast<AliVEvent*>(InputEvent());
312  if (!fVevent) {
313  printf("ERROR: fVevent not available\n");
314  return;
315  }
316 
317  TH1F *hNoEvents = (TH1F*)fListQAInfo->At(0);
318  hNoEvents->Fill(1);
319 
320  Double_t centrality = fVevent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
321  TH1F *hCentralityPassBefore = (TH1F*)fListQAInfo->At(1);
322  hCentralityPassBefore->Fill(centrality);
323 
324 
325  Bool_t pass = kFALSE;
326  CheckCentrality(fVevent,centrality,pass);
327 
328  if(!pass){ return;}
329 
330  hNoEvents->Fill(2);
331 
332  const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
333 
334  Double_t pVtxZ = -999;
335  pVtxZ = pVtx->GetZ();
336 
337  if(TMath::Abs(pVtxZ)>10) return;
338 
339  TH1F *histpVtxZ = (TH1F*)fListQAInfo->At(3);
340 
341  hNoEvents->Fill(3);
342 
343  if(histpVtxZ) histpVtxZ->Fill(pVtxZ);
344 
345  if(ntracks<2) return;
346 
347  // if(!pass) return;
348 
349  TH2F *HistTPCvsGlobalMultBeforeOutliers = (TH2F*)fListQAInfo->At(4);
350 
351  TH2F *HistTPCvsGlobalMultAfterOutliers = (TH2F*)fListQAInfo->At(5);
352 
353  Float_t multTPC(0.); // tpc mult estimate
354  Float_t multGlobal(0.); // global multiplicity
355 
356  const Int_t nGoodTracks = fVevent->GetNumberOfTracks();
357  if(!fData2011) { // cut on outliers
358 
359  for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
360  AliAODTrack* AODtrack =dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
361  if (!AODtrack) continue;
362  if (!(AODtrack->TestFilterBit(1))) continue;
363  if ((AODtrack->Pt() < .2) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.2)) continue;
364  multTPC++;
365  }//track loop
366 
367  for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
368  AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
369  if (!AODtrack) continue;
370  if (!(AODtrack->TestFilterBit(16))) continue;
371  if ((AODtrack->Pt() < .2) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.1)) continue;
372  Double_t b[2] = {-99., -99.};
373  Double_t bCov[3] = {-99., -99., -99.};
374  AliAODTrack copy(*AODtrack);
375  if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
376  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
377  multGlobal++;
378  } //track loop
379 
380  HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);
381 
382  if(multTPC < (-40.3+1.22*multGlobal) || multTPC > (32.1+1.59*multGlobal)){ pass = kFALSE;}
383 
384  if(!pass) return;
385  HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);
386 
387  hNoEvents->Fill(4);
388 
389  }
390 
391 
392  if(fData2011) { // cut on outliers
393 
394  for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
395  AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
396  if (!AODtrack) continue;
397  if (!(AODtrack->TestFilterBit(1))) continue;
398  if ((AODtrack->Pt() < .2) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.2)) continue;
399  multTPC++;
400  }
401  for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
402  AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
403  if (!AODtrack) continue;
404  if (!(AODtrack->TestFilterBit(16))) continue;
405  if ((AODtrack->Pt() < .2) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.1)) continue;
406  Double_t b[2] = {-99., -99.};
407  Double_t bCov[3] = {-99., -99., -99.};
408  AliAODTrack copy(*AODtrack);
409  if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
410  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
411  multGlobal++;
412 
413  } //track loop
414 
415  HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);
416 
417  if(multTPC < (-36.73 + 1.48*multGlobal) || multTPC > (62.87 + 1.78*multGlobal)){pass = kFALSE;}
418 
419  if(!pass) return;
420  HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);
421 
422  hNoEvents->Fill(5);
423 
424  }
425 
426  for(Int_t itrack = 0; itrack < ntracks; itrack++){
427 
428  AliAODTrack *track=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(itrack));
429  if(!track) continue;
430 
431  Float_t dcaXY = track->DCA();
432  Float_t dcaZ = track->ZAtDCA();
433 
434  TH2F* HistDCAbefore =(TH2F*)fListQAInfo->At(6);
435  HistDCAbefore->Fill(dcaZ,dcaXY);
436 
437  Double_t p = -999, pT = -999, phi = -999, eta = -999, dEdx =-999;
438  Double_t length = -999., beta =-999, tofTime = -999., tof = -999.;
439  Double_t c = TMath::C()*1.E-9;// m/ns
440 
441  //cout<<"track->GetFilterMap()= "<<track->GetFilterMap()<<endl;
442  if(!track->TestFilterBit(fFilterBit)) continue;
443 
444  p=track->P();
445  pT=track->Pt();
446  phi=track->Phi();
447  eta=track->Eta();
448  dEdx=track->GetDetPid()->GetTPCsignal();
449 
450  Float_t probMis = fPIDResponse->GetTOFMismatchProbability(track);
451  if (probMis < 0.01) { //if u want to reduce mismatch using also TPC
452 
453  //if ( (track->IsOn(AliAODTrack::kTOFin)) && (track->IsOn(AliAODTrack::kTIME)) && (track->IsOn(AliAODTrack::kTOFout))) {
454  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
455 
456  tofTime = track->GetTOFsignal();//in ps
457  length = track->GetIntegratedLength();
458 
459  tof = tofTime*1E-3; // ns
460  //cout<<"tof = "<<tof<<endl;
461  if (tof <= 0) continue;
462  //cout<<"length = "<<length<<endl;
463  if (length <= 0) continue;
464 
465  length = length*0.01; // in meters
466  tof = tof*c;
467  beta = length/tof;
468 
469  TH2F *HistBetavsPTOFbeforePID = (TH2F*)fListQAInfo->At(7);
470  HistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
471  }//TOF signal
472 
473  TH2F *HistdEdxvsPTPCbeforePID = (TH2F*)fListQAInfo->At(8);
474  HistdEdxvsPTPCbeforePID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
475 
476  //QA plot
477  TH1F *HistPhiDistBefore = (TH1F*)fListQAInfo->At(9);
478  HistPhiDistBefore->Fill(phi);
479  //
480  TH1F *HistEtaDistBefore = (TH1F*)fListQAInfo->At(10);
481  HistEtaDistBefore->Fill(eta);
482 
483 
484  if(pT<0.1) continue;
485  if(TMath::Abs(eta)>0.8) continue;
486 
487  Int_t TPCNcls = track->GetTPCNcls();
488 
489  if(TPCNcls<70 || dEdx<10) continue;
490 
491  // fill QA histograms
492 
493  TH2F* HistDCAAfter =(TH2F*)fListQAInfo->At(11);
494  HistDCAAfter->Fill(dcaZ,dcaXY);
495 
496  TH1F *HistPhiDistAfter = (TH1F*)fListQAInfo->At(12);
497  HistPhiDistAfter->Fill(phi);
498 
499  TH1F *HistEtaDistAfter = (TH1F*)fListQAInfo->At(13);
500  HistEtaDistAfter->Fill(eta);
501 
502 
503  Bool_t pWithinRange = kFALSE;
504  Int_t p_bin = -999;
505  Double_t pBins[100];
506  for(int b=0;b<100;b++){pBins[b] = 0.1*b;}
507  for(int i=0;i<100;i++){
508  if(p>pBins[i] && p<(pBins[i]+0.1)){
509  pWithinRange = kTRUE;
510  p_bin = i;
511  }
512  }
513 
514  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
515 
516  Double_t nsigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
517  Double_t nsigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
518 
519  int i = ispecie - AliPID::kPion;
520 
521  if(fPIDcuts && pWithinRange){// for pions, kaons and protons only
522  if(ispecie==AliPID::kPion || ispecie==AliPID::kKaon || ispecie==AliPID::kProton){
523  int index = 100*i+p_bin;
524 
525  if(fPurityFunction[index]->Eval(nsigmaTOF,nsigmaTPC)>fPurityLevel){
526  if(TMath::Sqrt(TMath::Power(nsigmaTPC,2)+TMath::Power(nsigmaTOF,2))<3){
527 
528  TH3 *hist1 = (TH3*)fListQAtpctof->At(ispecie);
529  if (hist1){
530  hist1->Fill(nsigmaTPC,nsigmaTOF,p);}
531 
532  //--------THnsparse---------
533  fvalueAll[0] = p;
534  fvalueAll[1] = eta;
535  fvalueAll[2] = nsigmaTOF;
536  fvalueAll[3] = nsigmaTPC;
537  THnSparseD *fSparseAll = (THnSparseD*)fListQAtpctof->At(ispecie+(AliPID::kSPECIESC)*3.0);
538  fSparseAll->Fill(fvalueAll);
539 
540  }
541  }
542  if(p_bin>7 && fPurityFunction[index]->Eval(nsigmaTOF,nsigmaTPC)>fPurityLevel){//p_bin>7
543  if(TMath::Sqrt(TMath::Power(nsigmaTPC,2)+TMath::Power(nsigmaTOF,2))<3){
544  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
545  TH2F *HistBetavsPTOFafterPID = (TH2F*)fListQAInfo->At(14);
546  HistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
547  }
548  }
549  TH2F *HistdEdxvsPTPCafterPID = (TH2F*)fListQAInfo->At(15);
550  HistdEdxvsPTPCafterPID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
551  }
552 
553  if(p_bin<8 && TMath::Sqrt(TMath::Power(nsigmaTPC,2)+TMath::Power(nsigmaTOF,2))<3){//p_bin<8
554  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
555  TH2F *HistBetavsPTOFafterPID = (TH2F*)fListQAInfo->At(14);
556  HistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
557  }
558  TH2F *HistdEdxvsPTPCafterPID = (TH2F*)fListQAInfo->At(15);
559  HistdEdxvsPTPCafterPID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
560  }
561  //====================for low momentum PID based on purity by using only TPC information
562  Double_t LowPtPIDTPCnsigLow_Pion[6] = {-3,-3,-3,-3,-3,-3};
563  Double_t LowPtPIDTPCnsigLow_Kaon[6] = {-3,-2,0,-1.8,-1.2,-0.8}; //for 0.4<Pt<0.5 the purity is lower than 0.7
564  Double_t LowPtPIDTPCnsigHigh_Pion[6] ={2.4,3,3,3,2,1.4};
565  Double_t LowPtPIDTPCnsigHigh_Kaon[6] ={3,2.2,0,-0.2,1,1.8}; //for 0.4<Pt<0.5 the purity is lower than 0.7
566 
567 
568  if(p_bin>7 && fPurityFunction[index]->Eval(nsigmaTOF,nsigmaTPC)>fPurityLevel){//p_bin>7
569  if(TMath::Sqrt(TMath::Power(nsigmaTPC,2)+TMath::Power(nsigmaTOF,2))<3){
570  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
571  TH2F *HistBetavsPTOFafterPID = (TH2F*)fListQAInfo->At(29);
572  HistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
573  }
574  }
575  TH2F *HistdEdxvsPTPCafterPID = (TH2F*)fListQAInfo->At(30);
576  HistdEdxvsPTPCafterPID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
577  }
578  if(p_bin<8){
579  if((ispecie==AliPID::kPion && nsigmaTPC>LowPtPIDTPCnsigLow_Pion[p_bin-2] && nsigmaTPC<LowPtPIDTPCnsigHigh_Pion[p_bin-2]) || (ispecie==AliPID::kKaon && nsigmaTPC>LowPtPIDTPCnsigLow_Kaon[p_bin-2] && nsigmaTPC<LowPtPIDTPCnsigHigh_Kaon[p_bin-2]) || (ispecie==AliPID::kProton && nsigmaTPC>-3 && nsigmaTPC<3)){
580  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
581  TH2F *HistBetavsPTOFafterPID = (TH2F*)fListQAInfo->At(29);
582  HistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
583  }
584  TH2F *HistdEdxvsPTPCafterPID = (TH2F*)fListQAInfo->At(30);
585  HistdEdxvsPTPCafterPID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
586  }
587  }
588 
589 
590  TH2F *hTPCnSigmavsP = (TH2F*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
591  if (hTPCnSigmavsP){
592  hTPCnSigmavsP->Fill(track->P()*track->Charge(),nsigmaTPC);}
593 
594  TH2F *hTOFnSigmavsP = (TH2F*)fListQAtpctof->At(ispecie+(AliPID::kSPECIESC)*2.0);
595  if (hTOFnSigmavsP){
596  hTOFnSigmavsP->Fill(track->P()*track->Charge(),nsigmaTOF);}
597 
598  //=======================With TPC+TOF nsigma method Only!==============================
599  if(fPurityFunction[index]->Eval(nsigmaTOF,nsigmaTPC)>fPurityLevel){
600  if(TMath::Sqrt(TMath::Power(nsigmaTPC,2)+TMath::Power(nsigmaTOF,2))<3){
601 
602  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
603  TH2F *HistBetavsPTOFafterPIDTPCTOF = (TH2F*)fListQAInfo->At(16);
604  HistBetavsPTOFafterPIDTPCTOF ->Fill(track->P()*track->Charge(),beta);
605  if(ispecie==AliPID::kPion){
606  TH2F *HistPion_BetavsPTOFafterPIDTPCTOF = (TH2F*)fListQAInfo->At(20);
607  HistPion_BetavsPTOFafterPIDTPCTOF ->Fill(track->P()*track->Charge(),beta);
608  }
609  if(ispecie==AliPID::kKaon){
610  TH2F *HistKaon_BetavsPTOFafterPIDTPCTOF = (TH2F*)fListQAInfo->At(22);
611  HistKaon_BetavsPTOFafterPIDTPCTOF ->Fill(track->P()*track->Charge(),beta);
612  }
613  if(ispecie==AliPID::kProton){
614  TH2F *HistProton_BetavsPTOFafterPIDTPCTOF = (TH2F*)fListQAInfo->At(24);
615  HistProton_BetavsPTOFafterPIDTPCTOF ->Fill(track->P()*track->Charge(),beta);
616  }
617  }
618 
619  TH2F *HistdEdxvsPTPCafterPIDTPCTOF = (TH2F*)fListQAInfo->At(17);
620  HistdEdxvsPTPCafterPIDTPCTOF -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
621  if(ispecie==AliPID::kPion){
622  TH2F *HistPion_dEdxvsPTPCafterPIDTPCTOF = (TH2F*)fListQAInfo->At(21);
623  HistPion_dEdxvsPTPCafterPIDTPCTOF -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
624  TH1F *HistPionEta = (TH1F*)fListQAInfo->At(26);
625  HistPionEta->Fill(eta);
626  }
627  if(ispecie==AliPID::kKaon){
628  TH2F *HistKaon_dEdxvsPTPCafterPIDTPCTOF = (TH2F*)fListQAInfo->At(23);
629  HistKaon_dEdxvsPTPCafterPIDTPCTOF -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
630  TH1F *HistKaonEta = (TH1F*)fListQAInfo->At(27);
631  HistKaonEta->Fill(eta);
632  }
633  if(ispecie==AliPID::kProton){
634  TH2F *HistProton_dEdxvsPTPCafterPIDTPCTOF = (TH2F*)fListQAInfo->At(25);
635  HistProton_dEdxvsPTPCafterPIDTPCTOF -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
636  TH1F *HistProtonEta = (TH1F*)fListQAInfo->At(28);
637  HistProtonEta->Fill(eta);
638  }
639  }
640  }
641  //======================With TPC nsigma Only!
642  if(TMath::Sqrt(TMath::Power(nsigmaTPC,2)+TMath::Power(nsigmaTOF,2))<3){
643  if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
644  TH2F *HistBetavsPTOFafterPIDTPConly = (TH2F*)fListQAInfo->At(18);
645  HistBetavsPTOFafterPIDTPConly ->Fill(track->P()*track->Charge(),beta);
646  }
647  TH2F *HistdEdxvsPTPCafterPIDTPConly = (TH2F*)fListQAInfo->At(19);
648  HistdEdxvsPTPCafterPIDTPConly -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
649  }
650  //========================================================================================
651 
652 
653  }
654  }
655  if(!fPIDcuts){
656  TH3 *hist1 = (TH3*)fListQAtpctof->At(ispecie);
657  if (hist1){
658  hist1->Fill(nsigmaTPC,nsigmaTOF,p);}
659 
660  TH2F *hTPCnSigmavsP = (TH2F*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
661  if (hTPCnSigmavsP){
662  hTPCnSigmavsP->Fill(track->P(),nsigmaTPC);}
663 
664  TH2F *hTOFnSigmavsP = (TH2F*)fListQAtpctof->At(ispecie+(AliPID::kSPECIESC)*2.0);
665  if (hTOFnSigmavsP){
666  hTOFnSigmavsP->Fill(track->P(),nsigmaTOF);}
667 
668  //--------THnsparse---------
669  fvalueAll[0] = p;
670  fvalueAll[1] = eta;
671  fvalueAll[2] = nsigmaTOF;
672  fvalueAll[3] = nsigmaTPC;
673  THnSparseD *fSparseAll = (THnSparseD*)fListQAtpctof->At(ispecie+(AliPID::kSPECIESC)*3.0);
674  fSparseAll->Fill(fvalueAll);
675 
676 
677  }
678  }
679  }//probMis
680 
681  }//track loop
682 
683  TH2F *HistTPCvsGlobalMultAfter = (TH2F*) fListQAInfo->At(31);
684  HistTPCvsGlobalMultAfter->Fill(multGlobal,multTPC);
685 }
686 //_________________________________________
687 void AliAnalysisTaskPIDconfig::CheckCentrality(AliVEvent* event,Double_t centvalue, Bool_t &centralitypass)
688 {
689  // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
690  if (!fUseCentrality) AliFatal("No centrality method set! FATAL ERROR!");
691  //cout << "Centrality evaluated-------------------------: " << centvalue <<endl;
692  if ((centvalue >= fCentralityPercentileMin) && (centvalue < fCentralityPercentileMax))
693  {
694  TH1F *hCentralityPass = (TH1F*)fListQAInfo->At(2);
695  hCentralityPass->Fill(centvalue);
696  //cout << "--------------Fill pass-------------------------"<<endl;
697  centralitypass = kTRUE;
698  }
699 
700 }
701 //______________________________________________________________________________
703 {
704  fPurityLevel = PurityLevel;
705  fPurityFunctionsFile = TFile::Open(Form("alien:///alice/cern.ch/user/n/nmohamma/PurityFunctions_%i-%icent.root",fCentralityPercentileMin,fCentralityPercentileMax));
706 
707  if((!fPurityFunctionsFile) || (!fPurityFunctionsFile->IsOpen())) {
708  printf("The purity functions file does not exist");
709  return;
710  }
711 
712  fPurityFunctionsList=(TDirectory*)fPurityFunctionsFile->Get("Filterbit1");
713  if(!fPurityFunctionsList){printf("The purity functions list is empty"); return;}
714 
715  TString species[3] = {"pion","kaon","proton"};
716  TList *Species_functions[3];
717  Int_t ispecie = 0;
718  for(ispecie = 0; ispecie < 3; ispecie++) {
719  Species_functions[ispecie] = (TList*)fPurityFunctionsList->Get(species[ispecie]);
720  if(!Species_functions[ispecie]) {
721  cout<<"Purity functions for species: "<<species[ispecie]<<" not found!!!"<<endl;
722  return;
723  }
724  }
725 
726  for(int i=0;i<180;i++){
727  int ispecie = i/60;
728  int iPbin = i%60;
729  fPurityFunction[i] = (TF2*)Species_functions[ispecie]->FindObject(Form("PurityFunction_%d%d",iPbin,iPbin+1));
730  cout<<fPurityFunction[i]->GetName()<<" - Bin: "<<i<<endl;
731  if(!fPurityFunction[i]){printf("Purity function does not exist"); return;}
732  }
733 }
734 //______________________________________________________________________________
736 {
737  //
738  // Create the qa objects for TPC + TOF combination
739 
740 
741  //TPC and TOF signal vs. momentum
742  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
743  fhistNsigmaP = new TH3F(Form("NsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),Form("TPC n#sigma vs. TOF n#sigma %s vs. p ;TPC n#sigma;TOF n#sigma;p [GeV]",AliPID::ParticleName(ispecie)),200,-20,20,200,-20,20,100,0.1,10);
745  }
746 
747  //TPC signal vs. momentum
748  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
749  fhistTPCnSigmavsP = new TH2F(Form("NsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),Form("TPC n#sigma %s vs. p ;p [GeV];TPC n#sigma",AliPID::ParticleName(ispecie)),100,0,10,125,-5,20);
751  }
752 
753  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
754  fhistTOFnSigmavsP = new TH2F(Form("NsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),Form("TOF n#sigma %s vs. p ;p [GeV];TOF n#sigma",AliPID::ParticleName(ispecie)),100,0,10,150,-10,20);
756  }
757 
758  Int_t binsv1[4]={100,20,200,200}; //p,eta,tofsig,tpcsig
759  Double_t xminv1[4]={0,-1,-20,-20};
760  Double_t xmaxv1[4]={10,1,20,20};
761 
762  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
763  fSparseAll = new THnSparseD (Form("fSparse_%s",AliPID::ParticleName(ispecie)),Form("fSparse_%s",AliPID::ParticleName(ispecie)),4,binsv1,xminv1,xmaxv1);
764  fSparseAll->GetAxis(0)->SetTitle("p_{T}");
765  fSparseAll->GetAxis(1)->SetTitle("#eta");
766  fSparseAll->GetAxis(2)->SetTitle("TOFn#sigma");
767  fSparseAll->GetAxis(3)->SetTitle("TPCn#sigma");
768 
770  }
771 
772 
773 
774 }
775 //______________________________________________________________________________
777 {
778  //event and track info
779  fNoEvents = new TH1F("number of events","no. of events",5,0.5,5.5);
780  fNoEvents->GetXaxis()->SetBinLabel(1,"RawEvents");
781  fNoEvents->GetXaxis()->SetBinLabel(2,"AfterCentralityCut");
782  fNoEvents->GetXaxis()->SetBinLabel(3,"AfterVTXZCut");
783  fNoEvents->GetXaxis()->SetBinLabel(4,"AfterTPCGlobalOutliersCut2010");
784  fNoEvents->GetXaxis()->SetBinLabel(5,"AfterTPCGlobalOutliersCut2011");
785 
786  fListQAInfo->Add(fNoEvents);
787 
788  fhistCentralityPassBefore = new TH1F("fcentralityPassBefore","centralityPassBefore", 100,0,100);
790 
791  fhistCentralityPassAfter = new TH1F("fcentralityPassAfter","centralityPassAfter", 100,0,100);
793 
794  fpVtxZ = new TH1F("pVtxZ","pVtxZ",100,-20,20);
795  fListQAInfo->Add(fpVtxZ);
796 
797  fTPCvsGlobalMultBeforeOutliers = new TH2F("TPC vs. Global Multiplicity Before","TPC vs. Global Multiplicity Before",500,0,6000,500,0,6000);
799 
800  fTPCvsGlobalMultAfterOutliers = new TH2F("TPC vs. Global Multiplicity After outliers","TPC vs. Global Multiplicity After outliers",500,0,6000,500,0,6000);
802 
803  fhistDCABefore = new TH2F("DCA xy vs z (before)","DCA before",200,0,10,200,0,10);
805 
806  fHistBetavsPTOFbeforePID = new TH2F("momentum vs beta before PID","momentum vs beta before PID",1000,-10.,10.,1000,0,1.2);
808 
809  fHistdEdxvsPTPCbeforePID = new TH2F("momentum vs dEdx before PID","momentum vs dEdx before PID",1000,-10.,10.,1000,0,1000);
811 
812  fhistPhiDistBefore = new TH1F("Phi Distribution Before Cuts","Phi Distribution Before Cuts",200,0,6.4);
814 
815  fhistEtaDistBefore = new TH1F("Eta Distribution Before Cuts","Eta Distribution Before Cuts",100,-2,2);
817 
818  fhistDCAAfter = new TH2F("DCA xy vs z (after)","DCA after",200,0,10,200,0,10);
820 
821  fhistPhiDistAfter = new TH1F("Phi Distribution After Cuts","Phi Distribution After Cuts",200,0,6.4);
823 
824  fhistEtaDistAfter = new TH1F("Eta Distribution After Cuts","Eta Distribution After Cuts",200,-10,10);
826 
827  fHistBetavsPTOFafterPID = new TH2F("momentum vs beta after PID","momentum vs beta after PID",1000,-10.,10.,1000,0,1.2);
829 
830  fHistdEdxvsPTPCafterPID = new TH2F("momentum vs dEdx after PID","momentum vs dEdx after PID",1000,-10.,10.,1000,0,1000);
832 
833  fHistBetavsPTOFafterPIDTPCTOF = new TH2F("momentum vs beta after PID TPC+TOF","momentum vs beta after PID TPC+TOF",1000,-10.,10.,1000,0,1.2);
835 
836  fHistdEdxvsPTPCafterPIDTPCTOF = new TH2F("momentum vs dEdx after PID TPC+TOF","momentum vs dEdx after PID TPC+TOF",1000,-10.,10.,1000,0,1000);
838 
839  fHistBetavsPTOFafterPIDTPConly = new TH2F("momentum vs beta after PID TPC only","momentum vs beta after PID TPC only",1000,-10.,10.,1000,0,1.2);
841 
842  fHistdEdxvsPTPCafterPIDTPConly = new TH2F("momentum vs dEdx after PID TPC only","momentum vs dEdx after PID TPC only",1000,-10.,10.,1000,0,1000);
844 
845  fHistPion_BetavsPTOFafterPIDTPCTOF = new TH2F("Pion momentum vs beta after PID TPC+TOF","Pion momentum vs beta after PID TPC+TOF",1000,-10.,10.,1000,0,1.2);
847 
848  fHistPion_dEdxvsPTPCafterPIDTPCTOF = new TH2F("Pion momentum vs dEdx after PID TPC+TOF","Pion momentum vs dEdx after PID TPC+TOF",1000,-10.,10.,1000,0,1000);
850 
851  fHistKaon_BetavsPTOFafterPIDTPCTOF = new TH2F("Kaon momentum vs beta after PID TPC+TOF","Kaon momentum vs beta after PID TPC+TOF",1000,-10.,10.,1000,0,1.2);
853 
854  fHistKaon_dEdxvsPTPCafterPIDTPCTOF = new TH2F("Kaon momentum vs dEdx after PID TPC+TOF","Kaon momentum vs dEdx after PID TPC+TOF",1000,-10.,10.,1000,0,1000);
856 
857  fHistProton_BetavsPTOFafterPIDTPCTOF = new TH2F("Proton momentum vs beta after PID TPC+TOF","Proton momentum vs beta after PID TPC+TOF",1000,-10.,10.,1000,0,1.2);
859 
860  fHistProton_dEdxvsPTPCafterPIDTPCTOF = new TH2F("Proton momentum vs dEdx after PID TPC+TOF","Proton momentum vs dEdx after PID TPC+TOF",1000,-10.,10.,1000,0,1000);
862 
863  fhistPionEtaDistAfter = new TH1F("Pion Eta Distribution After PID Cuts","Pion Eta Distribution After PID Cuts",100,-2,2);
865 
866  fhistKaonEtaDistAfter = new TH1F("Kaon Eta Distribution After PID Cuts","Kaon Eta Distribution After PID Cuts",100,-2,2);
868 
869  fhistProtonEtaDistAfter = new TH1F("Proton Eta Distribution After PID Cuts","Proton Eta Distribution PID After Cuts",100,-2,2);
871 
872  fHistBetavsPTOFafterPID_2 = new TH2F("momentum vs beta after PID (PID in low Pt TPC only with Purity>0.7)","momentum vs beta after PID (PID in low Pt TPC only with Purity>0.7)",1000,-10.,10.,1000,0,1.2);
874 
875  fHistdEdxvsPTPCafterPID_2 = new TH2F("momentum vs dEdx after PID (PID in low Pt TPC only with Purity>0.7)","momentum vs dEdx after PID (PID in low Pt TPC only with Purity>0.7)",1000,-10.,10.,1000,0,1000);
877 
878  fTPCvsGlobalMultAfter = new TH2F("TPC vs. Global Multiplicity After","TPC vs. Global Multiplicity After",500,0,6000,500,0,6000);
880 
881 }
882 
883 
TList * fListQAInfo
List with combined PID from TPC + TOF.
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
TH1F * fhistCentralityPassAfter
cen histo before
Definition: External.C:244
TH1F * fNoEvents
cen histo after
AliPIDResponse * fPIDResponse
aod
TH2F * fHistdEdxvsPTPCafterPIDTPCTOF
another hist
TH1F * fhistProtonEtaDistAfter
another hist
ClassImp(AliAnalysisTaskPIDconfig) AliAnalysisTaskPIDconfig
centrality
TH2F * fHistdEdxvsPTPCafterPID_2
another hist
TH1F * fhistPhiDistAfter
another hist
TCanvas * c
Definition: TestFitELoss.C:172
TH2F * fHistPion_dEdxvsPTPCafterPIDTPCTOF
another hist
TH1F * fhistPionEtaDistAfter
another hist
TH2F * fHistKaon_dEdxvsPTPCafterPIDTPCTOF
another hist
TH2F * fHistKaon_BetavsPTOFafterPIDTPCTOF
another hist
TH2F * fTPCvsGlobalMultAfterOutliers
another hist
TH1F * fhistEtaDistBefore
another hist
TH2F * fhistTPCnSigmavsP
another hist
TH1F * fhistCentralityPassBefore
list q ainfo
int Int_t
Definition: External.C:63
TH2F * fHistBetavsPTOFafterPIDTPCTOF
another hist
float Float_t
Definition: External.C:68
Bool_t fData2011
All species info.
TH2F * fHistdEdxvsPTPCafterPID
another hist
TH2F * fHistProton_dEdxvsPTPCafterPIDTPCTOF
another hist
TDirectory * fPurityFunctionsList
purity functions file
void CheckCentrality(AliVEvent *event, Double_t centrality, Bool_t &centralitypass)
TH2F * fHistBetavsPTOFafterPID
another hist
TObject * FindObject(int bin, const char *nameH, const TList *lst, Bool_t normPerEvent=kTRUE)
TH1F * fhistEtaDistAfter
another hist
TH2F * fHistBetavsPTOFafterPIDTPConly
another hist
TH2F * fHistBetavsPTOFafterPID_2
another hist
TH2F * fTPCvsGlobalMultAfter
another hist
TH2F * fHistdEdxvsPTPCbeforePID
another hist
TH2F * fHistProton_BetavsPTOFafterPIDTPCTOF
another hist
TH2F * fHistdEdxvsPTPCafterPIDTPConly
another hist
THnSparseD * fSparseAll
another hist
virtual void UserExec(Option_t *)
const char Option_t
Definition: External.C:48
TH2F * fHistBetavsPTOFbeforePID
another hist
TH1F * fhistKaonEtaDistAfter
another hist
void SetPIDPurityFunctions(Float_t PurityLevel)
bool Bool_t
Definition: External.C:53
TList * fListQAtpctof
List of all lists.
TList * fListQA
purity functions list
TH2F * fhistTOFnSigmavsP
another hist
TF2 * fPurityFunction[180]
All species info.
TH2F * fhistDCAAfter
dca after hist
TH2F * fTPCvsGlobalMultBeforeOutliers
another hist
TH2F * fHistPion_BetavsPTOFafterPIDTPCTOF
another hist
TH1F * fhistPhiDistBefore
another hist