AliPhysics  d565ceb (d565ceb)
AliAnalysisTaskCMEV0PID.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: AliAnalysisTaskCMEV0PID.cxx Rihan Haque 14/02/2018 $ */
17 
18 //-- general include---
19 #include "TChain.h"
20 #include "TTree.h"
21 #include "TGrid.h"
22 #include "TROOT.h"
23 #include "TObjArray.h"
24 #include "TMatrixDSym.h"
25 
26 #include "TMath.h"
27 #include "stdio.h"
28 #include "Riostream.h"
29 
30 //---- manager and handler---
31 #include "AliAnalysisManager.h"
32 #include "AliInputEventHandler.h"
33 
34 //---V0 and ZDC info---
35 #include "AliAODZDC.h"
36 #include "AliAODVZERO.h"
37 #include "AliAODVertex.h"
38 
39 //---AOD,ESD event--
40 #include "AliESDEvent.h"
41 #include "AliAODHeader.h"
42 #include "AliAODEvent.h"
43 #include "AliAODTrack.h"
44 
45 //----for PID-----
46 #include "AliPIDResponse.h"
47 #include "AliPIDCombined.h"
48 
49 //----- Vevent and tracks
50 #include "AliVEventHandler.h"
51 #include "AliVEvent.h"
52 #include "AliVTrack.h"
53 #include "AliVParticle.h"
54 
55 //----- must include-------
56 #include "AliMultSelection.h"
57 #include "AliAnalysisUtils.h"
58 #include "AliPhysicsSelection.h"
59 #include "AliFlowEventSimple.h"
60 #include "AliAnalysisTaskSE.h"
62 
63 
64 using std::cout;
65 using std::endl;
66 
67 
69 
70 
72  fVevent(NULL),
73  fESD(NULL),
74  fAOD(NULL),
75  fPIDResponse(NULL),
76  fMultSelection(NULL),
77  fAnalysisUtil(NULL),
78  fListHist(NULL),
79  mfileFBHijing(NULL),
80  fListFBHijing(NULL),
81  fListNUACorr(NULL),
82  fHistTaskConfigParameters(NULL),
83  fHistPileUpCount(NULL),
84  fHistEtaPtBefore(NULL),
85  fHistEtaPtAfter(NULL),
86  fHistTPCvsGlobalMultBefore(NULL),
87  fHistTPCvsGlobalMultAfter(NULL),
88  fHistTPCdEdxvsPBefore(NULL),
89  fHistTPCdEdxvsPAfter(NULL),
90  fHistTOFBetavsPBefore(NULL),
91  fHistTOFBetavsPAfter(NULL),
92  fHistTOFMassvsPtBefore(NULL),
93  fHistTOFMatchCount(NULL),
94  fHistTPCVsESDTrkBefore(NULL),
95  fHistTPCVsESDTrkAfter(NULL),
96  fHistTPConlyVsCL1Before(NULL),
97  fHistTPConlyVsV0MBefore(NULL),
98  fHistTPConlyVsCL1After(NULL),
99  fHistTPConlyVsV0MAfter(NULL),
100  fHistGlobalVsV0MBefore(NULL),
101  fHistGlobalVsV0MAfter(NULL),
102  fHistRawVsCorrMultFB(NULL),
103  hCentvsTPCmultCuts(NULL),
104  fSkipOutlierCut(0),
105  fFilterBit(1),
106  fNSigmaCut(2),
107  fMinPtCut(0.2),
108  fMaxPtCut(5.0),
109  fMinEtaCut(-0.8),
110  fMaxEtaCut(0.8),
111  fCentralityPercentMin(0),
112  fCentralityPercentMax(90),
113  fPileUpSlopeParm(3.43),
114  fPileUpConstParm(43),
115  bApplyMCcorr(kFALSE),
116  sPathOfMCFile("/alien"),
117  sNucleiTP("PbPb"),
118  fHistEventCount(NULL)
119 {
120  for(int i=0;i<3;i++){
121  fHistPtwithTPCNsigma[i]=NULL;
122  fHistPtwithTOFmasscut[i]=NULL;
123  fHistPtwithTOFSignal[i]=NULL;
124  fHistTOFnSigmavsPtAfter[i]=NULL;
125  fHistTPCnSigmavsPtAfter[i]=NULL;
126  fHistTPCTOFnSigmavsPtAfter[i]=NULL;
127  fHistTPCdEdxvsPtPIDAfter[i]=NULL;
128  }
129  for(int i=0;i<5;i++){
130  fHCorrectNUApos[i] = NULL;
131  fHCorrectNUAneg[i] = NULL;
132  }
133  for(int i=0;i<2;i++){
134  for(int j=0;j<3;j++){
135  fHist_Corr3p_EP_Norm_PN[i][j] = NULL;
136  fHist_Corr3p_EP_Norm_PP[i][j] = NULL;
137  fHist_Corr3p_EP_Norm_NN[i][j] = NULL;
138  fHist_Reso2n_EP_Norm_Det[i][j] = NULL;
139  }
140  }
141  for(int i=0;i<3;i++){
142  for(int j=0;j<5;j++){
143  fHist3DEtaPhiVz_Pos_Run[i][j]=NULL;
144  fHist3DEtaPhiVz_Neg_Run[i][j]=NULL;
145  }
146  }
147  for(int i=0;i<10;i++){
148  fFB_Efficiency_Cent[i] = NULL;
149  }
150 
151 
152  DefineInput(0,TChain::Class());
153  DefineOutput(1,TList::Class());
154 }
155 //______________________________empty constructor_______________________
157  fVevent(NULL),
158  fESD(NULL),
159  fAOD(NULL),
160  fPIDResponse(NULL),
161  fMultSelection(NULL),
162  fAnalysisUtil(NULL),
163  fListHist(NULL),
164  mfileFBHijing(NULL),
165  fListFBHijing(NULL),
166  fListNUACorr(NULL),
167  fHistTaskConfigParameters(NULL),
168  fHistPileUpCount(NULL),
169  fHistEtaPtBefore(NULL),
170  fHistEtaPtAfter(NULL),
171  fHistTPCvsGlobalMultBefore(NULL),
172  fHistTPCvsGlobalMultAfter(NULL),
173  fHistTPCdEdxvsPBefore(NULL),
174  fHistTPCdEdxvsPAfter(NULL),
175  fHistTOFBetavsPBefore(NULL),
176  fHistTOFBetavsPAfter(NULL),
177  fHistTOFMassvsPtBefore(NULL),
178  fHistTOFMatchCount(NULL),
179  fHistTPCVsESDTrkBefore(NULL),
180  fHistTPCVsESDTrkAfter(NULL),
181  fHistTPConlyVsCL1Before(NULL),
182  fHistTPConlyVsV0MBefore(NULL),
183  fHistTPConlyVsCL1After(NULL),
184  fHistTPConlyVsV0MAfter(NULL),
185  fHistGlobalVsV0MBefore(NULL),
186  fHistGlobalVsV0MAfter(NULL),
187  fHistRawVsCorrMultFB(NULL),
188  hCentvsTPCmultCuts(NULL),
189  fSkipOutlierCut(0),
190  fFilterBit(1),
191  fNSigmaCut(2),
192  fMinPtCut(0.2),
193  fMaxPtCut(5.0),
194  fMinEtaCut(-0.8),
195  fMaxEtaCut(0.8),
196  fCentralityPercentMin(0),
197  fCentralityPercentMax(90),
198  fPileUpSlopeParm(3.43),
199  fPileUpConstParm(43),
200  bApplyMCcorr(kFALSE),
201  sPathOfMCFile("/alien"),
202  sNucleiTP("PbPb"),
203  fHistEventCount(NULL)
204 {
205  for(int i=0;i<3;i++){
206  fHistPtwithTPCNsigma[i]=NULL;
207  fHistPtwithTOFmasscut[i]=NULL;
208  fHistPtwithTOFSignal[i]=NULL;
209  fHistTOFnSigmavsPtAfter[i]=NULL;
210  fHistTPCnSigmavsPtAfter[i]=NULL;
212  fHistTPCdEdxvsPtPIDAfter[i]=NULL;
213  }
214  for(int i=0;i<5;i++){
215  fHCorrectNUApos[i] = NULL;
216  fHCorrectNUAneg[i] = NULL;
217  }
218  for(int i=0;i<2;i++){
219  for(int j=0;j<3;j++){
220  fHist_Corr3p_EP_Norm_PN[i][j] = NULL;
221  fHist_Corr3p_EP_Norm_PP[i][j] = NULL;
222  fHist_Corr3p_EP_Norm_NN[i][j] = NULL;
223  fHist_Reso2n_EP_Norm_Det[i][j] = NULL;
224  }
225  }
226  for(int i=0;i<3;i++){
227  for(int j=0;j<5;j++){
228  fHist3DEtaPhiVz_Pos_Run[i][j]=NULL;
229  fHist3DEtaPhiVz_Neg_Run[i][j]=NULL;
230  }
231  }
232  for(int i=0;i<10;i++){
233  fFB_Efficiency_Cent[i] = NULL;
234  }
235 }
236 
237 //___________________________ destructor ___________________________
239 {
240  //Destructor
241  //if(fPIDResponse) delete fPIDResponse;
242  //if(fMultSelection) delete fMultSelection;
243  if(fListHist){
244  delete fListHist;
245  }
246  if(mfileFBHijing->IsOpen()){
247  mfileFBHijing->Close();
248  if(fListFBHijing) delete fListFBHijing;
249  }
250  for(int i=0;i<10;i++){
251  if(fFB_Efficiency_Cent[i])
252  delete fFB_Efficiency_Cent[i];
253  }
254  for(int i=0;i<5;i++){
255  if(fHCorrectNUApos[i]) delete fHCorrectNUApos[i];
256  if(fHCorrectNUAneg[i]) delete fHCorrectNUAneg[i];
257  }
258  if(fAnalysisUtil){
259  delete fAnalysisUtil; // its 'new' !!
260  }
261 }//---------------- sanity ------------------------
262 
263 
264 
266 {
267  //input hander
268  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
269  AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
270  if (!inputHandler) { printf("\n ...Input handler missing!!!...\n"); return; }
271 
272  //PileUp Multi-Vertex
273  fAnalysisUtil = new AliAnalysisUtils();
274  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
275  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
276 
277  //pid response object
278  fPIDResponse=inputHandler->GetPIDResponse();
279 
280  fListHist = new TList();
281  fListHist->SetOwner(kTRUE);
282 
283 
285 
286  if(!gGrid){
287  TGrid::Connect("alien://");
288  }
289 
291 
292 
293 
294  hCentvsTPCmultCuts = new TH2F("hCentvsTPCmultCuts","TPCmult Low,high",100,0,100,5,0,5);
295  hCentvsTPCmultCuts->GetYaxis()->SetBinLabel(1,"mean");
296  hCentvsTPCmultCuts->GetYaxis()->SetBinLabel(2,"sigma");
298 
300 
301  fHistTOFMatchCount = new TH2F("fHistTOFMatchCount","TofMatchFlag vs mismatch Prob",10,0,10,200,-5,5);
303 
304  fHistEtaPtBefore = new TH2F("fHistEtaPtBefore","Eta vs pT",100,-1.25,1.25,100,0,10);
306 
307  fHistEtaPtAfter = new TH2F("fHistEtaPtAfter","Eta vs pT",100,-1.25,1.25,100,0,10);
309 
310 
311  Int_t gMaxTPCFB1mult = 0;
312  Int_t gMaxGlobalmult = 0;
313  Int_t gMaxTPCcorrmult = 0;
314  Int_t gMaxESDtracks = 0;
315 
316  if(sNucleiTP=="pp"||sNucleiTP=="PP"){
317  gMaxGlobalmult = 200;
318  gMaxTPCFB1mult = 200;
319  gMaxTPCcorrmult = 500;
320  gMaxESDtracks = 1000;
321  fSkipOutlierCut = 1;
322  }
323  else if(sNucleiTP=="pPb"||sNucleiTP=="Pbp"||sNucleiTP=="PbP"||sNucleiTP=="PPb"){
324  gMaxGlobalmult = 400;
325  gMaxTPCFB1mult = 400;
326  gMaxTPCcorrmult = 500;
327  gMaxESDtracks = 2000;
328  fSkipOutlierCut = 1;
329  }
330  else{
331  gMaxGlobalmult = 4000;
332  gMaxTPCFB1mult = 4000;
333  gMaxTPCcorrmult = 5000;
334  gMaxESDtracks = 20000;
335  fSkipOutlierCut = 0;
336  }
337 
338 
339 
340  fHistTPCVsESDTrkBefore = new TH2F("fHistTPCVsESDTrkBefore","Before; TPC1; ESD trk",500,0,gMaxTPCcorrmult,200,0,gMaxESDtracks);
342  fHistTPCVsESDTrkAfter = new TH2F("fHistTPCVsESDTrkAfter"," After; TPC1; ESD trk",500,0,gMaxTPCcorrmult,200,0,gMaxESDtracks);
344 
345  fHistTPCvsGlobalMultBefore = new TH2F("fHistTPCvsGlobalMultBefore","Before; Global; TPC(fb1) ",200,0,gMaxGlobalmult,200,0,gMaxTPCFB1mult);
347  fHistTPCvsGlobalMultAfter = new TH2F("fHistTPCvsGlobalMultAfter"," After; Global; TPC(fb1) ",200,0,gMaxGlobalmult,200,0,gMaxTPCFB1mult);
349 
350  fHistGlobalVsV0MBefore = new TH2F("fHistGlobalVsV0MBefore","Before;Cent(V0M);Global",100,0,100,500,0,gMaxGlobalmult);
352  fHistGlobalVsV0MAfter = new TH2F("fHistGlobalVsV0MAfter"," After; Cent(V0M);Global",100,0,100,500,0,gMaxGlobalmult);
354 
355  fHistTPConlyVsCL1Before = new TH2F("fHistTPConlyVsCL1Before","Before;Cent(CL1); TPC(FB1)",100,0,100,250,0,gMaxTPCcorrmult);
357  fHistTPConlyVsCL1After = new TH2F("fHistTPConlyVsCL1After","After; Cent(CL1); TPC(FB1) ",100,0,100,250,0,gMaxTPCcorrmult);
359 
360  fHistTPConlyVsV0MBefore = new TH2F("fHistTPConlyVsV0MBefore","Before;Cent(V0M); TPC(FB1)",100,0,100,250,0,gMaxTPCcorrmult);
362  fHistTPConlyVsV0MAfter = new TH2F("fHistTPConlyVsV0MAfter","After; Cent(V0M); TPC(FB1) ",100,0,100,250,0,gMaxTPCcorrmult);
364 
365  fHistRawVsCorrMultFB = new TH2F("fHistRawVsCorrMultFB",Form("FB%d;Mult_{raw};Mult_{corr}",fFilterBit),gMaxTPCFB1mult,0,gMaxTPCFB1mult,gMaxTPCcorrmult,0,gMaxTPCcorrmult);
367 
368 
369  fHistTPCdEdxvsPBefore = new TH2F("fHistTPCdEdxvsPBefore","Before; p (GeV/c); dEdx (arb)",200,-5,5,200,0,250);
371  fHistTPCdEdxvsPAfter = new TH2F("fHistTPCdEdxvsPAfter"," After; p (GeV/c); dEdx (arb)",200,-5,5, 200,0,250);
373 
374 
375  fHistTOFBetavsPBefore = new TH2F("fHistTOFBetavsPBefore","Before; p (GeV/c); beta ",200,-5,5,100,0.0,1.2);
377  fHistTOFBetavsPAfter = new TH2F("fHistTOFBetavsPAfter"," After; p (GeV/c); beta ",200,-5,5,100,0.0,1.2);
379 
380 
381  fHistTOFMassvsPtBefore = new TH2F("fHistTOFMassvsPtBefore","Before; p_{T}(GeV/c); m^{2}(GeV^{2}/c^{4})",200,-5,5,500,-0.5,4.5);
383 
384 
385  //---------------- PID Histograms ---------------------
386  //Dont forget to add histograms in the List. !!!
387 
388 
389  char const *gSpecies[3] = {"Pion","Kaon","proton"};
390 
391  for(int i=0;i<3;i++){
392  fHistPtwithTPCNsigma[i] = new TH1F(Form("fHistPtwithTPCNsigma_%s",gSpecies[i]), Form("%s;p_{T}(GeV/c))",gSpecies[i]),200,-5,5);
394  fHistPtwithTOFmasscut[i] = new TH1F(Form("fHistPtwithTOFmasscut_%s",gSpecies[i]),Form("%s;p_{T}(GeV/c))",gSpecies[i]),200,-5,5);
396  fHistPtwithTOFSignal[i] = new TH1F(Form("fHistPtwithTOFSignal_%s", gSpecies[i]),Form("%s;p_{T}(GeV/c))",gSpecies[i]),200,-5,5);
398 
399  fHistTOFnSigmavsPtAfter[i] = new TH2F(Form("fHistTOFnSigmavsPtAfter_%s",gSpecies[i]),Form("%s;p_{T}(GeV/c);n#sigma_{TOF}",gSpecies[i]),200,-5,5,400,-10.0,10.0);
401  fHistTPCnSigmavsPtAfter[i] = new TH2F(Form("fHistTPCnSigmavsPtAfter_%s",gSpecies[i]),Form("%s;p_{T}(GeV/c);n#sigma_{TPC}",gSpecies[i]),200,-5,5,400,-10.0,10.0);
403 
404  fHistTPCTOFnSigmavsPtAfter[i] = new TH3F(Form("fHistTPCTOFnSigmavsPtAfter_%s",gSpecies[i]),Form("%s; p_{T}(GeV/c); n#sigma_{TPC}; n#sigma_{TOF}",gSpecies[i]),100,0,5,400,-10,10,400,-10,10);
406 
407  fHistTPCdEdxvsPtPIDAfter[i] = new TH2F(Form("fHistTPCdEdxvsPtAfter_%s",gSpecies[i]),"AfterCut; p_{T} (GeV/c); dEdx (arb)",400,0,10,200,0,250);
409  }// PID histograms done
410 
411 
412 
413  Double_t centRange[11] = {0,5,10,20,30,40,50,60,70,80,90};
414 
415  //------------------- CME 3p correlator Charged hadrons (EP method) ------------------
416  for(int i=0;i<2;i++){
417  for(int j=0;j<3;j++){
418  //Detector: 0 = V0A, 1 = V0C, 3 = Q-cumulant
419  fHist_Corr3p_EP_Norm_PN[i][j] = new TProfile(Form("fHist_Corr3p_EP_Norm_PosNeg_Mag%d_Det%d",i,j+1),"opposit charge correlator",10,centRange,"");
420  //fHist_Corr3p_EP_Norm_PN[i][j]->Sumw2();
422  fHist_Corr3p_EP_Norm_PP[i][j] = new TProfile(Form("fHist_Corr3p_EP_Norm_PosPos_Mag%d_Det%d",i,j+1),"pos-pos charge correlator",10,centRange,"");
423  //fHist_Corr3p_EP_Norm_PP[i][j]->Sumw2();
425  fHist_Corr3p_EP_Norm_NN[i][j] = new TProfile(Form("fHist_Corr3p_EP_Norm_NegNeg_Mag%d_Det%d",i,j+1),"neg-neg charge correlator",10,centRange,"");
426  //fHist_Corr3p_EP_Norm_NN[i][j]->Sumw2();
428  }
429  //EP Resolution:
430  for(int j=0;j<3;j++){
431  //Det: 0 = v0c-v0a, 1 = v0a-TPC, 2 = v0c-TPC,
432  fHist_Reso2n_EP_Norm_Det[i][j] = new TProfile(Form("fHist_Reso2n_EP_Norm_Mag%d_DetComb%d",i,j+1),"Event plane Resolution",10,centRange,"");
433  //fHist_Reso2n_EP_Norm_Det[i][j]->Sumw2();
435  }
436  }
437  //-------------------------------------------------------------------------------
438 
439 
440 
441 
442  //-------------------------- Define NUA Hist for PID -----------------------------
443  Int_t gCentForNUA[6] = {0,5,10,20,40,90};
444  Char_t name[100];
445  Char_t title[100];
446 
447  for(int i=0;i<3;i++){
448  for(int j=0;j<5;j++){
449  sprintf(name,"fHistEtaPhiVz_%s_Pos_Cent%d_Run%d",gSpecies[i],j,1);
450  sprintf(title,"eta,phi,Vz %sPos, Cent%d-%d, FB %d",gSpecies[i],gCentForNUA[j],gCentForNUA[j+1],fFilterBit);
451  fHist3DEtaPhiVz_Pos_Run[i][j] = new TH3F(name,title,10,-10,10,50,0,6.283185,16,-0.8,0.8);
453 
454  sprintf(name,"fHistEtaPhiVz_%s_Neg_Cent%d_Run%d",gSpecies[i],j,1);
455  sprintf(title,"eta,phi,Vz %sNeg, Cent%d-%d, FB %d",gSpecies[i],gCentForNUA[j],gCentForNUA[j+1],fFilterBit);
456  fHist3DEtaPhiVz_Neg_Run[i][j] = new TH3F(name,title,10,-10,10,50,0,6.283185,16,-0.8,0.8);
458  }
459  }
460  //---------------------------------------------------------------------------------
461 
462 
463 
464 
465  PostData(1,fListHist);
466  std::cout<<"........ UserCreateOutputObject called ... fFilterBit = "<<fFilterBit<<" Cmax = "<<fCentralityPercentMax <<"\n"<<std::endl;
467 }
468 
469 
470 
471 
472 
473 
474 
475 //______________________________________________________________________
477  //printf("info: UserExec is called.... 1 \n");
478 
479  Float_t stepCount = 0.5;
480 
481  fHistEventCount->Fill(stepCount); //1
482  stepCount++;
483 
484  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
485  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
486 
487  if(!(fESD || fAOD)){ printf("ERROR: fESD & fAOD not available\n"); return; }
488 
489  fVevent = dynamic_cast<AliVEvent*>(InputEvent());
490 
491  if (!fVevent) { printf("ERROR: fVevent not available\n"); return; }
492 
493  fHistEventCount->Fill(stepCount); //2
494  stepCount++;
495 
496 
497 
498  //--------- Check if I have PID response object --------
499  if(!fPIDResponse){
500  printf("\n\n...... PIDResponse object not found..... \n\n");
501  return;
502  }
503 
504 
505 
506  //-------------- Vtx cuts ---------------
507  const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
508 
509  Double_t pVtxZ = -999;
510  pVtxZ = pVtx->GetZ();
511 
512  if(TMath::Abs(pVtxZ)>10.) return;
513 
514  fHistEventCount->Fill(stepCount); //3
515  stepCount++;
516 
517  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
518 
519  if(!fMultSelection) { printf("\n\n **WARNING** \n::UserExec() AliMultSelection object not found.\n\n"); exit(1); }
520 
521  Float_t centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
522  Float_t centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
523 //Float_t centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
524 //Float_t centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
525 
526 
527 
528  Float_t centrality = centrV0M;
529 
531  return;
532  }
533 
534  fHistEventCount->Fill(stepCount); //4
535  stepCount++;
536 
537 
538  Int_t ntracks=fAOD->GetNumberOfTracks();
539 
540  if(ntracks<2) return; // Check this cut....!!!
541 
542 
543  fHistEventCount->Fill(stepCount); //5
544  stepCount++;
545 
546 
547  Int_t cent10bin = -1;
548 
549  cent10bin = GetCentralityScaled0to10(centrality); //Centrality in 0-10 scale
550 
551 
552 
553  Int_t cForNUA = 0; //Centrality array index for NUA correcion
554 
555  if(centrality<5.0) {
556  cForNUA = 0;
557  }
558  else if(centrality>=5.0 && centrality<10){
559  cForNUA = 1; // 1=5-10,
560  }
561  else if(centrality>=10.0 && centrality<20) {
562  cForNUA = 2; // 2 = 10-20,
563  }
564  else if(centrality>=20 && centrality<40){
565  cForNUA = 3; // 3=20-40
566  }
567  else if(centrality>=40){
568  cForNUA = 4; // 4=40-90
569  }
570 
571 
572 
573  //Variables for MC tracking correction
574  Int_t ptBinMC = 1;
575  Float_t ptWgtMC = 1.0;
576  Float_t ptTrk = 0.1;
577 
578 
579  //-------------- Track loop for outlier and PileUp cut -------------------
580 
581  Bool_t bIsPileup=kFALSE;
582 
583  Int_t isPileup = fAOD->IsPileupFromSPD(3);
584 
585  if(isPileup != 0) {
586  fHistPileUpCount->Fill(0.5);
587  bIsPileup=kTRUE;
588  }
589  else if(PileUpMultiVertex(fAOD)) {
590  fHistPileUpCount->Fill(1.5);
591  bIsPileup=kTRUE;
592  }
593  else if(((AliAODHeader*)fAOD->GetHeader())->GetRefMultiplicityComb08() < 0) {
594  fHistPileUpCount->Fill(2.5);
595  bIsPileup=kTRUE;
596  }
597  else if(fAOD->IsIncompleteDAQ()) {
598  fHistPileUpCount->Fill(3.5);
599  bIsPileup=kTRUE;
600  }
601  else if(fabs(centrV0M-centrCL1)> 5.0) {//default: 7.5
602  fHistPileUpCount->Fill(4.5);
603  bIsPileup=kTRUE;
604  }
605 
606  // check vertex consistency
607  const AliAODVertex* vtTrc = fAOD->GetPrimaryVertex();
608  const AliAODVertex* vtSPD = fAOD->GetPrimaryVertexSPD();
609 
610  if(vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
611  fHistPileUpCount->Fill(5.5);
612  bIsPileup=kTRUE;
613  }
614 
615  double covTrc[6], covSPD[6];
616  vtTrc->GetCovarianceMatrix(covTrc);
617  vtSPD->GetCovarianceMatrix(covSPD);
618 
619  double dz = vtTrc->GetZ() - vtSPD->GetZ();
620 
621  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
622  double errTrc = TMath::Sqrt(covTrc[5]);
623  double nsigTot = dz/errTot;
624  double nsigTrc = dz/errTrc;
625 
626  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
627  fHistPileUpCount->Fill(6.5);
628  bIsPileup=kTRUE;
629  }
630 
631  //---------------- a dobrin --------------
632 
633  Bool_t bIsOutLier=kTRUE;
634 
635  Float_t multTPC = 0; // tpc mult estimate
636  Float_t RefMultRaw = 0; // tpc mult estimate
637  Float_t RefMultCorr = 0; // tpc mult estimate
638  Float_t RefMultRawFB = 0;
639  Float_t RefMultCorrFB= 0;
640 
641  Float_t multTPCAll = 0; // tpc mult estimate
642  Float_t multGlobal = 0; // global multiplicity
643 
644 
645 //const Int_t nGoodTracks = fVevent->GetNumberOfTracks();
646  const Int_t nGoodTracks = fAOD->GetNumberOfTracks();
647 
648  for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { //-------------------------
649 
650  AliAODTrack* AODtrack =dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
651  if(!AODtrack) continue;
652 
653  if(AODtrack->TestFilterBit(128)) multTPCAll++; //A. Dobrin, no track cuts
654 
655  if(!(AODtrack->TestFilterBit(1))) continue;
656  if((AODtrack->Pt() < .2) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.2)) continue;
657  multTPC++;
658 
659  ptTrk = AODtrack->Pt();
660 
661  if(ptTrk >= 0.2 && ptTrk <= 5.0) {
662  if(fFB_Efficiency_Cent[cent10bin]){
663  ptBinMC = fFB_Efficiency_Cent[cent10bin]->FindBin(ptTrk);
664  ptWgtMC = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBinMC);
665  //cout<<iTrack<<"cent = "<<cent10bin<<" pt = "<<ptTrk<<"\t bin = "<<ptBinMC<<"\t Wgt = "<<ptWgtMC<<endl;
666  RefMultRaw++;
667  RefMultCorr += ptWgtMC;
668  if((AODtrack->TestFilterBit(fFilterBit))){
669  RefMultRawFB++;
670  RefMultCorrFB += ptWgtMC;
671  }
672  }
673  }
674 
675  Double_t b[2] = {-99., -99.};
676  Double_t bCov[3] = {-99., -99., -99.};
677  AliAODTrack copy(*AODtrack);
678  if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
679  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
680  multGlobal++;
681 
682  }//--- track loop outlier/PileUp ----
683 
684  Int_t multEsd = ((AliAODHeader*)fAOD->GetHeader())->GetNumberOfESDTracks();
685 
686  Float_t multESDTPCDiff = (Float_t) multEsd - fPileUpSlopeParm*multTPCAll;
687 
688  if(multESDTPCDiff > fPileUpConstParm) {
689  fHistPileUpCount->Fill(7.5);
690  bIsPileup=kTRUE;
691  }
692  else if(bIsPileup==kFALSE) {
693  if(!fMultSelection->GetThisEventIsNotPileup()) bIsPileup=kTRUE;
694  if(!fMultSelection->GetThisEventIsNotPileupMV()) bIsPileup=kTRUE;
695  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()) bIsPileup=kTRUE;
696  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()) bIsPileup=kTRUE;
697  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()) bIsPileup=kTRUE;
698  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()) bIsPileup=kTRUE;
699  if(!fMultSelection->GetThisEventHasGoodVertex2016()) bIsPileup=kTRUE;
700  if(bIsPileup) fHistPileUpCount->Fill(9.5);
701  }
702  //-----------------------------------------------------------------
703 
704 
705 
706  fHistTPCVsESDTrkBefore->Fill(multTPCAll,multEsd); //A. Dobrin
707 
708  fHistTPCvsGlobalMultBefore->Fill(multGlobal,multTPC);
709 
710 
711 
712 //if(multTPC < (-40.3+1.22*multGlobal) || multTPC > (32.1+1.59*multGlobal)) { bIsOutLier = kFALSE;} //Run-1 or 2011
713  if(multTPC < (-20.0+1.15*multGlobal) || multTPC > (200.+1.40*multGlobal)) { bIsOutLier = kFALSE;}
714 
715  fHistEventCount->Fill(stepCount); //6
716  stepCount++;
717 
718 
719  fHistTPConlyVsCL1Before->Fill(centrCL1,multTPCAll);
720  fHistTPConlyVsV0MBefore->Fill(centrV0M,multTPCAll);
721  fHistGlobalVsV0MBefore->Fill(centrV0M, multGlobal);
722 
723 
724 
725 
726 
727 
728 
729  if(!fSkipOutlierCut && bIsOutLier) return; //outlier TPC vs Global
730 
731  fHistTPCvsGlobalMultAfter->Fill(multGlobal,multTPC);
732 
733  fHistEventCount->Fill(stepCount); //7
734  stepCount++;
735 
736 
737  if(!fSkipOutlierCut && bIsPileup) return; //PileUp A. Dobrin
738 
739  fHistTPCVsESDTrkAfter->Fill(multTPCAll,multEsd);
740 
741  fHistEventCount->Fill(stepCount); //8
742  stepCount++;
743 
744 
745 
746 
747 
748 
749 
750 
751 
752 
753  Int_t icentBin = centrality;
754  icentBin++;
755 
756  Float_t TPCmultLowLimit = hCentvsTPCmultCuts->GetBinContent(icentBin,1);
757  Float_t TPCmultHighLimit = hCentvsTPCmultCuts->GetBinContent(icentBin,1);
758 
759  TPCmultLowLimit -= 5.0 * hCentvsTPCmultCuts->GetBinContent(icentBin,2); //mean - 5sigma
760  TPCmultHighLimit += 5.0 * hCentvsTPCmultCuts->GetBinContent(icentBin,2); //mean + 5sigma
761  //std::cout<<" Cent = "<<centrality<<"\t icent = "<<icentBin<<" low = "<<TPCmultLowLimit<<"\t high = "<<TPCmultHighLimit<<std::endl;
762 
763  if(!fSkipOutlierCut){
764  if(multTPC<TPCmultLowLimit || multTPC>TPCmultHighLimit) return;
765  }
766 
767  fHistEventCount->Fill(stepCount); //9
768  stepCount++;
769 
770 
771 
772  fHistTPConlyVsCL1After->Fill(centrCL1,multTPCAll);
773  fHistTPConlyVsV0MAfter->Fill(centrV0M,multTPCAll);
774  fHistGlobalVsV0MAfter->Fill(centrV0M, multGlobal);
775 
776 
777  // MC corrected Refmult:
778  fHistRawVsCorrMultFB->Fill(RefMultRawFB,RefMultCorrFB); // FB set by AddTask..
779 
780  //--------------------------------------------------------
781 
782 
783 
784 
785 
786 
787 
788 
789 
790 
791 
792 
793 
794 
795 
796 
797  //if(nGoodTracks!=ntracks) std::cout<<" Event with ntracks = "<<ntracks<<"\t good Tracks = "<<nGoodTracks<<std::endl;
798 
799 
800 
801 
802 
803 
804 
805 
806 
807 
808 
809 
810 
811 
812 
813 
814 
815 
816 
817 
818 
819 
820 
821 
822  //--------- PID works begin ----------
823  Double_t PDGmassPion = 0.13957;
824  Double_t PDGmassKaon = 0.49368;
825  Double_t PDGmassProton = 0.93827;
826 
827  PDGmassProton *= PDGmassProton;
828  PDGmassPion *= PDGmassPion;
829  PDGmassKaon *= PDGmassKaon;
830 
831  Double_t mass=0,mom = -999, pT = -999, phi = -999, eta = -999, dEdx =-999;
832  Double_t length = -999., beta =-999, tofTime = -999., tof = -999.;
833  Double_t dPhi1,dPhi2,dPt1,ptw1,dPt2,ptw2,dEta1,dEta2;
834  Double_t nSigTOFpion, nSigTPCpion;
835  Double_t nSigTOFkaon, nSigTPCkaon;
836  Double_t nSigTOFproton, nSigTPCproton;
837  Double_t c = TMath::C()*1.E-9;// velocity of light m/ns
838  Float_t dcaXY, dcaZ, WgtEP, w1NUA, w2NUA;
839  Float_t probMis;
840  Int_t TOFmatch=0;
841  Int_t charge,ptBin,iBinNUA,dChrg1,dChrg2;
842 
843 
844  //--------- Track Loop for PID studies ----------------
845 
846  for(Int_t itrack = 0; itrack < ntracks; itrack++) {
847 
848  AliAODTrack *track=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(itrack));
849  if(!track) continue;
850 
851  dcaXY = track->DCA();
852  dcaZ = track->ZAtDCA();
853 
854  if(!track->TestFilterBit(fFilterBit)) continue;
855 
856  mom=track->P();
857  pT=track->Pt();
858  phi=track->Phi();
859  eta=track->Eta();
860  dEdx=track->GetDetPid()->GetTPCsignal();
861  charge = track->Charge();
862 
863 
864 
865 
866 
867  //-------------- Check TOF status ------------------
868  AliPIDResponse::EDetPidStatus status;
869  status = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,track);
870  TOFmatch = 0;
871  if(status==AliPIDResponse::kDetPidOk){
872  TOFmatch++;
873  }
874  probMis = fPIDResponse->GetTOFMismatchProbability(track);
875  fHistTOFMatchCount->Fill(TOFmatch,probMis);
876 
877  mass = -9.9;
878  beta = -0.5;
879 
880  if(TOFmatch>0 && probMis < 0.01) {
881  //This conditions are called when detector status is checked above :
882  //if((track->IsOn(AliAODTrack::kTOFin)) && (track->IsOn(AliAODTrack::kTIME)) && (track->IsOn(AliAODTrack::kTOFout))) {
883  //if((track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid))) { //Naghmeh used it
884  tofTime = track->GetTOFsignal(); // in pico seconds
885  length = track->GetIntegratedLength();
886  tof = tofTime*1.E-3; // ns
887  if (tof <= 0) tof = 9999;
888  length = length*0.01; // in meters
889  tof = tof*c;
890  beta = length/tof;
891  mass = mom*mom*(1./(beta*beta) - 1);
892  }//------------ TOF signal -------------------------
893 
894 
895  //QA histograms:
896  fHistEtaPtBefore->Fill(eta,pT);
897  fHistTPCdEdxvsPBefore->Fill(mom*charge,dEdx);
898  fHistTOFBetavsPBefore->Fill(mom*charge,beta);
899  fHistTOFMassvsPtBefore->Fill(pT*charge,mass);
900 
901 
902 
903 
904  //-------- Apply Default track cuts for analysis: ---------
905  if(pT < fMinPtCut || pT > fMaxPtCut) continue;
906  if(eta < fMinEtaCut || eta > fMaxEtaCut) continue;
907  if(!(track->TestFilterBit(fFilterBit))) continue;
908  //-----------------------------------------------------
909 
910 
911 
912 
913 
914  //============= charged hadron analysis: ============
915 
916  dPt1 = pT;
917  dPhi1 = phi;
918  dEta1 = eta;
919  dChrg1=charge;
920  //------ get MC weight and NUA for track 1--------------
921  if(fFB_Efficiency_Cent[cent10bin]){
922  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt1);
923  ptw1 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
924  }
925  else{ ptw1 = 1.0; }
926 
927 
928  if(dChrg1>0){
929  if(fHCorrectNUApos[cForNUA]){
930  iBinNUA = fHCorrectNUApos[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
931  w1NUA = fHCorrectNUApos[cForNUA]->GetBinContent(iBinNUA);
932  }
933  else{ w1NUA = 1.0; }
934  }
935  else{
936  if(fHCorrectNUAneg[cForNUA]){
937  iBinNUA = fHCorrectNUAneg[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
938  w1NUA = fHCorrectNUAneg[cForNUA]->GetBinContent(iBinNUA);
939  }
940  else{ w1NUA = 1.0; }
941  }
942 
943  //-----------------------------------------------------------
944 
945 
946  /*
947  //---2nd track loop (nested)---
948  for(Int_t jtrack = 0; jtrack < ntracks; jtrack++) {
949 
950  if(jtrack==itrack) continue;
951 
952  AliAODTrack *track2=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(jtrack));
953  if(!track2) continue;
954 
955  if(!track2->TestFilterBit(fFilterBit)) continue;
956 
957  dPt2 = track2->Pt();
958  dPhi2 = track2->Phi();
959  dEta2 = track2->Eta();
960  dChrg2= track2->Charge();
961 
962  //------ get MC weight and NUA for track 2--------------
963  if(fFB_Efficiency_Cent[cent10bin]){
964  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt2);
965  ptw2 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
966  }
967  else{ ptw1 = 1.0; }
968 
969  if(dChrg2>0){
970  if(fHCorrectNUApos[cForNUA]){
971  iBinNUA = fHCorrectNUApos[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
972  w2NUA = fHCorrectNUApos[cForNUA]->GetBinContent(iBinNUA);
973  }
974  else{ w2NUA = 1.0; }
975  }
976  else{
977  if(fHCorrectNUAneg[cForNUA]){
978  iBinNUA = fHCorrectNUAneg[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
979  w2NUA = fHCorrectNUAneg[cForNUA]->GetBinContent(iBinNUA);
980  }
981  else{ w2NUA = 1.0; }
982  }
983 
984  WgtEP = ptw1*ptw2*w1NUA*w2NUA;
985  //-----------------------------------------------------------
986 
987 
988 
989  if(dChrg1!=dChrg2) {
990  fHist_Corr3p_EP_Norm_PN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*Psi2V0A),WgtEP);
991  fHist_Corr3p_EP_Norm_PN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*Psi2V0C),WgtEP);
992  }
993  else if(dChrg1>0 && dChrg2>0 && skipPairHBT==0) {
994  fHist_Corr3p_EP_Norm_PP[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*Psi2V0A),WgtEP);
995  fHist_Corr3p_EP_Norm_PP[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*Psi2V0C),WgtEP);
996  }
997 
998  else if(dChrg1<0 && dChrg2<0 && skipPairHBT==0){
999  fHist_Corr3p_EP_Norm_NN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*Psi2V0A),WgtEP);
1000  fHist_Corr3p_EP_Norm_NN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*Psi2V0C),WgtEP);
1001  }
1002 
1003 
1004 
1005  }//-------- 2nd track loop ends ------------------
1006 
1007  */
1008 
1009 
1010 
1011 
1012 
1013 
1014 
1015 
1016 
1017 
1018 
1019 
1020 
1021 
1022 
1023 
1024  //============ PID business starts here =============
1025 
1026  fHistEtaPtAfter->Fill(eta,pT);
1027  fHistTPCdEdxvsPAfter->Fill(mom*charge,dEdx);
1028 
1029  if(TOFmatch>0 && probMis < 0.01){
1030  fHistTOFBetavsPAfter->Fill(mom*charge,beta);
1031  }
1032 
1033  nSigTOFpion=fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
1034  nSigTOFkaon=fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
1035  nSigTOFproton=fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
1036 
1037  fHistTOFMatchCount->Fill(TOFmatch+2,nSigTOFpion);
1038  fHistTOFMatchCount->Fill(TOFmatch+4,nSigTOFkaon);
1039  fHistTOFMatchCount->Fill(TOFmatch+6,nSigTOFproton);
1040 
1041  if(!TOFmatch || probMis > 0.01){ // I dont want mismatched track in my signal distribution
1042  nSigTOFpion = -9.99;
1043  nSigTOFkaon = -9.99;
1044  nSigTOFproton = -9.99;
1045  }
1046 
1047  nSigTPCpion = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
1048  nSigTPCkaon = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
1049  nSigTPCproton = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
1050 
1051  //0=pi, 1=K, 2=Proton
1052  fHistTPCTOFnSigmavsPtAfter[0]->Fill(pT,nSigTPCpion,nSigTOFpion);
1053  fHistTPCTOFnSigmavsPtAfter[1]->Fill(pT,nSigTPCkaon,nSigTOFkaon);
1054  fHistTPCTOFnSigmavsPtAfter[2]->Fill(pT,nSigTPCproton,nSigTOFproton);
1055 
1056 
1057 
1058 
1059 
1060  /*
1061  if(TMath::Abs(nSigTPCpion)<=fNSigmaCut){
1062  fHistTOFnSigmavsPtAfter[0]->Fill(pT*charge,nSigTOFpion);
1063  fHistPtwithTPCNsigma[0]->Fill(pT*charge);
1064  fHistTPCdEdxvsPtPIDAfter[0]->Fill(pT,dEdx);
1065  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
1066  fHistPtwithTOFSignal[0]->Fill(pT*charge);
1067  }
1068  }
1069  if(TMath::Abs(nSigTPCkaon)<=fNSigmaCut){
1070  fHistPtwithTPCNsigma[1]->Fill(pT*charge);
1071  fHistTOFnSigmavsPtAfter[1]->Fill(pT*charge,nSigTOFkaon);
1072  fHistTPCdEdxvsPtPIDAfter[1]->Fill(pT,dEdx);
1073  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
1074  fHistPtwithTOFSignal[1]->Fill(pT*charge);
1075  }
1076  }
1077  if(TMath::Abs(nSigTPCproton)<=fNSigmaCut){
1078  fHistPtwithTPCNsigma[2]->Fill(pT*charge);
1079  fHistTOFnSigmavsPtAfter[2]->Fill(pT*charge,nSigTOFproton);
1080  fHistTPCdEdxvsPtPIDAfter[2]->Fill(pT,dEdx);
1081  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
1082  fHistPtwithTOFSignal[2]->Fill(pT*charge);
1083  }
1084  }
1085  */
1086 
1087 
1088 
1089  // nSigmaTOF distribution for circular cut
1090 
1091  if(TMath::Sqrt(nSigTPCpion*nSigTPCpion+nSigTOFpion*nSigTOFpion)<=fNSigmaCut){
1092  fHistTOFnSigmavsPtAfter[0]->Fill(pT*charge,nSigTOFpion);
1093  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
1094  fHistPtwithTOFSignal[0]->Fill(pT*charge);
1095  }
1096  }
1097 
1098  if(TMath::Sqrt(nSigTPCkaon*nSigTPCkaon+nSigTOFkaon*nSigTOFkaon)<=fNSigmaCut){
1099  fHistTOFnSigmavsPtAfter[1]->Fill(pT*charge,nSigTOFkaon);
1100  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
1101  fHistPtwithTOFSignal[1]->Fill(pT*charge);
1102  }
1103  }
1104 
1105  if(TMath::Sqrt(nSigTPCproton*nSigTPCproton+nSigTOFproton*nSigTOFproton)<=fNSigmaCut){
1106  fHistTOFnSigmavsPtAfter[2]->Fill(pT*charge,nSigTOFproton);
1107  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
1108  fHistPtwithTOFSignal[2]->Fill(pT*charge);
1109  }
1110  }
1111 
1112 
1113 
1114 
1115 
1116  if(TMath::Abs(nSigTPCpion)<=fNSigmaCut){
1117  fHistPtwithTPCNsigma[0]->Fill(pT*charge);
1118  fHistTPCdEdxvsPtPIDAfter[0]->Fill(pT,dEdx);
1119  }
1120  if(TMath::Abs(nSigTPCkaon)<=fNSigmaCut){
1121  fHistPtwithTPCNsigma[1]->Fill(pT*charge);
1122  fHistTPCdEdxvsPtPIDAfter[1]->Fill(pT,dEdx);
1123  }
1124  if(TMath::Abs(nSigTPCproton)<=fNSigmaCut){
1125  fHistPtwithTPCNsigma[2]->Fill(pT*charge);
1126  fHistTPCdEdxvsPtPIDAfter[2]->Fill(pT,dEdx);
1127  }
1128 
1129 
1130 
1131  //============== Fill NUA Histograms for Pion ---------------------
1132  if(pT<0.6 && TMath::Abs(nSigTPCpion)<=2.5){
1133  if(charge>0){
1134  fHist3DEtaPhiVz_Pos_Run[0][cForNUA]->Fill(pVtxZ,phi,eta);
1135  }
1136  else if(charge<0){
1137  fHist3DEtaPhiVz_Neg_Run[0][cForNUA]->Fill(pVtxZ,phi,eta);
1138  }
1139  }
1140  else if(pT>=0.6 && pT<=2.0 && TMath::Abs(nSigTPCpion)<=2.5 && TMath::Abs(nSigTOFpion)<=2.0 ){
1141  if(charge>0){
1142  fHist3DEtaPhiVz_Pos_Run[0][cForNUA]->Fill(pVtxZ,phi,eta);
1143  }
1144  else if(charge<0){
1145  fHist3DEtaPhiVz_Neg_Run[0][cForNUA]->Fill(pVtxZ,phi,eta);
1146  }
1147  }
1148  /*
1149  else if(pT>=2.0 && pT<=4.0 && TMath::Abs(nSigTPCpion)<=1.0 && TMath::Abs(nSigTOFpion)<=1.0){
1150  if(charge>0){
1151  fHist3DEtaPhiVz_Pos_Run[0][cForNUA]->Fill(pVtxZ,phi,eta);
1152  }
1153  else if(charge<0){
1154  fHist3DEtaPhiVz_Neg_Run[0][cForNUA]->Fill(pVtxZ,phi,eta);
1155  }
1156  }*/
1157 
1158 
1159 
1160  //============== Fill NUA Histograms for Kaon ---------------------
1161  if(pT<0.6 && TMath::Abs(nSigTPCkaon)<=2.5){
1162  if(charge>0){
1163  fHist3DEtaPhiVz_Pos_Run[1][cForNUA]->Fill(pVtxZ,phi,eta);
1164  }
1165  else if(charge<0){
1166  fHist3DEtaPhiVz_Neg_Run[1][cForNUA]->Fill(pVtxZ,phi,eta);
1167  }
1168  }
1169  else if(pT>=0.6 && pT<=2.0 && TMath::Abs(nSigTPCkaon)<=2.5 && TMath::Abs(nSigTOFkaon)<=2.0){
1170  if(charge>0){
1171  fHist3DEtaPhiVz_Pos_Run[1][cForNUA]->Fill(pVtxZ,phi,eta);
1172  }
1173  else if(charge<0){
1174  fHist3DEtaPhiVz_Neg_Run[1][cForNUA]->Fill(pVtxZ,phi,eta);
1175  }
1176  }
1177  /*
1178  else if(pT>=3.0 && pT<=4.0 && TMath::Abs(nSigTPCkaon)<=1.0 && TMath::Abs(nSigTOFkaon)<=1.0){
1179  if(charge>0){
1180  fHist3DEtaPhiVz_Pos_Run[1][cForNUA]->Fill(pVtxZ,phi,eta);
1181  }
1182  else if(charge<0){
1183  fHist3DEtaPhiVz_Neg_Run[1][cForNUA]->Fill(pVtxZ,phi,eta);
1184  }
1185  }*/
1186 
1187 
1188 
1189  //============== Fill NUA Histograms for proton ---------------------
1190  if(pT<0.9 && TMath::Abs(nSigTPCproton)<=2.5){
1191  if(charge>0){
1192  fHist3DEtaPhiVz_Pos_Run[2][cForNUA]->Fill(pVtxZ,phi,eta);
1193  }
1194  else if(charge<0){
1195  fHist3DEtaPhiVz_Neg_Run[2][cForNUA]->Fill(pVtxZ,phi,eta);
1196  }
1197  }
1198  else if(pT>=0.9 && pT<=3.5 && TMath::Abs(nSigTPCproton)<=2.5 && TMath::Abs(nSigTOFproton)<=2.5){
1199  if(charge>0){
1200  fHist3DEtaPhiVz_Pos_Run[2][cForNUA]->Fill(pVtxZ,phi,eta);
1201  }
1202  else if(charge<0){
1203  fHist3DEtaPhiVz_Neg_Run[2][cForNUA]->Fill(pVtxZ,phi,eta);
1204  }
1205  }
1206  /*
1207  else if(pT>=3.5 && pT<=4.0 && TMath::Abs(nSigTPCproton)<=1.5 && TMath::Abs(nSigTOFproton)<=1.0){
1208  if(charge>0){
1209  fHist3DEtaPhiVz_Pos_Run[2][cForNUA]->Fill(pVtxZ,phi,eta);
1210  }
1211  else if(charge<0){
1212  fHist3DEtaPhiVz_Neg_Run[2][cForNUA]->Fill(pVtxZ,phi,eta);
1213  }
1214  }*/
1215 
1216 
1217 
1218  if(!TOFmatch || probMis > 0.01 || beta>0.2) continue;
1219 
1220  // nSigmaTPC distribution for Fixed nSigmaTOF cut
1221 
1222  if(TMath::Abs(nSigTOFpion)<=fNSigmaCut){
1223  fHistTPCnSigmavsPtAfter[0]->Fill(pT*charge,nSigTPCpion);
1224  fHistPtwithTOFmasscut[0]->Fill(pT*charge);
1225  }
1226  if(TMath::Abs(nSigTOFkaon)<=fNSigmaCut){
1227  fHistTPCnSigmavsPtAfter[1]->Fill(pT*charge,nSigTPCkaon);
1228  fHistPtwithTOFmasscut[1]->Fill(pT*charge);
1229  }
1230  if(TMath::Abs(nSigTOFproton)<=fNSigmaCut){
1231  fHistTPCnSigmavsPtAfter[2]->Fill(pT*charge,nSigTPCproton);
1232  fHistPtwithTOFmasscut[2]->Fill(pT*charge);
1233  }
1234 
1235  }//track loop ends
1236 
1237 
1238 
1239 
1240 
1241 
1242  PostData(1,fListHist);
1243 
1244  fHistEventCount->Fill(stepCount); //10
1245  stepCount++;
1246 
1247 }//================ UserExec ==============
1248 
1249 
1250 
1251 
1252 
1253 
1254 
1255 double AliAnalysisTaskCMEV0PID::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1256 {
1257  // calculate sqrt of weighted distance to other vertex
1258  if (!v0 || !v1) {
1259  AliDebug(2,"\n\n ::GetWDist => One of vertices is not valid\n\n");
1260  return 0;
1261  }
1262  static TMatrixDSym vVb(3);
1263  double dist = -1;
1264  double dx = v0->GetX()-v1->GetX();
1265  double dy = v0->GetY()-v1->GetY();
1266  double dz = v0->GetZ()-v1->GetZ();
1267  double cov0[6],cov1[6];
1268  v0->GetCovarianceMatrix(cov0);
1269  v1->GetCovarianceMatrix(cov1);
1270  vVb(0,0) = cov0[0]+cov1[0];
1271  vVb(1,1) = cov0[2]+cov1[2];
1272  vVb(2,2) = cov0[5]+cov1[5];
1273  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1274  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1275  vVb.InvertFast();
1276  if (!vVb.IsValid()) {
1277  AliDebug(2,"Singular Matrix\n");
1278  return dist;
1279  }
1280  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1281  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1282  return dist>0 ? TMath::Sqrt(dist) : -1;
1283 }
1284 
1285 
1287  { // check for multi-vertexer pile-up
1288  const int kMinPlpContrib = 5;
1289  const double kMaxPlpChi2 = 5.0;
1290  const double kMinWDist = 15;
1291 
1292  const AliVVertex* vtPrm = 0;
1293  const AliVVertex* vtPlp = 0;
1294 
1295  int nPlp = 0;
1296 
1297  if(!(nPlp=faod->GetNumberOfPileupVerticesTracks()))
1298  return kFALSE;
1299 
1300  vtPrm = faod->GetPrimaryVertex();
1301  if(vtPrm == faod->GetPrimaryVertexSPD())
1302  return kTRUE; // there are pile-up vertices but no primary
1303 
1304  //int bcPrim = vtPrm->GetBC();
1305 
1306  for(int ipl=0;ipl<nPlp;ipl++) {
1307  vtPlp = (const AliVVertex*)faod->GetPileupVertexTracks(ipl);
1308  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1309  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1310  //int bcPlp = vtPlp->GetBC();
1311  //if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2)
1312  // return kTRUE; // pile-up from other BC
1313 
1314  double wDst = GetWDist(vtPrm,vtPlp);
1315  if (wDst<kMinWDist) continue;
1316 
1317  return kTRUE; // pile-up: well separated vertices
1318  }
1319  return kFALSE;
1320 }
1321 
1322 
1323 
1325 
1326  if(bApplyMCcorr){
1327  if(!gGrid){
1328  TGrid::Connect("alien://");
1329  }
1330 
1331  if(!mfileFBHijing){
1332 
1333  mfileFBHijing = TFile::Open(sMCfilePath,"READ");
1334 
1335  fListFBHijing = dynamic_cast<TList*>(mfileFBHijing->FindObjectAny("fMcEffiHij"));
1336 
1337  if(!fListFBHijing){
1338  std::cout<<"\n\n !!!!**** Warning: FB Efficiency File/List not found *****\n\n"<<std::endl;
1339  //exit(1);
1340  }
1341  else if(fListFBHijing) {
1342  for(int i=0;i<10;i++) {
1343  fFB_Efficiency_Cent[i] = (TH1D *) fListFBHijing->FindObject(Form("eff_unbiased_%d",i));
1344  //std::cout<<" input MC hist"<<i<<" = "<<fFB_Efficiency_Cent[i]->GetName()<<std::endl;
1345  }
1346  }
1347  }
1348  }
1349  else{ // if MC efficiency Not used/ file not found, then use weight = 1.
1350  for(int i=0;i<10;i++){
1351  fFB_Efficiency_Cent[i] = new TH1D(Form("eff_unbiased_%d",i),"",1,0,50.);
1352  fFB_Efficiency_Cent[i]->SetBinContent(1,1.0);
1353  }
1354  if(bApplyMCcorr){ printf("\n\n!!***** Warning *****!!\n MC correction File not found, using Efficiency = 1.0 !!\n\n");}
1355  }
1356 }
1357 
1358 
1359 
1360 
1362 
1363  Int_t cIndex = 0;
1364 
1365  if(fCent<5.0) {
1366  cIndex = 0;
1367  }
1368  else if(fCent>=5.0 && fCent<10){
1369  cIndex = 1;
1370  }
1371  else if(fCent>=10.0) {
1372  cIndex = abs(fCent/10.0)+1;
1373  }
1374  return cIndex;
1375 }
1376 
1378  //std::cout<<" centrality outlier function called "<<std::endl;
1379  Float_t fMeanTPC[100] = {2902.95,2758.33,2642.78,2536.67,2435.37,2340.06,2248.44,2163.71,2080.49,2001.54,1925.86,1852.64,1781.97,1715.56,1650.53,1587.23,1527.51,1468.19,1412.73,1357.86,1305.35,1254.33,1205.57,1157.28,1111.53,1066.42,1023.15,981.594,940.795,901.766,863.651,826.183,790.53,756.358,722.654,690.513,659.443,628.807,599.748,571.664,544.446,518.042,492.369,468.072,444.694,422.487,400.104,379.129,359.147,339.62,320.817,302.788,285.791,269.015,253.688,238.671,224.039,209.932,196.915,184.647,172.76,161.381,150.395,140.288,131.033,121.58,113.112,104.938,97.3078,90.2178,83.5974,77.2645,70.7126,65.4424,60.1404,55.5644,50.8314,46.3761,43.024,38.625,35.3435,32.2304,29.4192,26.821,24.3303,21.9332,19.4215,16.7163,14.9414,13.1092,0.};
1380 
1381  Float_t fSigmaTPC[100] = {122.209,107.901,103.452,100.498,97.7403,94.7845,93.2543,90.0548,88.1106,85.7382,84.0812,82.2978,80.3817,78.6002,77.3448,75.5086,73.6842,71.9733,70.3447,69.1999,67.878,66.3511,65.0406,63.4866,62.4409,60.7899,59.1328,58.426,56.8618,55.8871,54.1031,53.4959,52.0482,51.0441,49.6218,48.7646,47.5166,46.5247,45.0727,44.4311,43.4531,42.0404,41.0238,40.1384,39.2588,38.2461,36.5951,36.0552,35.3727,33.7883,32.7167,32.4486,31.3709,30.3444,29.505,28.5139,27.4471,26.5359,25.9506,25.127,24.3797,23.2985,22.279,21.4698,20.781,20.8193,19.9509,18.8036,17.9145,16.961,16.7375,15.852,14.9324,14.7663,13.5969,13.4533,12.3067,12.7835,11.7283,10.6758,10.6676,10.6492,9.04614,8.89065,8.66093,8.50997,7.98812,6.91087,7.12045,7.29593,0.};
1382 
1383  for(int i=0;i<90;i++) {
1384  hCentvsTPCmultCuts->SetBinContent(i+1,1,fMeanTPC[i]);
1385  hCentvsTPCmultCuts->SetBinContent(i+1,2,fSigmaTPC[i]);
1386  }
1387 }
1388 
1389 
1390 
1391 
1392 
1393 
1395 
1396  fHistTaskConfigParameters = new TH1F("fHistTaskConfigParameters","Task parameters",20,0,20);
1397  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(1,"FilterBit");
1398  fHistTaskConfigParameters->SetBinContent(1,fFilterBit);
1399  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(2,"n#sigmaTPC");
1400  fHistTaskConfigParameters->SetBinContent(2,fNSigmaCut);
1401  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(3,"MinPt");
1402  fHistTaskConfigParameters->SetBinContent(3,fMinPtCut);
1403  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(4,"MaxPt");
1404  fHistTaskConfigParameters->SetBinContent(4,fMaxPtCut);
1405  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(5,"MinEta");
1406  fHistTaskConfigParameters->SetBinContent(5,fMinEtaCut);
1407  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(6,"MaxEta");
1408  fHistTaskConfigParameters->SetBinContent(6,fMaxEtaCut);
1409 
1410  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(11,"CentralityMin");
1412  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(12,"CentralityMax");
1414 
1415  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(13,"VertexMin(cm)");
1416  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(14,"VertexMax(cm)");
1417 
1419 
1420 
1421 
1422  fHistPileUpCount = new TH1F("fHistPileUpCount", "fHistPileUpCount", 15, 0., 15.);
1423  fHistPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
1424  fHistPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
1425  fHistPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultComb08");
1426  fHistPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
1427  fHistPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>5.0");
1428  fHistPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
1429  fHistPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
1430  Int_t puConst = fPileUpConstParm;
1431  fHistPileUpCount->GetXaxis()->SetBinLabel(8,Form("multESDTPCDif>%d",puConst));
1432  fHistPileUpCount->GetXaxis()->SetBinLabel(9,Form("multGlobTPCDif>%d",puConst));
1433  fHistPileUpCount->GetXaxis()->SetBinLabel(10,"PileUpMultSelTask");
1435 
1436 
1437 
1438 
1439  fHistEventCount = new TH1F("fHistEventCount","Event counts",15,0,15);
1440  fHistEventCount->GetXaxis()->SetBinLabel(1,"Called UserExec()");
1441  fHistEventCount->GetXaxis()->SetBinLabel(2,"Called Exec()");
1442  fHistEventCount->GetXaxis()->SetBinLabel(3,"AOD Exist");
1443  fHistEventCount->GetXaxis()->SetBinLabel(4,"Vz<10");
1444  fHistEventCount->GetXaxis()->SetBinLabel(5,Form("%2.0f<Cent<%2.0f",fCentralityPercentMin,fCentralityPercentMax));
1445  fHistEventCount->GetXaxis()->SetBinLabel(6,"noAODtrack>2 ");
1446  fHistEventCount->GetXaxis()->SetBinLabel(7,"TPC vs Global");
1447  fHistEventCount->GetXaxis()->SetBinLabel(8,"TPC128 vs ESD");
1448  fHistEventCount->GetXaxis()->SetBinLabel(9,"Cent vs TPC");
1449  fHistEventCount->GetXaxis()->SetBinLabel(10,"Survived Events");
1450  fListHist->Add(fHistEventCount);
1451 
1452  //fHistEventCount->Fill(1);
1453 
1454 }
1455 
1456 
1457 
1458 
1459 
1460 
1461 
1462 
1463 
1464 
1466 {
1467  if(fListNUACorr){
1468  for(int i=0;i<5;i++){
1469  fHCorrectNUApos[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Pos_Cent%d_Run%d",i,run));
1470  fHCorrectNUAneg[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Neg_Cent%d_Run%d",i,run));
1471  }
1472  }
1473  else {
1474  printf("\n\n ******** could not open NUA Histograms for run %d, Use Wgt = 1.0 *********\n\n",run);
1475  for(int i=0;i<5;i++){
1476  fHCorrectNUApos[i] = new TH3D(Form("fHCorrectNUApos_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
1477  fHCorrectNUAneg[i] = new TH3D(Form("fHCorrectNUAneg_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
1478  fHCorrectNUApos[i]->SetBinContent(1,1,1,1.0);
1479  fHCorrectNUAneg[i]->SetBinContent(1,1,1,1.0);
1480  //exit(1);
1481  }
1482  }
1483 }
1484 
Int_t charge
TProfile * fHist_Corr3p_EP_Norm_NN[2][3]
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
AliAnalysisUtils * fAnalysisUtil
Double_t mass
AliMultSelection * fMultSelection
PID response Handler.
centrality
char Char_t
Definition: External.C:18
TH1F * fHistPtwithTPCNsigma[3]
last in the list
TCanvas * c
Definition: TestFitELoss.C:172
TH2F * fHistEtaPtAfter
Eta-Pt acceptance.
AliPIDResponse * fPIDResponse
aod
Int_t GetCentralityScaled0to10(Float_t fCent)
int Int_t
Definition: External.C:63
TProfile * fHist_Corr3p_EP_Norm_PN[2][3]
5 centrality bin, read NUA from file
float Float_t
Definition: External.C:68
TProfile * fHist_Reso2n_EP_Norm_Det[2][3]
Definition: External.C:252
Definition: External.C:212
TH3F * fHist3DEtaPhiVz_Neg_Run[3][5]
3 particle 5 centrality bin
TH2F * fHistTPCvsGlobalMultBefore
Eta-Pt acceptance.
TList * fListHist
Event selection.
virtual void UserExec(Option_t *)
TH1F * fHistPileUpCount
Task input parameters FB / cut values etc.
Bool_t PileUpMultiVertex(const AliAODEvent *faod)
TH3D * fHCorrectNUAneg[5]
5 centrality bin, read NUA from file
double GetWDist(const AliVVertex *v0, const AliVVertex *v1)
TH1D * fFB_Efficiency_Cent[10]
3 particle 5 centrality bin
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void SetupMCcorrectionMap(TString sMCfilePath)
TProfile * fHist_Corr3p_EP_Norm_PP[2][3]