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