fix missing initialisation in named constructor
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskPIDconfig.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: 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
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"
50 #include "AliPIDCombined.h"
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
74 using std::cout;
75 using std::endl;
76
77
78 ClassImp(AliAnalysisTaskPIDconfig)
79 //ClassImp()
80 //___________________________________________________________________
81 AliAnalysisTaskPIDconfig::AliAnalysisTaskPIDconfig():
82 AliAnalysisTaskSE(),
83 fVevent(0),
84 fESD(0),
85 fAOD(0),
86 fPIDResponse(0),
87 fTriggerSelection(0),
88 fCentralityPercentileMin(0.),
89 fCentralityPercentileMax(5.),
90 fFilterBit(128),
91 fDCAxyCut(-1),
92 fDCAzCut(-1),
93 fData2011(kFALSE),
94 fTriggerMB(kTRUE),
95 fTriggerCentral(kFALSE),
96 fUseCentrality(kTRUE),
97 fCutTPCmultiplicityOutliersAOD(kFALSE),
98 fPIDcuts(kFALSE),
99 fCentralityEstimator("V0M"),
100 fContoursFile(0),
101 fCutContourList(0),
102 fListQA(0x0),
103 fListQAtpctof(0x0),
104 fListQAInfo(0x0),
105 fhistCentralityPass(0),
106 fNoEvents(0),
107 fpVtxZ(0),
108 fhistDCABefore(0),
109 fhistDCAAfter(0),
110 fhistPhiDistBefore(0),
111 fhistPhiDistAfter(0),
112 fhistEtaDistBefore(0),
113 fhistEtaDistAfter(0),
114 fTPCvsGlobalMultBeforeOutliers(0),
115 fTPCvsGlobalMultAfterOutliers(0),
116 fTPCvsGlobalMultAfter(0),
117 fHistBetavsPTOFbeforePID(0),
118 fHistdEdxvsPTPCbeforePID(0),
119 fhistNsigmaP(0),
120 fhistNsigmaPTPCTOF(0),
121 fhistTPCnSigmavsP(0),
122 fHistBetavsPTOFafterPID(0),
123 fHistdEdxvsPTPCafterPID(0),
124 fHistBetavsPTOFafterPIDTPCTOF(0),
125 fHistdEdxvsPTPCafterPIDTPCTOF(0),
126 fHistBetavsPTOFafterPIDTPConly(0),
127 fHistdEdxvsPTPCafterPIDTPConly(0)
128 //fSparseSpecies(0),
129 //fvalueSpecies(0)
130 {
131     for(int i=0;i<150;i++){
132         fCutContour[i]= NULL;
133         fCutGraph[i]=NULL;
134     }
135 }
136
137
138 //___________________________________________________________________
139
140 AliAnalysisTaskPIDconfig::AliAnalysisTaskPIDconfig(const char *name):
141 AliAnalysisTaskSE(name),
142 fVevent(0),
143 fESD(0),
144 fAOD(0),
145 fPIDResponse(0),
146 fTriggerSelection(0),
147 fCentralityPercentileMin(0.),
148 fCentralityPercentileMax(5.),
149 fFilterBit(128),
150 fDCAxyCut(-1),
151 fDCAzCut(-1),
152 fData2011(kFALSE),
153 fTriggerMB(kTRUE),
154 fTriggerCentral(kFALSE),
155 fUseCentrality(kTRUE),
156 fCutTPCmultiplicityOutliersAOD(kFALSE),
157 fPIDcuts(kFALSE),
158 fCentralityEstimator("V0M"),
159 fContoursFile(0),
160 fCutContourList(0),
161 fListQA(0x0),
162 fListQAtpctof(0x0),
163 fListQAInfo(0x0),
164 fhistCentralityPass(0),
165 fNoEvents(0),
166 fpVtxZ(0),
167 fhistDCABefore(0),
168 fhistDCAAfter(0),
169 fhistPhiDistBefore(0),
170 fhistPhiDistAfter(0),
171 fhistEtaDistBefore(0),
172 fhistEtaDistAfter(0),
173 fTPCvsGlobalMultBeforeOutliers(0),
174 fTPCvsGlobalMultAfterOutliers(0),
175 fTPCvsGlobalMultAfter(0),
176 fHistBetavsPTOFbeforePID(0),
177 fHistdEdxvsPTPCbeforePID(0),
178 fhistNsigmaP(0),
179 fhistNsigmaPTPCTOF(0),
180 fhistTPCnSigmavsP(0),
181 fHistBetavsPTOFafterPID(0),
182 fHistdEdxvsPTPCafterPID(0),
183 fHistBetavsPTOFafterPIDTPCTOF(0),
184 fHistdEdxvsPTPCafterPIDTPCTOF(0),
185 fHistBetavsPTOFafterPIDTPConly(0),
186 fHistdEdxvsPTPCafterPIDTPConly(0)
187 //fSparseSpecies(0),
188 //fvalueSpecies(0)
189 {
190     //Default Constructor
191     //fCutContour[150]=NULL;
192     //fCutGraph[150]=NULL;
193     DefineInput(0,TChain::Class());
194     DefineOutput(1,TList::Class());
195 }
196
197 //_____________________________________________________________________
198 AliAnalysisTaskPIDconfig::~AliAnalysisTaskPIDconfig()
199 {
200     //Destructor
201     
202     fContoursFile->Close();
203     for(int i=0;i<150;i++){
204         delete fCutContour[i];
205         delete fCutGraph[i];
206     }
207     
208     
209     //   delete fPID;
210     //  delete fPIDqa;
211     //   delete fTrackCuts;
212     // delete fSparseSpecies;
213     //delete []fvalueSpecies;
214     
215     
216 }
217 //______________________________________________________________________
218 void AliAnalysisTaskPIDconfig::UserCreateOutputObjects()
219 {
220     //
221     // Create the output QA objects
222     //
223     
224     AliLog::SetClassDebugLevel("AliAnalysisTaskPIDconfig",10);
225     
226     //input hander
227     AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
228     AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
229     if (!inputHandler) AliFatal("Input handler needed");
230     
231     //pid response object
232     fPIDResponse=inputHandler->GetPIDResponse();
233     if (!fPIDResponse) AliError("PIDResponse object was not created");
234     
235     if(fPIDcuts){ GetPIDContours(); cout<<"********** PID cut contours retrieved **********"<<endl;}
236     //
237     fListQA=new TList;
238     fListQA->SetOwner();
239     
240     fListQAtpctof=new TList;
241     fListQAtpctof->SetOwner();
242     fListQAtpctof->SetName("PID_TPC_TOF");
243     
244     fListQAInfo=new TList;
245     fListQAInfo->SetOwner();
246     fListQAInfo->SetName("Event_Track_Info");
247     
248     fListQA->Add(fListQAtpctof);
249     fListQA->Add(fListQAInfo);
250     
251     SetupTPCTOFqa();
252     SetupEventInfo();
253     
254     PostData(1,fListQA);
255 }
256 //______________________________________________________________________
257 void AliAnalysisTaskPIDconfig::UserExec(Option_t*){
258     //Main loop
259     //Called for each event
260     
261     // create pointer to event
262     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
263     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
264     
265     Int_t ntracks=fAOD->GetNumberOfTracks();
266     
267     
268     if(!(fESD || fAOD)){
269         printf("ERROR: fESD & fAOD not available\n");
270         return;
271     }
272     fVevent = dynamic_cast<AliVEvent*>(InputEvent());
273     if (!fVevent) {
274         printf("ERROR: fVevent not available\n");
275         return;
276     }
277     
278     Bool_t pass = kFALSE;
279     CheckCentrality(fVevent,pass);
280     
281     if(!pass){ return;}
282     
283     const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
284     
285     Double_t pVtxZ = -999;
286     pVtxZ = pVtx->GetZ();
287     
288     if(TMath::Abs(pVtxZ)>10) return;
289     
290     TH1F *hNoEvents = (TH1F*)fListQAInfo->At(1);
291     TH1F *histpVtxZ = (TH1F*)fListQAInfo->At(2);
292     
293     if(hNoEvents) hNoEvents->Fill(0);
294     if(histpVtxZ) histpVtxZ->Fill(pVtxZ);
295     
296     if(ntracks<2) return;
297     
298     // if(!pass) return;
299     
300     TH2F *HistTPCvsGlobalMultBeforeOutliers = (TH2F*)fListQAInfo->At(3);
301     
302     TH2F *HistTPCvsGlobalMultAfterOutliers = (TH2F*)fListQAInfo->At(4);
303     
304     Float_t multTPC(0.); // tpc mult estimate
305     Float_t multGlobal(0.); // global multiplicity
306     
307     const Int_t nGoodTracks = fVevent->GetNumberOfTracks();
308     if(!fData2011) { // cut on outliers
309         
310         for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
311             AliAODTrack* AODtrack =dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
312             if (!AODtrack) continue;
313             if (!(AODtrack->TestFilterBit(1))) continue;
314             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;
315             multTPC++;
316         }//track loop
317         
318         for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
319             AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
320             if (!AODtrack) continue;
321             if (!(AODtrack->TestFilterBit(16))) continue;
322             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;
323             Double_t b[2] = {-99., -99.};
324             Double_t bCov[3] = {-99., -99., -99.};
325             AliAODTrack copy(*AODtrack);
326             if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
327             if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
328             multGlobal++;
329         } //track loop
330         
331         HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);
332         
333         if(multTPC < (-40.3+1.22*multGlobal) || multTPC > (32.1+1.59*multGlobal)){ pass = kFALSE;}
334         
335         if(!pass) return;
336         HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);
337         
338     }
339     
340     
341     if(fData2011) { // cut on outliers
342         //Float_t multTPC(0.); // tpc mult estimate
343         //Float_t multGlob(0.); // global multiplicity
344         for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
345             AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
346             if (!AODtrack) continue;
347             if (!(AODtrack->TestFilterBit(1))) continue;
348             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;
349             multTPC++;
350         }
351         for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
352             AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
353             if (!AODtrack) continue;
354             if (!(AODtrack->TestFilterBit(16))) continue;
355             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;
356             Double_t b[2] = {-99., -99.};
357             Double_t bCov[3] = {-99., -99., -99.};
358             AliAODTrack copy(*AODtrack);
359             if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
360             if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
361             multGlobal++;
362             
363         } //track loop
364         
365         HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);
366         
367         if(multTPC < (-36.73 + 1.48*multGlobal) || multTPC > (62.87 + 1.78*multGlobal)){pass = kFALSE;}
368         
369         if(!pass) return;
370         HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);
371         
372     }
373     
374     for(Int_t itrack = 0; itrack < ntracks; itrack++){
375         
376         AliAODTrack *track=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(itrack));
377         if(!track) continue;
378         
379         Float_t dcaXY = track->DCA();
380         Float_t dcaZ  = track->ZAtDCA();
381         
382         TH2F* HistDCAbefore =(TH2F*)fListQAInfo->At(5);
383         HistDCAbefore->Fill(dcaZ,dcaXY);
384         
385         Double_t p = -999, pT = -999, phi = -999, eta = -999, dEdx =-999;
386         Double_t length = -999., beta =-999, tofTime = -999., tof = -999.;
387         Double_t c = TMath::C()*1.E-9;// m/ns
388         
389         //cout<<"track->GetFilterMap()= "<<track->GetFilterMap()<<endl;
390         if(!track->TestFilterBit(fFilterBit)) continue;
391         
392         p=track->P();
393         pT=track->Pt();
394         phi=track->Phi();
395         eta=track->Eta();
396         dEdx=track->GetDetPid()->GetTPCsignal();
397         
398         Float_t probMis = fPIDResponse->GetTOFMismatchProbability(track);
399         if (probMis < 0.01) { //if u want to reduce mismatch using also TPC
400             
401             //if ( (track->IsOn(AliAODTrack::kTOFin)) && (track->IsOn(AliAODTrack::kTIME)) &&  (track->IsOn(AliAODTrack::kTOFout))) {
402             if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
403                 
404                 tofTime = track->GetTOFsignal();//in ps
405                 length = track->GetIntegratedLength();
406                 
407                 tof = tofTime*1E-3; // ns
408                 //cout<<"tof = "<<tof<<endl;
409                 if (tof <= 0) continue;
410                 //cout<<"length = "<<length<<endl;
411                 if (length <= 0) continue;
412                 
413                 length = length*0.01; // in meters
414                 tof = tof*c;
415                 beta = length/tof;
416                 
417                 TH2F *HistBetavsPTOFbeforePID = (TH2F*)fListQAInfo->At(6);
418                 HistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
419             }//TOF signal
420             
421             TH2F *HistdEdxvsPTPCbeforePID = (TH2F*)fListQAInfo->At(7);
422             HistdEdxvsPTPCbeforePID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
423         
424             //QA plot
425             TH1F *HistPhiDistBefore = (TH1F*)fListQAInfo->At(8);
426             HistPhiDistBefore->Fill(phi);
427             //
428             TH1F *HistEtaDistBefore = (TH1F*)fListQAInfo->At(9);
429             HistEtaDistBefore->Fill(eta);
430         
431         
432             if(pT<0.1) continue;
433             if(TMath::Abs(eta)>0.8) continue;
434         
435             Int_t TPCNcls = track->GetTPCNcls();
436         
437             if(TPCNcls<70 || dEdx<10) continue;
438         
439             // fill QA histograms
440         
441             TH2F* HistDCAAfter =(TH2F*)fListQAInfo->At(10);
442             HistDCAAfter->Fill(dcaZ,dcaXY);
443         
444             TH1F *HistPhiDistAfter = (TH1F*)fListQAInfo->At(11);
445             HistPhiDistAfter->Fill(phi);
446         
447             TH1F *HistEtaDistAfter = (TH1F*)fListQAInfo->At(12);
448             HistEtaDistAfter->Fill(eta);
449         
450         
451             Bool_t pWithinRange = kFALSE;
452             Int_t p_bin = -999;
453             Double_t pBins[50];
454             for(int b=0;b<50;b++){pBins[b] = 0.1*b;}
455             for(int i=0;i<50;i++){
456                 if(p>pBins[i] && p<(pBins[i]+0.1)){
457                     pWithinRange = kTRUE;
458                     p_bin = i;
459                 }
460             }
461         
462             for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
463
464                 Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
465                 Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
466             
467                 int i = ispecie - AliPID::kPion;
468             
469                 if(fPIDcuts && pWithinRange){// for pions, kaons and protons only
470                     if(ispecie==AliPID::kPion || ispecie==AliPID::kKaon || ispecie==AliPID::kProton){
471                         int index = 50*i+p_bin;
472                        
473                         if(p_bin>7 && fCutContour[index]->IsInside(nSigmaTOF,nSigmaTPC)){//p_bin>7
474                             TH3 *hist1 = (TH3*)fListQAtpctof->At(ispecie);
475                             if (hist1){
476                                 hist1->Fill(nSigmaTPC,nSigmaTOF,p);}
477                             
478                             if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
479                                 TH2F *HistBetavsPTOFafterPID = (TH2F*)fListQAInfo->At(13);
480                                 HistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
481                             }
482                             TH2F *HistdEdxvsPTPCafterPID = (TH2F*)fListQAInfo->At(14);
483                             HistdEdxvsPTPCafterPID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
484                         }
485                         
486                         if(p_bin<8 && nSigmaTPC<3 && nSigmaTPC>-3){//p_bin<8
487                             if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
488                                 TH2F *HistBetavsPTOFafterPID = (TH2F*)fListQAInfo->At(13);
489                                 HistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
490                             }
491                             TH2F *HistdEdxvsPTPCafterPID = (TH2F*)fListQAInfo->At(14);
492                             HistdEdxvsPTPCafterPID -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
493                         }
494                         
495                         TH2F *hTPCnSigmavsP = (TH2F*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
496                         if (hTPCnSigmavsP){
497                             hTPCnSigmavsP->Fill(track->P()*track->Charge(),nSigmaTPC);}
498                         
499                         //=======================With TPC+TOF nsigma method Only!==============================
500                         if(fCutContour[index]->IsInside(nSigmaTOF,nSigmaTPC)){
501                             TH3 *hist2 = (TH3*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC);
502                             if (hist2){
503                                 hist2->Fill(nSigmaTPC,nSigmaTOF,p);}
504                             
505                             if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
506                                 TH2F *HistBetavsPTOFafterPIDTPCTOF = (TH2F*)fListQAInfo->At(15);
507                                 HistBetavsPTOFafterPIDTPCTOF ->Fill(track->P()*track->Charge(),beta);
508                             }
509                             TH2F *HistdEdxvsPTPCafterPIDTPCTOF = (TH2F*)fListQAInfo->At(16);
510                             HistdEdxvsPTPCafterPIDTPCTOF -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
511                         }
512                         //======================With TPC nsigma Only!
513                         if(nSigmaTPC<3 && nSigmaTPC>-3){
514                             if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
515                                 TH2F *HistBetavsPTOFafterPIDTPConly = (TH2F*)fListQAInfo->At(17);
516                                 HistBetavsPTOFafterPIDTPConly ->Fill(track->P()*track->Charge(),beta);
517                             }
518                             TH2F *HistdEdxvsPTPCafterPIDTPConly = (TH2F*)fListQAInfo->At(18);
519                             HistdEdxvsPTPCafterPIDTPConly -> Fill(track->P()*track->Charge(),dEdx); //TPC signal
520                         }
521                         //========================================================================================
522                         
523                         
524                         
525                     }
526                 }
527                 if(!fPIDcuts){
528                     TH3 *hist1 = (TH3*)fListQAtpctof->At(ispecie);
529                     if (hist1){
530                         hist1->Fill(nSigmaTPC,nSigmaTOF,p);}
531                     
532                     TH2F *hTPCnSigmavsP = (TH2F*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
533                     if (hTPCnSigmavsP){
534                         hTPCnSigmavsP->Fill(nSigmaTPC,p);}
535                     
536                 }
537             }
538         }//probMis
539         
540     }//track loop
541     
542     TH2F *HistTPCvsGlobalMultAfter = (TH2F*) fListQAInfo->At(19);
543     HistTPCvsGlobalMultAfter->Fill(multGlobal,multTPC);
544     
545 }
546 //_________________________________________
547 void AliAnalysisTaskPIDconfig::CheckCentrality(AliVEvent* event, Bool_t &centralitypass)
548 {
549     // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
550     if (!fUseCentrality) AliFatal("No centrality method set! FATAL ERROR!");
551     Double_t centvalue = event->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
552     //cout << "Centrality evaluated-------------------------: " << centvalue <<endl;
553     if ((centvalue >= fCentralityPercentileMin) && (centvalue < fCentralityPercentileMax))
554     {
555         TH1F *hCentralityPass = (TH1F*)fListQAInfo->At(0);
556         hCentralityPass->Fill(centvalue);
557         //cout << "--------------Fill pass-------------------------"<<endl;
558         centralitypass = kTRUE;
559     }
560     
561 }
562 //______________________________________________________________________________
563 void AliAnalysisTaskPIDconfig::GetPIDContours()
564 {
565     fContoursFile = new TFile(Form("$ALICE_ROOT/PWGCF/FLOW/database/PIDCutContours_%i-%i.root",fCentralityPercentileMin,fCentralityPercentileMax));
566     
567     fCutContourList=(TDirectory*)fContoursFile->Get("Filterbit1");
568     if(!fCutContourList){printf("The contour file is empty"); return;}
569
570     Double_t pBinning[50];
571     for(int b=0;b<50;b++){pBinning[b]=b;}
572     TString species[3] = {"pion","kaon","proton"};
573     
574     for(int i=0;i<150;i++){
575         int ispecie = i/50;
576         int iPbin = i%50;
577         TList *Species_contours = (TList*)fCutContourList->Get(species[ispecie]);
578         //if(Species_contours){cout<<"Species_contours exists"<<endl;}
579         
580         TString Graph_Name = "contourlines_";
581         Graph_Name += species[ispecie];
582         Graph_Name += Form("%.f%.f-%i%icent",pBinning[iPbin],pBinning[iPbin]+1,fCentralityPercentileMin,fCentralityPercentileMax);
583         //cout<<Graph_Name<<endl;
584         fCutGraph[i] = (TGraph*)Species_contours->FindObject(Graph_Name);
585         
586         if(!fCutGraph[i]){cout<<"Contour Graph does not exist"<<endl; continue;}
587         
588         fCutContour[i] = new TCutG(Graph_Name.Data(),fCutGraph[i]->GetN(),fCutGraph[i]->GetX(),fCutGraph[i]->GetY());
589         
590     }
591     
592 }
593 //______________________________________________________________________________
594 void AliAnalysisTaskPIDconfig::SetupTPCTOFqa()
595 {
596     //
597     // Create the qa objects for TPC + TOF combination
598     
599     
600     //TPC and TOF signal vs. momentum
601     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
602         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);
603         fListQAtpctof->Add(fhistNsigmaP);
604     }
605     
606     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
607         fhistNsigmaPTPCTOF = new TH3F(Form("NsigmaP_TPC_TOF_%s_TPC+TOFMethod",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);
608         fListQAtpctof->Add(fhistNsigmaPTPCTOF);
609     }
610
611     //TPC signal vs. momentum
612     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
613         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);
614         fListQAtpctof->Add(fhistTPCnSigmavsP);
615     }
616     
617 }
618 //______________________________________________________________________________
619 void AliAnalysisTaskPIDconfig::SetupEventInfo()
620 {
621     //event and track info
622     
623     fhistCentralityPass = new TH1F("fcentralityPass","centralityPass", 100,0,100);
624     fListQAInfo->Add(fhistCentralityPass);
625     
626     fNoEvents = new TH1F("number of events","no. of events",1,0,1);
627     fListQAInfo->Add(fNoEvents);
628     
629     fpVtxZ = new TH1F("pVtxZ","pVtxZ",100,-20,20);
630     fListQAInfo->Add(fpVtxZ);
631     
632     fTPCvsGlobalMultBeforeOutliers = new TH2F("TPC vs. Global Multiplicity Before","TPC vs. Global Multiplicity Before",500,0,6000,500,0,6000);
633     fListQAInfo->Add(fTPCvsGlobalMultBeforeOutliers);
634     
635     fTPCvsGlobalMultAfterOutliers = new TH2F("TPC vs. Global Multiplicity After outliers","TPC vs. Global Multiplicity After outliers",500,0,6000,500,0,6000);
636     fListQAInfo->Add(fTPCvsGlobalMultAfterOutliers);
637     
638     fhistDCABefore = new TH2F("DCA xy vs z (before)","DCA before",200,0,10,200,0,10);
639     fListQAInfo->Add(fhistDCABefore);
640     
641     fHistBetavsPTOFbeforePID = new TH2F("momentum vs beta before PID","momentum vs beta before PID",1000,-10.,10.,1000,0,1.2);
642     fListQAInfo->Add(fHistBetavsPTOFbeforePID);
643     
644     fHistdEdxvsPTPCbeforePID = new TH2F("momentum vs dEdx before PID","momentum vs dEdx before PID",1000,-10.,10.,1000,0,1000);
645     fListQAInfo->Add(fHistdEdxvsPTPCbeforePID);
646     
647     fhistPhiDistBefore = new TH1F("Phi Distribution Before Cuts","Phi Distribution Before Cuts",200,0,6.4);
648     fListQAInfo->Add(fhistPhiDistBefore);
649     
650     fhistEtaDistBefore = new TH1F("Eta Distribution Before Cuts","Eta Distribution Before Cuts",200,-10,10);
651     fListQAInfo->Add(fhistEtaDistBefore);
652     
653     fhistDCAAfter = new TH2F("DCA xy vs z (after)","DCA after",200,0,10,200,0,10);
654     fListQAInfo->Add(fhistDCAAfter);
655     
656     fhistPhiDistAfter = new TH1F("Phi Distribution After Cuts","Phi Distribution After Cuts",200,0,6.4);
657     fListQAInfo->Add(fhistPhiDistAfter);
658     
659     fhistEtaDistAfter = new TH1F("Eta Distribution After Cuts","Eta Distribution After Cuts",200,-10,10);
660     fListQAInfo->Add(fhistEtaDistAfter);
661     
662     fHistBetavsPTOFafterPID = new TH2F("momentum vs beta after PID","momentum vs beta after PID",1000,-10.,10.,1000,0,1.2);
663     fListQAInfo->Add(fHistBetavsPTOFafterPID);
664     
665     fHistdEdxvsPTPCafterPID = new TH2F("momentum vs dEdx after PID","momentum vs dEdx after PID",1000,-10.,10.,1000,0,1000);
666     fListQAInfo->Add(fHistdEdxvsPTPCafterPID);
667     
668     fHistBetavsPTOFafterPIDTPCTOF = new TH2F("momentum vs beta after PID TPC+TOF","momentum vs beta after PID",1000,-10.,10.,1000,0,1.2);
669     fListQAInfo->Add(fHistBetavsPTOFafterPID);
670     
671     fHistdEdxvsPTPCafterPIDTPCTOF = new TH2F("momentum vs dEdx after PID TPC+TOF","momentum vs dEdx after PID",1000,-10.,10.,1000,0,1000);
672     fListQAInfo->Add(fHistdEdxvsPTPCafterPID);
673     
674     fHistBetavsPTOFafterPIDTPConly = new TH2F("momentum vs beta after PID TPC only","momentum vs beta after PID",1000,-10.,10.,1000,0,1.2);
675     fListQAInfo->Add(fHistBetavsPTOFafterPID);
676     
677     fHistdEdxvsPTPCafterPIDTPConly = new TH2F("momentum vs dEdx after PID TPC only","momentum vs dEdx after PID",1000,-10.,10.,1000,0,1000);
678     fListQAInfo->Add(fHistdEdxvsPTPCafterPID);
679     
680     fTPCvsGlobalMultAfter = new TH2F("TPC vs. Global Multiplicity After","TPC vs. Global Multiplicity After",500,0,6000,500,0,6000);
681     fListQAInfo->Add(fTPCvsGlobalMultAfter);
682     
683 }
684
685