AliPhysics  64bcec2 (64bcec2)
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 // #include "vector" not supported by ROOT6 yet.
31 
32 //---- manager and handler---
33 #include "AliAnalysisManager.h"
34 #include "AliInputEventHandler.h"
35 
36 //---V0 and ZDC info---
37 #include "AliAODZDC.h"
38 #include "AliAODVZERO.h"
39 #include "AliAODVertex.h"
40 
41 //---AOD,ESD event--
42 #include "AliESDEvent.h"
43 #include "AliAODHeader.h"
44 #include "AliAODEvent.h"
45 #include "AliAODTrack.h"
46 
47 //----for PID-----
48 #include "AliPIDResponse.h"
49 #include "AliPIDCombined.h"
50 
51 //----- Vevent and tracks
52 #include "AliVEventHandler.h"
53 #include "AliVEvent.h"
54 #include "AliVTrack.h"
55 #include "AliVParticle.h"
56 
57 //----- must include-------
58 #include "AliMultSelection.h"
59 #include "AliAnalysisUtils.h"
60 #include "AliPhysicsSelection.h"
61 #include "AliFlowEventSimple.h"
62 #include "AliAnalysisTaskSE.h"
64 
65 using namespace std;
66 
67 //using std::cout;
68 //using std::endl;
69 
70 
72 
73 
75  fVevent(NULL),
76  fESD(NULL),
77  fAOD(NULL),
78  fPIDResponse(NULL),
79  fMultSelection(NULL),
80  fAnalysisUtil(NULL),
81  fListHist(NULL),
82  mfileFBHijing(NULL),
83  fListFBHijing(NULL),
84  fListNUACorr(NULL),
85  fListV0MCorr(NULL),
86  fHistTaskConfigParameters(NULL),
87  fHistPileUpCount(NULL),
88  fHistMultSelPUCount(NULL),
89  fHistEtaPtBefore(NULL),
90  fHistEtaPtAfter(NULL),
91  fHistTPCvsGlobalMultBefore(NULL),
92  fHistTPCvsGlobalMultAfter(NULL),
93  fHistTPCdEdxvsPBefore(NULL),
94  fHistTPCdEdxvsPAfter(NULL),
95  fHistTOFBetavsPBefore(NULL),
96  fHistTOFBetavsPAfter(NULL),
97  fHistTOFMassvsPtBefore(NULL),
98  fHistTOFMatchCount(NULL),
99  fHistTPCVsESDTrkBefore(NULL),
100  fHistTPCVsESDTrkAfter(NULL),
101  fHistTPConlyVsCL1Before(NULL),
102  fHistTPConlyVsV0MBefore(NULL),
103  fHistTPConlyVsCL1After(NULL),
104  fHistTPConlyVsV0MAfter(NULL),
105  fHistGlobalVsV0MBefore(NULL),
106  fHistGlobalVsV0MAfter(NULL),
107  fHistRawVsCorrMultFB(NULL),
108  hCentvsTPCmultCuts(NULL),
109  fHV0AEventPlaneVsCent(NULL),
110  fHV0CEventPlaneVsCent(NULL),
111  fHTPCAEventPlaneVsCent(NULL),
112  fHTPCCEventPlaneVsCent(NULL),
113  fV0MultChVsRun(NULL),
114  fCentDistBefore(NULL),
115  fCentDistAfter(NULL),
116  fHCorrectV0M(NULL),
117  fHAvgerageQnV0A(NULL),
118  fHAvgerageQnV0C(NULL),
119  fV0AQ2xVsCentRun(NULL),
120  fV0AQ2yVsCentRun(NULL),
121  fV0CQ2xVsCentRun(NULL),
122  fV0CQ2yVsCentRun(NULL),
123  fV0AQ3xVsCentRun(NULL),
124  fV0AQ3yVsCentRun(NULL),
125  fV0CQ3xVsCentRun(NULL),
126  fV0CQ3yVsCentRun(NULL),
127  fTPCAQ2xVsCentRun(NULL),
128  fTPCAQ2yVsCentRun(NULL),
129  fTPCCQ2xVsCentRun(NULL),
130  fTPCCQ2yVsCentRun(NULL),
131  fTPCAQ3xVsCentRun(NULL),
132  fTPCAQ3yVsCentRun(NULL),
133  fTPCCQ3xVsCentRun(NULL),
134  fTPCCQ3yVsCentRun(NULL),
135  fTPCAQ4xVsCentRun(NULL),
136  fTPCAQ4yVsCentRun(NULL),
137  fTPCCQ4xVsCentRun(NULL),
138  fTPCCQ4yVsCentRun(NULL),
139  fFilterBit(1),
140  gN(1),
141  gM(1),
142  gPsiN(2),
143  fOldRunNum(111),
144  fEventCount(0),
145  fNSigmaCut(2),
146  fMinPtCut(0.2),
147  fMaxPtCut(5.0),
148  fMinEtaCut(-0.8),
149  fMaxEtaCut(0.8),
150  fCentralityPercentMin(0),
151  fCentralityPercentMax(90),
152  fPileUpSlopeParm(3.43),
153  fPileUpConstParm(43),
154  bApplyMCcorr(kFALSE),
155  bV0MGainCorr(kFALSE),
156  bSkipPileUpCut(kFALSE),
157  bFillNUAHistPID(kFALSE),
158  sPathOfMCFile("/alien"),
159  sNucleiTP("PbPb"),
160  fHistEventCount(NULL)
161 {
162  for(int i=0;i<3;i++){
163  fHistPtwithTPCNsigma[i]=NULL;
164  fHistPtwithTOFmasscut[i]=NULL;
165  fHistPtwithTOFSignal[i]=NULL;
166  fHistTOFnSigmavsPtAfter[i]=NULL;
167  fHistTPCnSigmavsPtAfter[i]=NULL;
168  fHistTPCTOFnSigmavsPtAfter[i]=NULL;
169  fHistTPCdEdxvsPtPIDAfter[i]=NULL;
170  }
171  for(int i=0;i<5;i++){
172  fHCorrectNUApos[i] = NULL;
173  fHCorrectNUAneg[i] = NULL;
174  }
175  for(int i=0;i<5;i++){ // for PID
176  fHCorrectNUAposPion[i] = NULL;
177  fHCorrectNUAnegPion[i] = NULL;
178  fHCorrectNUAposKaon[i] = NULL;
179  fHCorrectNUAnegKaon[i] = NULL;
180  fHCorrectNUAposProton[i] = NULL;
181  fHCorrectNUAnegProton[i] = NULL;
182  }
183 
184  for(int i=0;i<2;i++){
185  for(int j=0;j<4;j++){
186  fHist_Corr3p_EP_Norm_PN[i][j] = NULL;
187  fHist_Corr3p_EP_Norm_PP[i][j] = NULL;
188  fHist_Corr3p_EP_Norm_NN[i][j] = NULL;
189  }
190  for(int j=0;j<4;j++) {
191  fHist_Reso2n_EP_Norm_Det[i][j] = NULL;
192  }
193  //PID:
194  for(int j=0;j<4;j++){
195  fHist_Corr3p_Pion_EP_Norm_PN[i][j] = NULL;
196  fHist_Corr3p_Pion_EP_Norm_PP[i][j] = NULL;
197  fHist_Corr3p_Pion_EP_Norm_NN[i][j] = NULL;
198  fHist_Corr3p_Kaon_EP_Norm_PN[i][j] = NULL;
199  fHist_Corr3p_Kaon_EP_Norm_PP[i][j] = NULL;
200  fHist_Corr3p_Kaon_EP_Norm_NN[i][j] = NULL;
201  fHist_Corr3p_Proton_EP_Norm_PN[i][j] = NULL;
202  fHist_Corr3p_Proton_EP_Norm_PP[i][j] = NULL;
203  fHist_Corr3p_Proton_EP_Norm_NN[i][j] = NULL;
204  }
205  }
206  for(int i=0;i<4;i++){
207  for(int j=0;j<5;j++){
208  fHist3DEtaPhiVz_Pos_Run[i][j]=NULL;
209  fHist3DEtaPhiVz_Neg_Run[i][j]=NULL;
210  }
211  }
212  for(int i=0;i<10;i++){
213  fFB_Efficiency_Cent[i] = NULL;
214  }
215 
216 
217  DefineInput(0,TChain::Class());
218  DefineOutput(1,TList::Class());
219 }
220 //______________________________empty constructor_______________________
223  fVevent(NULL),
224  fESD(NULL),
225  fAOD(NULL),
226  fPIDResponse(NULL),
227  fMultSelection(NULL),
228  fAnalysisUtil(NULL),
229  fListHist(NULL),
230  mfileFBHijing(NULL),
231  fListFBHijing(NULL),
232  fListNUACorr(NULL),
233  fListV0MCorr(NULL),
234  fHistTaskConfigParameters(NULL),
235  fHistPileUpCount(NULL),
236  fHistMultSelPUCount(NULL),
237  fHistEtaPtBefore(NULL),
238  fHistEtaPtAfter(NULL),
239  fHistTPCvsGlobalMultBefore(NULL),
240  fHistTPCvsGlobalMultAfter(NULL),
241  fHistTPCdEdxvsPBefore(NULL),
242  fHistTPCdEdxvsPAfter(NULL),
243  fHistTOFBetavsPBefore(NULL),
244  fHistTOFBetavsPAfter(NULL),
245  fHistTOFMassvsPtBefore(NULL),
246  fHistTOFMatchCount(NULL),
247  fHistTPCVsESDTrkBefore(NULL),
248  fHistTPCVsESDTrkAfter(NULL),
249  fHistTPConlyVsCL1Before(NULL),
250  fHistTPConlyVsV0MBefore(NULL),
251  fHistTPConlyVsCL1After(NULL),
252  fHistTPConlyVsV0MAfter(NULL),
253  fHistGlobalVsV0MBefore(NULL),
254  fHistGlobalVsV0MAfter(NULL),
255  fHistRawVsCorrMultFB(NULL),
256  hCentvsTPCmultCuts(NULL),
257  fHV0AEventPlaneVsCent(NULL),
258  fHV0CEventPlaneVsCent(NULL),
259  fHTPCAEventPlaneVsCent(NULL),
260  fHTPCCEventPlaneVsCent(NULL),
261  fV0MultChVsRun(NULL),
262  fCentDistBefore(NULL),
263  fCentDistAfter(NULL),
264  fHCorrectV0M(NULL),
265  fHAvgerageQnV0A(NULL),
266  fHAvgerageQnV0C(NULL),
267  fV0AQ2xVsCentRun(NULL),
268  fV0AQ2yVsCentRun(NULL),
269  fV0CQ2xVsCentRun(NULL),
270  fV0CQ2yVsCentRun(NULL),
271  fV0AQ3xVsCentRun(NULL),
272  fV0AQ3yVsCentRun(NULL),
273  fV0CQ3xVsCentRun(NULL),
274  fV0CQ3yVsCentRun(NULL),
275  fTPCAQ2xVsCentRun(NULL),
276  fTPCAQ2yVsCentRun(NULL),
277  fTPCCQ2xVsCentRun(NULL),
278  fTPCCQ2yVsCentRun(NULL),
279  fTPCAQ3xVsCentRun(NULL),
280  fTPCAQ3yVsCentRun(NULL),
281  fTPCCQ3xVsCentRun(NULL),
282  fTPCCQ3yVsCentRun(NULL),
283  fTPCAQ4xVsCentRun(NULL),
284  fTPCAQ4yVsCentRun(NULL),
285  fTPCCQ4xVsCentRun(NULL),
286  fTPCCQ4yVsCentRun(NULL),
287  fFilterBit(1),
288  gN(1),
289  gM(1),
290  gPsiN(2),
291  fOldRunNum(111),
292  fEventCount(0),
293  fNSigmaCut(2),
294  fMinPtCut(0.2),
295  fMaxPtCut(5.0),
296  fMinEtaCut(-0.8),
297  fMaxEtaCut(0.8),
298  fCentralityPercentMin(0),
299  fCentralityPercentMax(90),
300  fPileUpSlopeParm(3.43),
301  fPileUpConstParm(43),
302  bApplyMCcorr(kFALSE),
303  bV0MGainCorr(kFALSE),
304  bSkipPileUpCut(kFALSE),
305  bFillNUAHistPID(kFALSE),
306  sPathOfMCFile("/alien"),
307  sNucleiTP("PbPb"),
308  fHistEventCount(NULL)
309 {
310  for(int i=0;i<3;i++){
311  fHistPtwithTPCNsigma[i]=NULL;
312  fHistPtwithTOFmasscut[i]=NULL;
313  fHistPtwithTOFSignal[i]=NULL;
314  fHistTOFnSigmavsPtAfter[i]=NULL;
315  fHistTPCnSigmavsPtAfter[i]=NULL;
317  fHistTPCdEdxvsPtPIDAfter[i]=NULL;
318  }
319  for(int i=0;i<5;i++){
320  fHCorrectNUApos[i] = NULL;
321  fHCorrectNUAneg[i] = NULL;
322  }
323  for(int i=0;i<5;i++){ // for PID
324  fHCorrectNUAposPion[i] = NULL;
325  fHCorrectNUAnegPion[i] = NULL;
326  fHCorrectNUAposKaon[i] = NULL;
327  fHCorrectNUAnegKaon[i] = NULL;
328  fHCorrectNUAposProton[i] = NULL;
329  fHCorrectNUAnegProton[i] = NULL;
330  }
331 
332  for(int i=0;i<2;i++){
333  for(int j=0;j<4;j++){
334  fHist_Corr3p_EP_Norm_PN[i][j] = NULL;
335  fHist_Corr3p_EP_Norm_PP[i][j] = NULL;
336  fHist_Corr3p_EP_Norm_NN[i][j] = NULL;
337  }
338  for(int j=0;j<4;j++) {
339  fHist_Reso2n_EP_Norm_Det[i][j] = NULL;
340  }
341  //PID:
342  for(int j=0;j<4;j++){
343  fHist_Corr3p_Pion_EP_Norm_PN[i][j] = NULL;
344  fHist_Corr3p_Pion_EP_Norm_PP[i][j] = NULL;
345  fHist_Corr3p_Pion_EP_Norm_NN[i][j] = NULL;
346  fHist_Corr3p_Kaon_EP_Norm_PN[i][j] = NULL;
347  fHist_Corr3p_Kaon_EP_Norm_PP[i][j] = NULL;
348  fHist_Corr3p_Kaon_EP_Norm_NN[i][j] = NULL;
349  fHist_Corr3p_Proton_EP_Norm_PN[i][j] = NULL;
350  fHist_Corr3p_Proton_EP_Norm_PP[i][j] = NULL;
351  fHist_Corr3p_Proton_EP_Norm_NN[i][j] = NULL;
352  }
353  }
354  for(int i=0;i<4;i++){
355  for(int j=0;j<5;j++){
356  fHist3DEtaPhiVz_Pos_Run[i][j]=NULL;
357  fHist3DEtaPhiVz_Neg_Run[i][j]=NULL;
358  }
359  }
360  for(int i=0;i<10;i++){
361  fFB_Efficiency_Cent[i] = NULL;
362  }
363 }
364 
365 //___________________________ destructor ___________________________
367 {
368  //Destructor
369  //if(fPIDResponse) delete fPIDResponse;
370  //if(fMultSelection) delete fMultSelection;
371 
372  if(fListHist){
373  delete fListHist;
374  }
375 
376  if(fHCorrectV0M) delete fHCorrectV0M;
377  if(fHAvgerageQnV0A) delete fHAvgerageQnV0A;
378  if(fHAvgerageQnV0C) delete fHAvgerageQnV0C;
379 
380  if(mfileFBHijing->IsOpen()){
381  mfileFBHijing->Close();
382  if(fListFBHijing) delete fListFBHijing;
383  }
384  for(int i=0;i<10;i++){
385  if(fFB_Efficiency_Cent[i])
386  delete fFB_Efficiency_Cent[i];
387  }
388  for(int i=0;i<5;i++){
389  if(fHCorrectNUApos[i]) delete fHCorrectNUApos[i];
390  if(fHCorrectNUAneg[i]) delete fHCorrectNUAneg[i];
391  }
392  for(int i=0;i<5;i++){ // for PID
393  if(fHCorrectNUAposPion[i]) delete fHCorrectNUAposPion[i];
394  if(fHCorrectNUAnegPion[i]) delete fHCorrectNUAnegPion[i];
395  if(fHCorrectNUAposKaon[i]) delete fHCorrectNUAposKaon[i];
396  if(fHCorrectNUAnegKaon[i]) delete fHCorrectNUAnegKaon[i];
399  }
400 
401 
402  if(fAnalysisUtil) delete fAnalysisUtil; // its 'new' !!
403 
404 }//---------------- sanity ------------------------
405 
406 
407 
409 {
410 
411  //std::cout<<"\n..UserCreateOutputObject called.. with isCorr = "<<isCorr<<"\n ....check if succeeded...\n"<<endl;
412 
413  //input hander
414  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
415  AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
416  if (!inputHandler) { printf("\n ...Input handler missing!!!...\n"); return; }
417 
418  //PileUp Multi-Vertex
419  fAnalysisUtil = new AliAnalysisUtils();
420  fAnalysisUtil->SetUseMVPlpSelection(kTRUE);
421  fAnalysisUtil->SetUseOutOfBunchPileUp(kTRUE);
422 
423  //pid response object
424  fPIDResponse=inputHandler->GetPIDResponse();
425 
426  fListHist = new TList();
427  fListHist->SetOwner(kTRUE);
428 
429 
431 
432  if(!gGrid){
433  TGrid::Connect("alien://");
434  }
435 
437 
438 
439 
440  hCentvsTPCmultCuts = new TH2F("hCentvsTPCmultCuts","TPCmult Low,high",100,0,100,5,0,5);
441  hCentvsTPCmultCuts->GetYaxis()->SetBinLabel(1,"mean");
442  hCentvsTPCmultCuts->GetYaxis()->SetBinLabel(2,"sigma");
444 
446 
447  fHistTOFMatchCount = new TH2F("fHistTOFMatchCount","TofMatchFlag vs mismatch Prob",10,0,10,200,-5,5);
449 
450  fHistEtaPtBefore = new TH2F("fHistEtaPtBefore","Eta vs pT",100,-1.25,1.25,100,0,10);
452 
453  fHistEtaPtAfter = new TH2F("fHistEtaPtAfter","Eta vs pT",100,-1.25,1.25,100,0,10);
455 
456 
457  Int_t gMaxTPCFB1mult = 0;
458  Int_t gMaxGlobalmult = 0;
459  Int_t gMaxTPCcorrmult = 0;
460  Int_t gMaxESDtracks = 0;
461 
462  if(sNucleiTP=="pp"||sNucleiTP=="PP"){
463  gMaxGlobalmult = 200;
464  gMaxTPCFB1mult = 200;
465  gMaxTPCcorrmult = 500;
466  gMaxESDtracks = 1000;
467  //fSkipOutlierCut = 1;
468  }
469  else if(sNucleiTP=="pPb"||sNucleiTP=="Pbp"||sNucleiTP=="PbP"||sNucleiTP=="PPb"){
470  gMaxGlobalmult = 400;
471  gMaxTPCFB1mult = 400;
472  gMaxTPCcorrmult = 500;
473  gMaxESDtracks = 2000;
474  //fSkipOutlierCut = 1;
475  }
476  else{
477  gMaxGlobalmult = 4000;
478  gMaxTPCFB1mult = 4000;
479  gMaxTPCcorrmult = 5000;
480  gMaxESDtracks = 20000;
481  //fSkipOutlierCut = 0;
482  }
483 
484  //if(bSkipPileUpCut) { fSkipOutlierCut = 1;}
485 
486 
487  fHistTPCVsESDTrkBefore = new TH2F("fHistTPCVsESDTrkBefore","Before; TPC1; ESD trk",500,0,gMaxTPCcorrmult,200,0,gMaxESDtracks);
489  fHistTPCVsESDTrkAfter = new TH2F("fHistTPCVsESDTrkAfter"," After; TPC1; ESD trk",500,0,gMaxTPCcorrmult,200,0,gMaxESDtracks);
491 
492  fHistTPCvsGlobalMultBefore = new TH2F("fHistTPCvsGlobalMultBefore","Before; Global; TPC(fb1) ",200,0,gMaxGlobalmult,200,0,gMaxTPCFB1mult);
494  fHistTPCvsGlobalMultAfter = new TH2F("fHistTPCvsGlobalMultAfter"," After; Global; TPC(fb1) ",200,0,gMaxGlobalmult,200,0,gMaxTPCFB1mult);
496 
497  fHistGlobalVsV0MBefore = new TH2F("fHistGlobalVsV0MBefore","Before;Cent(V0M);Global",100,0,100,500,0,gMaxGlobalmult);
499  fHistGlobalVsV0MAfter = new TH2F("fHistGlobalVsV0MAfter"," After; Cent(V0M);Global",100,0,100,500,0,gMaxGlobalmult);
501 
502  fHistTPConlyVsCL1Before = new TH2F("fHistTPConlyVsCL1Before","Before;Cent(CL1); TPC(FB1)",100,0,100,250,0,gMaxTPCcorrmult);
504  fHistTPConlyVsCL1After = new TH2F("fHistTPConlyVsCL1After","After; Cent(CL1); TPC(FB1) ",100,0,100,250,0,gMaxTPCcorrmult);
506 
507  fHistTPConlyVsV0MBefore = new TH2F("fHistTPConlyVsV0MBefore","Before;Cent(V0M); TPC(FB1)",100,0,100,250,0,gMaxTPCcorrmult);
509  fHistTPConlyVsV0MAfter = new TH2F("fHistTPConlyVsV0MAfter","After; Cent(V0M); TPC(FB1) ",100,0,100,250,0,gMaxTPCcorrmult);
511 
512  fHistRawVsCorrMultFB = new TH2F("fHistRawVsCorrMultFB",Form("FB%d;Mult_{raw};Mult_{corr}",fFilterBit),gMaxTPCFB1mult,0,gMaxTPCFB1mult,gMaxTPCcorrmult,0,gMaxTPCcorrmult);
514 
515 
516  fHistTPCdEdxvsPBefore = new TH2F("fHistTPCdEdxvsPBefore","Before; p (GeV/c); dEdx (arb)",200,-5,5,200,0,250);
518  fHistTPCdEdxvsPAfter = new TH2F("fHistTPCdEdxvsPAfter"," After; p (GeV/c); dEdx (arb)",200,-5,5, 200,0,250);
520 
521 
522  fHistTOFBetavsPBefore = new TH2F("fHistTOFBetavsPBefore","Before; p (GeV/c); beta ",200,-5,5,100,0.0,1.2);
524  fHistTOFBetavsPAfter = new TH2F("fHistTOFBetavsPAfter"," After; p (GeV/c); beta ",200,-5,5,100,0.0,1.2);
526 
527 
528  fHistTOFMassvsPtBefore = new TH2F("fHistTOFMassvsPtBefore","Before; p_{T}(GeV/c); m^{2}(GeV^{2}/c^{4})",200,-5,5,500,-0.5,4.5);
530 
531  fCentDistBefore = new TH1F("fCentDistBefore","no Cut; Cent (%); Events ",100,0,100);
533 
534  fCentDistAfter = new TH1F("fCentDistAfter","with Cut; Cent (%); Events ",100,0,100);
536 
537 
538  //---------------- PID Histograms ---------------------
539  //Dont forget to add histograms in the List. !!!
540 
541 
542  const char *gSpecies[4] = {"Pion","Kaon","proton","Charge"};
543 
544  for(int i=0;i<3;i++){
545  fHistPtwithTPCNsigma[i] = new TH1F(Form("fHistPtwithTPCNsigma_%s",gSpecies[i]), Form("%s;p_{T}(GeV/c))",gSpecies[i]),200,-5,5);
547  fHistPtwithTOFmasscut[i] = new TH1F(Form("fHistPtwithTOFmasscut_%s",gSpecies[i]),Form("%s;p_{T}(GeV/c))",gSpecies[i]),200,-5,5);
549  fHistPtwithTOFSignal[i] = new TH1F(Form("fHistPtwithTOFSignal_%s", gSpecies[i]),Form("%s;p_{T}(GeV/c))",gSpecies[i]),200,-5,5);
551 
552  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);
554 
555  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);
557 
558  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);
560 
561  fHistTPCdEdxvsPtPIDAfter[i] = new TH2F(Form("fHistTPCdEdxvsPtAfter_%s",gSpecies[i]),"AfterCut; p_{T} (GeV/c); dEdx (arb)",400,0,10,200,0,250);
563  }// PID histograms done
564 
565 
566 
567  Double_t centRange[11] = {0,5,10,20,30,40,50,60,70,80,90};
568  const char *gDetForEP[4] = {"V0A","V0C","TPC-A","TPC-C"};
569 
570  //------------------- CME 3p correlator Charged hadrons (EP method) ------------------
571  for(int i=0;i<2;i++){
572  for(int j=0;j<4;j++){
573  //Detector: 0 = V0A, 1 = V0C, 3 = TPCA, 4 = TPCC
574  fHist_Corr3p_EP_Norm_PN[i][j] = new TProfile(Form("fHist_Corr3p_EP_Norm_PosNeg_Mag%d_Det%d",i,j+1),Form("US, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
575  //fHist_Corr3p_EP_Norm_PN[i][j]->Sumw2();
577  fHist_Corr3p_EP_Norm_PP[i][j] = new TProfile(Form("fHist_Corr3p_EP_Norm_PosPos_Mag%d_Det%d",i,j+1),Form("P-P, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
578  //fHist_Corr3p_EP_Norm_PP[i][j]->Sumw2();
580  fHist_Corr3p_EP_Norm_NN[i][j] = new TProfile(Form("fHist_Corr3p_EP_Norm_NegNeg_Mag%d_Det%d",i,j+1),Form("N-N, #Psi_{2}, %s",gDetForEP[j]),10,centRange,"");
581  //fHist_Corr3p_EP_Norm_NN[i][j]->Sumw2();
583  }
584  //EP Resolution:
585  for(int j=0;j<4;j++){
586  //Det: 0 = v0c-v0a, 1 = v0a-TPC, 2 = v0c-TPC, 3 =TPC-A TPC-C
587  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,"");
588  //fHist_Reso2n_EP_Norm_Det[i][j]->Sumw2();
590  }
591  //----------- PID -------------------
592  for(int j=0;j<4;j++){ //Detector: 0 = V0A, 1 = V0C, 3 = TPCA, 4 = TPCC
593  //----------> Pion:
594  fHist_Corr3p_Pion_EP_Norm_PN[i][j] = new TProfile(Form("fHist_Corr3p_Pion_EP_Norm_PosNeg_Mag%d_Det%d",i,j+1),Form("US, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
595  //fHist_Corr3p_Pion_EP_Norm_PN[i][j]->Sumw2();
597  fHist_Corr3p_Pion_EP_Norm_PP[i][j] = new TProfile(Form("fHist_Corr3p_Pion_EP_Norm_PosPos_Mag%d_Det%d",i,j+1),Form("P-P, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
598  //fHist_Corr3p_Pion_EP_Norm_PP[i][j]->Sumw2();
600  fHist_Corr3p_Pion_EP_Norm_NN[i][j] = new TProfile(Form("fHist_Corr3p_Pion_EP_Norm_NegNeg_Mag%d_Det%d",i,j+1),Form("N-N, #Psi_{2}, %s",gDetForEP[j]),10,centRange,"");
601  //fHist_Corr3p_Pion_EP_Norm_NN[i][j]->Sumw2();
603  //----------> Kaon:
604  fHist_Corr3p_Kaon_EP_Norm_PN[i][j] = new TProfile(Form("fHist_Corr3p_Kaon_EP_Norm_PosNeg_Mag%d_Det%d",i,j+1),Form("US, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
605  //fHist_Corr3p_Kaon_EP_Norm_PN[i][j]->Sumw2();
607  fHist_Corr3p_Kaon_EP_Norm_PP[i][j] = new TProfile(Form("fHist_Corr3p_Kaon_EP_Norm_PosPos_Mag%d_Det%d",i,j+1),Form("P-P, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
608  //fHist_Corr3p_Kaon_EP_Norm_PP[i][j]->Sumw2();
610  fHist_Corr3p_Kaon_EP_Norm_NN[i][j] = new TProfile(Form("fHist_Corr3p_Kaon_EP_Norm_NegNeg_Mag%d_Det%d",i,j+1),Form("N-N, #Psi_{2}, %s",gDetForEP[j]),10,centRange,"");
611  //fHist_Corr3p_Kaon_EP_Norm_NN[i][j]->Sumw2();
613  //----------> Proton:
614  fHist_Corr3p_Proton_EP_Norm_PN[i][j] = new TProfile(Form("fHist_Corr3p_Proton_EP_Norm_PosNeg_Mag%d_Det%d",i,j+1),Form("US, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
615  //fHist_Corr3p_Proton_EP_Norm_PN[i][j]->Sumw2();
617  fHist_Corr3p_Proton_EP_Norm_PP[i][j] = new TProfile(Form("fHist_Corr3p_Proton_EP_Norm_PosPos_Mag%d_Det%d",i,j+1),Form("P-P, #Psi_{2} %s",gDetForEP[j]),10,centRange,"");
618  //fHist_Corr3p_Proton_EP_Norm_PP[i][j]->Sumw2();
620  fHist_Corr3p_Proton_EP_Norm_NN[i][j] = new TProfile(Form("fHist_Corr3p_Proton_EP_Norm_NegNeg_Mag%d_Det%d",i,j+1),Form("N-N, #Psi_{2}, %s",gDetForEP[j]),10,centRange,"");
621  //fHist_Corr3p_Proton_EP_Norm_NN[i][j]->Sumw2();
623  }//Det loop
624  }//magfield loop
625 
626 
627  Double_t truncPi = 3.1416;
628 
629  fHV0AEventPlaneVsCent = new TH2F("fHV0AEventPlaneVsCent",Form("Psi %d from V0A", gPsiN), 10,centRange,50,-0.0,truncPi);
631 
632  fHV0CEventPlaneVsCent = new TH2F("fHV0CEventPlaneVsCent",Form("Psi %d from V0C", gPsiN), 10,centRange,50,-0.0,truncPi);
634 
635  fHTPCAEventPlaneVsCent = new TH2F("fHTPCAEventPlaneVsCent",Form("Psi %d from V0A",gPsiN),10,centRange,50,-0.0,truncPi);
637 
638  fHTPCCEventPlaneVsCent = new TH2F("fHTPCCEventPlaneVsCent",Form("Psi %d from V0A",gPsiN),10,centRange,50,-0.0,truncPi);
640  //-------------------------------------------------------------------------------
641 
642 
643  //---- to store NUA and calib histograms -----
644  TList *fListNUACalib = new TList();
645  fListNUACalib->SetName("fListNUACalib");
646  fListNUACalib->SetOwner(kTRUE);
647 
648 
649  //----------------- V0 Calibration hist: ---------------------
650  fV0MultChVsRun = new TH2F("fV0MultChVsRun","1-32 V0C, 33-64 V0A",64,0,64,20,0,100);
651  fListNUACalib->Add(fV0MultChVsRun);
652 
653  const char *sCorrect[2]={"w Corr","w/o Corr"};
654  Int_t isCorr = 1;
655 
656  if(fListV0MCorr){
657  isCorr = 0;
658  }
659 
660 
661  fV0AQ2xVsCentRun = new TProfile("fV0ACos2nVsCentRun",Form("<Cos2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
662  fListNUACalib->Add(fV0AQ2xVsCentRun);
663  fV0AQ2yVsCentRun = new TProfile("fV0ASin2nVsCentRun",Form("<Sin2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
664  fListNUACalib->Add(fV0AQ2yVsCentRun);
665  fV0CQ2xVsCentRun = new TProfile("fV0CCos2nVsCentRun",Form("<Cos2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
666  fListNUACalib->Add(fV0CQ2xVsCentRun);
667  fV0CQ2yVsCentRun = new TProfile("fV0CSin2nVsCentRun",Form("<Sin2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
668  fListNUACalib->Add(fV0CQ2yVsCentRun);
669 
670  fV0AQ3xVsCentRun = new TProfile("fV0ACos3nVsCentRun",Form("<Cos3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
671  fListNUACalib->Add(fV0AQ3xVsCentRun);
672  fV0AQ3yVsCentRun = new TProfile("fV0ASin3nVsCentRun",Form("<Sin3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
673  fListNUACalib->Add(fV0AQ3yVsCentRun);
674  fV0CQ3xVsCentRun = new TProfile("fV0CCos3nVsCentRun",Form("<Cos3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
675  fListNUACalib->Add(fV0CQ3xVsCentRun);
676  fV0CQ3yVsCentRun = new TProfile("fV0CSin3nVsCentRun",Form("<Sin3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
677  fListNUACalib->Add(fV0CQ3yVsCentRun);
678 
679  isCorr = 1;
680  if(fListNUACorr){
681  isCorr = 0;
682  }
683  //------------------- TPC Qvector Recentering Histograms --------------
684  fTPCAQ2xVsCentRun = new TProfile("fTPCACos2nVsCentRun",Form("<Cos2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
685  fListNUACalib->Add(fTPCAQ2xVsCentRun);
686  fTPCAQ2yVsCentRun = new TProfile("fTPCASin2nVsCentRun",Form("<Sin2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
687  fListNUACalib->Add(fTPCAQ2yVsCentRun);
688  fTPCCQ2xVsCentRun = new TProfile("fTPCCCos2nVsCentRun",Form("<Cos2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
689  fListNUACalib->Add(fTPCCQ2xVsCentRun);
690  fTPCCQ2yVsCentRun = new TProfile("fTPCCSin2nVsCentRun",Form("<Sin2> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
691  fListNUACalib->Add(fTPCCQ2yVsCentRun);
692 
693  fTPCAQ3xVsCentRun = new TProfile("fTPCACos3nVsCentRun",Form("<Cos3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
694  fListNUACalib->Add(fTPCAQ3xVsCentRun);
695  fTPCAQ3yVsCentRun = new TProfile("fTPCASin3nVsCentRun",Form("<Sin3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
696  fListNUACalib->Add(fTPCAQ3yVsCentRun);
697  fTPCCQ3xVsCentRun = new TProfile("fTPCCCos3nVsCentRun",Form("<Cos3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
698  fListNUACalib->Add(fTPCCQ3xVsCentRun);
699  fTPCCQ3yVsCentRun = new TProfile("fTPCCSin3nVsCentRun",Form("<Sin3> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
700  fListNUACalib->Add(fTPCCQ3yVsCentRun);
701 
702  fTPCAQ4xVsCentRun = new TProfile("fTPCACos4nVsCentRun",Form("<Cos4> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
703  fListNUACalib->Add(fTPCAQ4xVsCentRun);
704  fTPCAQ4yVsCentRun = new TProfile("fTPCASin4nVsCentRun",Form("<Sin4> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
705  fListNUACalib->Add(fTPCAQ4yVsCentRun);
706  fTPCCQ4xVsCentRun = new TProfile("fTPCCCos4nVsCentRun",Form("<Cos4> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
707  fListNUACalib->Add(fTPCCQ4xVsCentRun);
708  fTPCCQ4yVsCentRun = new TProfile("fTPCCSin4nVsCentRun",Form("<Sin4> vs cent (%s)",sCorrect[isCorr]),90,0,90,"");
709  fListNUACalib->Add(fTPCCQ4yVsCentRun);
710 
711  //-------------------------------------------------------------------------------
712 
713 
714 
715 
716  //-------------------------- Define NUA Hist for PID -----------------------------
717  Int_t gCentForNUA[6] = {0,5,10,20,40,90};
718  Char_t name[100];
719  Char_t title[100];
720 
721  for(int i=0;i<4;i++){
722  for(int j=0;j<5;j++){
723  sprintf(name,"fHistEtaPhiVz_%s_Pos_Cent%d_Run%d",gSpecies[i],j,1);
724  sprintf(title,"eta,phi,Vz %sPos, Cent%d-%d, FB %d",gSpecies[i],gCentForNUA[j],gCentForNUA[j+1],fFilterBit);
725  fHist3DEtaPhiVz_Pos_Run[i][j] = new TH3F(name,title,10,-10,10,50,0,6.283185,16,-0.8,0.8);
726  fListNUACalib->Add(fHist3DEtaPhiVz_Pos_Run[i][j]);
727 
728  sprintf(name,"fHistEtaPhiVz_%s_Neg_Cent%d_Run%d",gSpecies[i],j,1);
729  sprintf(title,"eta,phi,Vz %sNeg, Cent%d-%d, FB %d",gSpecies[i],gCentForNUA[j],gCentForNUA[j+1],fFilterBit);
730  fHist3DEtaPhiVz_Neg_Run[i][j] = new TH3F(name,title,10,-10,10,50,0,6.283185,16,-0.8,0.8);
731  fListNUACalib->Add(fHist3DEtaPhiVz_Neg_Run[i][j]);
732  }
733  }
734  //---------------------------------------------------------------------------------
735 
736 
737  fListHist->Add(fListNUACalib);
738 
739  PostData(1,fListHist);
740  cout<<"\n.........UserCreateOutputObject called.........\n fFilterBit = "<<fFilterBit<<" CentMax = "<<fCentralityPercentMax;
741  cout<<" PU C = "<<fPileUpConstParm<<" PsiN = "<<gPsiN<<"\n\n"<<endl;
742 
743  //Reset the counter:
744  //watch.Reset();
745 
746 }
747 
748 
749 
750 
751 
752 
753 
754 
755 
756 
757 
758 
759 
760 
761 
762 
763 
764 //______________________________________________________________________
766  //printf("info: UserExec is called.... 1 \n");
767 
768  //watch.Start(kTRUE); //debug only
769 
770  if(fEventCount==100) return;
771 
772  Float_t stepCount = 0.5;
773 
774  fHistEventCount->Fill(stepCount); //1
775  stepCount++;
776 
777  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
778  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
779 
780  if(!(fESD || fAOD)){ printf("ERROR: fESD & fAOD not available\n"); return; }
781 
782  fVevent = dynamic_cast<AliVEvent*>(InputEvent());
783 
784  if (!fVevent) { printf("ERROR: fVevent not available\n"); return; }
785 
786  fHistEventCount->Fill(stepCount); //2
787  stepCount++;
788 
789 
790 
791 
792 
793 
794 
795  //--------- Check if I have PID response object --------
796  if(!fPIDResponse){
797  printf("\n\n...... PIDResponse object not found..... \n\n");
798  return;
799  }
800 
801 
802 
803 
804 
805 
806 
807 
808  //-------------- Vtx cuts ---------------
809  const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
810 
811  Double_t pVtxZ = -999;
812  pVtxZ = pVtx->GetZ();
813 
814  if(TMath::Abs(pVtxZ)>10.) return;
815 
816  fHistEventCount->Fill(stepCount); //3
817  stepCount++;
818 
819  fMultSelection = (AliMultSelection*) InputEvent()->FindListObject("MultSelection");
820 
821  if(!fMultSelection) { printf("\n\n **WARNING** \n::UserExec() AliMultSelection object not found.\n\n"); exit(1); }
822 
823  Float_t centrV0M = fMultSelection->GetMultiplicityPercentile("V0M");
824  Float_t centrCL1 = fMultSelection->GetMultiplicityPercentile("CL1");
825 //Float_t centrCL0 = fMultSelection->GetMultiplicityPercentile("CL0");
826 //Float_t centrTRK = fMultSelection->GetMultiplicityPercentile("TRK");
827 
828  Float_t centrality = centrV0M;
829 
831  return;
832  }
833 
834  fHistEventCount->Fill(stepCount); //4
835  stepCount++;
836 
837  fCentDistBefore->Fill(centrality);
838 
839 
840 
841 
842 
843 
844 
845 
846  Int_t ntracks=fAOD->GetNumberOfTracks();
847  if(ntracks<2) return; // Check this cut....!!!
848 
849  fHistEventCount->Fill(stepCount); //5
850  stepCount++;
851 
852 
853 
854 
855 
856 
857  Int_t cent10bin = -1;
858 
859  cent10bin = GetCentralityScaled0to10(centrality); //Centrality in 0-10 scale
860 
861 
862 //Centrality array index for NUA correcion
863  Int_t cForNUA = 0;
864 
865  if(centrality<5.0) {
866  cForNUA = 0;
867  }
868  else if(centrality>=5.0 && centrality<10){
869  cForNUA = 1; // 1=5-10,
870  }
871  else if(centrality>=10.0 && centrality<20) {
872  cForNUA = 2; // 2 = 10-20,
873  }
874  else if(centrality>=20 && centrality<40){
875  cForNUA = 3; // 3=20-40
876  }
877  else if(centrality>=40){
878  cForNUA = 4; // 4=40-90
879  }
880 
881 
882 
883 
884 
885 
886 
887 
888  //---------- Magnetic field --------
889  Double_t fMagField = fAOD->GetMagneticField();
890 
891  const Int_t QAindex = (fMagField > 0) ? 1 : 0;
892  //---------------------------------
893 
894 
895 
896 
897 
898 
899  //Load NUA and V0M correction map run by run:
900  Int_t runNumber = fAOD->GetRunNumber();
901 
902  if(runNumber!=fOldRunNum) {
903 
904  GetNUACorrectionHist(runNumber);
905 
906  if(bV0MGainCorr) {
907  GetV0MCorrectionHist(runNumber);
908  }
909 
910  fOldRunNum = runNumber;
911  }
912  //------------------------------------------
913 
914 
915 
916 
917 
918 
919 
920  //----- Event Plane variables:-------
921  Double_t PsiNV0A = 0;
922  Double_t PsiNV0C = 0;
923 
924  Double_t PsiNTPCA = 0; // eta <0
925  Double_t PsiNTPCC = 0; // eta >0
926  Double_t PsiNTPCF = 0; // Full TPC
927 
928  Double_t sumTPCQn2x[5] = {0,0,0}; //[0]= eta<0; [1]= eta>0; [2]= -0.8 < eta < 0.8
929  Double_t sumTPCQn2y[5] = {0,0,0};
930  Double_t sumTPCQn3x[5] = {0,0,0}; //[0]= eta<0; [1]= eta>0; [2]= -0.8 < eta < 0.8
931  Double_t sumTPCQn3y[5] = {0,0,0};
932  Double_t sumTPCQn4x[5] = {0,0,0}; //[0]= eta<0; [1]= eta>0; [2]= -0.8 < eta < 0.8
933  Double_t sumTPCQn4y[5] = {0,0,0};
934  //------------------------------------
935 
936 
937 
938 
939 
940 
941 
942  //Variables for MC tracking correction
943  Int_t ptBinMC = 1;
944  Int_t iBinNUA = 1;
945  Float_t ptWgtMC = 1.0;
946  Float_t WgtNUA = 1.0;
947  Float_t ptTrk = 0.1;
948  Float_t dEdx = 0;
949 
950 
951  //-------------- Track loop for outlier and PileUp cut -------------------
952 
953  //---------------- a dobrin --------------
954 
955  Bool_t bIsPileup=kFALSE;
956 
957  Int_t isPileup = fAOD->IsPileupFromSPD(3);
958 
959  if(isPileup != 0) {
960  fHistPileUpCount->Fill(0.5);
961  bIsPileup=kTRUE;
962  }
963  else if(PileUpMultiVertex(fAOD)) {
964  fHistPileUpCount->Fill(1.5);
965  bIsPileup=kTRUE;
966  }
967  else if(((AliAODHeader*)fAOD->GetHeader())->GetRefMultiplicityComb08() < 0) {
968  fHistPileUpCount->Fill(2.5);
969  bIsPileup=kTRUE;
970  }
971  else if(fAOD->IsIncompleteDAQ()) {
972  fHistPileUpCount->Fill(3.5);
973  bIsPileup=kTRUE;
974  }
975  else if(fabs(centrV0M-centrCL1)> 5.0) {//default: 7.5
976 //else if(fabs(centrV0M-centrCL1)> 7.5) {//default: 7.5
977  fHistPileUpCount->Fill(4.5);
978  bIsPileup=kTRUE;
979  }
980 
981  // check vertex consistency
982  const AliAODVertex* vtTrc = fAOD->GetPrimaryVertex();
983  const AliAODVertex* vtSPD = fAOD->GetPrimaryVertexSPD();
984 
985  if(vtTrc->GetNContributors() < 2 || vtSPD->GetNContributors()<1) {
986  fHistPileUpCount->Fill(5.5);
987  bIsPileup=kTRUE;
988  }
989 
990  double covTrc[6], covSPD[6];
991  vtTrc->GetCovarianceMatrix(covTrc);
992  vtSPD->GetCovarianceMatrix(covSPD);
993 
994  double dz = vtTrc->GetZ() - vtSPD->GetZ();
995 
996  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
997  double errTrc = TMath::Sqrt(covTrc[5]);
998  double nsigTot = dz/errTot;
999  double nsigTrc = dz/errTrc;
1000 
1001  if(TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) {
1002  fHistPileUpCount->Fill(6.5);
1003  bIsPileup=kTRUE;
1004  }
1005 
1006  Float_t multTPC = 0; // tpc mult estimate
1007  //Float_t RefMultRaw = 0; // tpc mult estimate
1008  //Float_t RefMultCorr = 0; // tpc mult estimate
1009  Float_t RefMultRawFB = 0;
1010  Float_t RefMultCorrFB= 0;
1011 
1012  Float_t multTPCAll = 0; // tpc mult estimate
1013  Float_t multGlobal = 0; // global multiplicity
1014 
1015  Float_t multEtaNeg, multEtaPos, multEtaFull;
1016 
1017  Int_t ChTrk;
1018  Double_t etaTrk, phiTrk; // Never define eta as float, always double.!!
1019 
1020  Double_t fMaxPtEP = 5.0;
1021  Double_t fMinPtEP = 0.2;
1022  Double_t fMaxEtaEP = 0.8;
1023  Double_t fMinEtaEP = -0.8;
1024 
1025  multEtaNeg=0;
1026  multEtaPos=0;
1027  multEtaFull=0;
1028 
1029 
1030 //const Int_t nGoodTracks = fVevent->GetNumberOfTracks();
1031  const Int_t nGoodTracks = fAOD->GetNumberOfTracks();
1032 
1033  for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { //-------------------------
1034 
1035  AliAODTrack* AODtrack =dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
1036  if(!AODtrack) continue;
1037 
1038  if(AODtrack->TestFilterBit(128)) multTPCAll++; //A. Dobrin, no track cuts
1039 
1040  ptTrk = AODtrack->Pt();
1041  etaTrk = AODtrack->Eta();
1042  phiTrk = AODtrack->Phi();
1043  ChTrk = AODtrack->Charge();
1044  dEdx = AODtrack->GetDetPid()->GetTPCsignal();
1045 
1046  //cuts for Outlier as in FlowEvent Task
1047  if(!(AODtrack->TestFilterBit(1))) continue;
1048  if((ptTrk < 0.2) || (TMath::Abs(etaTrk) > 0.8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.1)) continue;
1049 
1050  multTPC++;
1051  //RefMultRaw++;
1052  //RefMultCorr += ptWgtMC;
1053 
1054 
1055 
1056  //cuts for EP calculation:
1057  if((ptTrk <= fMaxPtEP) && (ptTrk >= fMinPtEP) && (etaTrk <= fMaxEtaEP) && (etaTrk >= fMinEtaEP) && (dEdx >= 10.0) && (AODtrack->DCA() <= fDCAxyMax) && (AODtrack->ZAtDCA() <= fDCAzMax) && ( TMath::Abs(ChTrk) > 0))
1058  {
1059 
1060  ptWgtMC = 1.0;
1061 
1062  if(fFB_Efficiency_Cent[cent10bin]){
1063  ptBinMC = fFB_Efficiency_Cent[cent10bin]->FindBin(ptTrk); //Charge independent MC correction atm.
1064  ptWgtMC = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBinMC);
1065  }
1066 
1067  RefMultRawFB++;
1068  RefMultCorrFB += ptWgtMC;
1069 
1070  if(ptTrk >= fMinPtCut && ptTrk <= fMaxPtCut) {
1071 
1072  //------> Get NUA weights for EP <----------
1073  if(ChTrk>0){
1074  if(fHCorrectNUApos[cForNUA]){
1075  iBinNUA = fHCorrectNUApos[cForNUA]->FindBin(pVtxZ,phiTrk,etaTrk);
1076  WgtNUA = fHCorrectNUApos[cForNUA]->GetBinContent(iBinNUA);
1077  }
1078  else{ WgtNUA = 1.0; }
1079  }
1080  else{
1081  if(fHCorrectNUAneg[cForNUA]){
1082  iBinNUA = fHCorrectNUAneg[cForNUA]->FindBin(pVtxZ,phiTrk,etaTrk);
1083  WgtNUA = fHCorrectNUAneg[cForNUA]->GetBinContent(iBinNUA);
1084  }
1085  else{ WgtNUA = 1.0; }
1086  }
1087 
1088  if(etaTrk < -0.05){
1089  sumTPCQn2x[0] += WgtNUA*TMath::Cos(gPsiN*phiTrk);
1090  sumTPCQn2y[0] += WgtNUA*TMath::Sin(gPsiN*phiTrk);
1091  sumTPCQn3x[0] += WgtNUA*TMath::Cos(3*phiTrk);
1092  sumTPCQn3y[0] += WgtNUA*TMath::Sin(3*phiTrk);
1093  sumTPCQn4x[0] += WgtNUA*TMath::Cos(4*phiTrk);
1094  sumTPCQn4y[0] += WgtNUA*TMath::Sin(4*phiTrk);
1095  multEtaNeg++;
1096  }
1097  else if(etaTrk > 0.05){
1098  sumTPCQn2x[1] += WgtNUA*TMath::Cos(gPsiN*phiTrk);
1099  sumTPCQn2y[1] += WgtNUA*TMath::Sin(gPsiN*phiTrk);
1100  sumTPCQn3x[1] += WgtNUA*TMath::Cos(3*phiTrk);
1101  sumTPCQn3y[1] += WgtNUA*TMath::Sin(3*phiTrk);
1102  sumTPCQn4x[1] += WgtNUA*TMath::Cos(4*phiTrk);
1103  sumTPCQn4y[1] += WgtNUA*TMath::Sin(4*phiTrk);
1104  multEtaPos++;
1105  }
1106  sumTPCQn2x[3] += WgtNUA*TMath::Cos(gPsiN*phiTrk);
1107  sumTPCQn2y[3] += WgtNUA*TMath::Sin(gPsiN*phiTrk);
1108  sumTPCQn3x[3] += WgtNUA*TMath::Cos(3*phiTrk);
1109  sumTPCQn3y[3] += WgtNUA*TMath::Sin(3*phiTrk);
1110  sumTPCQn4x[3] += WgtNUA*TMath::Cos(4*phiTrk);
1111  sumTPCQn4y[3] += WgtNUA*TMath::Sin(4*phiTrk);
1112  multEtaFull++;
1113  }
1114  }
1115 
1116 
1117 
1118  Double_t b[2] = {-99., -99.};
1119  Double_t bCov[3] = {-99., -99., -99.};
1120  AliAODTrack copy(*AODtrack);
1121  if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
1122  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1123  multGlobal++;
1124 
1125  }//--- track loop outlier/PileUp ----
1126 
1127  Int_t multEsd = ((AliAODHeader*)fAOD->GetHeader())->GetNumberOfESDTracks();
1128 
1129  Float_t multESDTPCDiff = (Float_t) multEsd - fPileUpSlopeParm*multTPCAll;
1130 
1131  if(multESDTPCDiff > fPileUpConstParm) {
1132  fHistPileUpCount->Fill(7.5);
1133  bIsPileup=kTRUE;
1134  }
1135  else if(bIsPileup==kFALSE) {
1136  if(!fMultSelection->GetThisEventIsNotPileup()){
1137  fHistMultSelPUCount->Fill(0.5);
1138  bIsPileup=kTRUE;
1139  }
1140  if(!fMultSelection->GetThisEventIsNotPileupMV()){
1141  fHistMultSelPUCount->Fill(1.5);
1142  bIsPileup=kTRUE;
1143  }
1144  if(!fMultSelection->GetThisEventIsNotPileupInMultBins()){
1145  fHistMultSelPUCount->Fill(2.5);
1146  bIsPileup=kTRUE;
1147  }
1148  if(!fMultSelection->GetThisEventHasNoInconsistentVertices()){
1149  fHistMultSelPUCount->Fill(2.5);
1150  bIsPileup=kTRUE;
1151  }
1152  if(!fMultSelection->GetThisEventPassesTrackletVsCluster()){
1153  fHistMultSelPUCount->Fill(2.5);
1154  bIsPileup=kTRUE;
1155  }
1156  if(!fMultSelection->GetThisEventIsNotIncompleteDAQ()){
1157  fHistMultSelPUCount->Fill(2.5);
1158  bIsPileup=kTRUE;
1159  }
1160  if(!fMultSelection->GetThisEventHasGoodVertex2016()){
1161  fHistMultSelPUCount->Fill(2.5);
1162  bIsPileup=kTRUE;
1163  }
1164  if(bIsPileup) fHistPileUpCount->Fill(9.5);
1165  }
1166  //-----------------------------------------------------------------
1167 
1168 
1169 
1170  fHistTPCVsESDTrkBefore->Fill(multTPCAll,multEsd); //A. Dobrin
1171 
1172  fHistTPCvsGlobalMultBefore->Fill(multGlobal,multTPC);
1173 
1174 
1175  Bool_t bIsOutLier=kFALSE;
1176 
1177  if(multTPC < (-20.0+1.15*multGlobal) || multTPC > (200.+1.40*multGlobal)) { bIsOutLier = kTRUE;}
1178 
1179  fHistEventCount->Fill(stepCount); //6
1180  stepCount++;
1181 
1182 
1183  fHistTPConlyVsCL1Before->Fill(centrCL1,multTPCAll);
1184  fHistTPConlyVsV0MBefore->Fill(centrV0M,multTPCAll);
1185  fHistGlobalVsV0MBefore->Fill(centrV0M, multGlobal);
1186 
1187 
1188 
1189  //if bSkipPileUpCut is kTRUE then don't apply PileUp removal.
1190  if(!bSkipPileUpCut && bIsOutLier) return; //outlier TPC vs Global
1191 
1192  fHistTPCvsGlobalMultAfter->Fill(multGlobal,multTPC);
1193 
1194  fHistEventCount->Fill(stepCount); //7
1195  stepCount++;
1196 
1197  //cout<<"After Outlier = "<<fSkipOutlierCut<<" multTPC = "<<multTPC<<" multGlobal = "<<multGlobal<<endl;
1198 
1199 
1200  if(!bSkipPileUpCut && bIsPileup) return; //PileUp A. Dobrin
1201 
1202  fHistTPCVsESDTrkAfter->Fill(multTPCAll,multEsd);
1203 
1204  fHistEventCount->Fill(stepCount); //8
1205  stepCount++;
1206 
1207  //cout<<"After PU cut = "<<fSkipOutlierCut<<" multTPC = "<<multTPC<<" multGlobal = "<<multGlobal<<endl;
1208 
1209 
1210 
1211 
1212 
1213  /*
1214  Int_t icentBin = centrality;
1215  icentBin++;
1216 
1217  Float_t TPCmultLowLimit = hCentvsTPCmultCuts->GetBinContent(icentBin,1);
1218  Float_t TPCmultHighLimit = hCentvsTPCmultCuts->GetBinContent(icentBin,1);
1219 
1220  TPCmultLowLimit -= 5.0 * hCentvsTPCmultCuts->GetBinContent(icentBin,2); //mean - 5sigma
1221  TPCmultHighLimit += 5.0 * hCentvsTPCmultCuts->GetBinContent(icentBin,2); //mean + 5sigma
1222  //std::cout<<" Cent = "<<centrality<<"\t icent = "<<icentBin<<" low = "<<TPCmultLowLimit<<"\t high = "<<TPCmultHighLimit<<std::endl;
1223 
1224  //centrality outlier
1225  if(!bSkipPileUpCut){ if(multTPC<TPCmultLowLimit || multTPC>TPCmultHighLimit) return; }
1226  */
1227 
1228 
1229 
1230 
1231 
1232  fHistEventCount->Fill(stepCount); //9
1233  stepCount++;
1234 
1235 
1236 
1237  fHistTPConlyVsCL1After->Fill(centrCL1,multTPCAll);
1238  fHistTPConlyVsV0MAfter->Fill(centrV0M,multTPCAll);
1239  fHistGlobalVsV0MAfter->Fill(centrV0M, multGlobal);
1240 
1241 
1242  // MC corrected Refmult:
1243  fHistRawVsCorrMultFB->Fill(RefMultRawFB,RefMultCorrFB); // FB set by AddTask..
1244 
1245 
1246  Float_t EvtCent = centrality;
1247 
1248 
1249  if(multEtaNeg<2 || multEtaPos<2) return; //Minimum 2 tracks in each eta
1250 
1251  fHistEventCount->Fill(stepCount); //10
1252  stepCount++;
1253 
1254  //--------------------------------------------------------
1255 
1256 
1257 
1258 
1259 
1260 
1261 
1262 
1263 
1264 
1265 
1266 
1267 
1268 
1269 
1270 
1271 
1272 
1273  //--------------- cent CL1 <= 90 cut ------------------
1274  Int_t icentV0Qn = centrCL1; // cent CL1 used for V0 calibration.
1275  icentV0Qn += 1;
1276 
1277  if(icentV0Qn>90) return;
1278 
1279  fHistEventCount->Fill(stepCount); //11
1280  stepCount++;
1281 
1282 
1283 
1284 
1285  //-------- V0M info ---------------
1286  const AliAODVZERO *fAODV0 = fAOD->GetVZEROData();
1287 
1288  //do v0m recentering
1289  Double_t QxanCor = 0, QyanCor = 0;
1290  Double_t QxcnCor = 0, QycnCor = 0;
1291 
1292  Double_t Qxan3 = 0., Qyan3 = 0.;
1293  Double_t Qxcn3 = 0., Qycn3 = 0.;
1294  Double_t Qxan2 = 0., Qyan2 = 0.;
1295  Double_t Qxcn2 = 0., Qycn2 = 0.;
1296 
1297  Double_t phiV0;
1298  Float_t fMultv0 = 0;
1299  Float_t sumMa = 0;
1300  Float_t sumMc = 0;
1301 
1302  for(int iV0 = 0; iV0 < 64; iV0++) { //0-31 is V0C, 32-63 VOA
1303 
1304  fMultv0 = fAODV0->GetMultiplicity(iV0);
1305 
1306  if(fHCorrectV0M){
1307  fMultv0 = fMultv0 * fHCorrectV0M->GetBinContent(iV0+1); // Gain Correction
1308  //cout<<"info: run = "<<runNumber<<" cent = "<<centrCL1<<"\t channel = "<<iV0<<" gain = "<<fHCorrectV0M->GetBinContent(iV0+1)<<endl;
1309  }
1310 
1311  fV0MultChVsRun->Fill(iV0+0.5,centrCL1,fMultv0);
1312  //fV0MultChVsRun->Fill(iV0+0.5,runindex,fMultv0);
1313 
1314  phiV0 = TMath::PiOver4()*(0.5 + iV0 % 8);
1315 
1316  if(iV0 < 32){
1317  Qxcn2 += TMath::Cos(2*phiV0) * fMultv0;
1318  Qycn2 += TMath::Sin(2*phiV0) * fMultv0;
1319  Qxcn3 += TMath::Cos(3*phiV0) * fMultv0;
1320  Qycn3 += TMath::Sin(3*phiV0) * fMultv0;
1321  sumMc += fMultv0;
1322  }
1323  else if(iV0 >= 32){
1324  Qxan2 += TMath::Cos(2*phiV0) * fMultv0;
1325  Qyan2 += TMath::Sin(2*phiV0) * fMultv0;
1326  Qxan3 += TMath::Cos(3*phiV0) * fMultv0;
1327  Qyan3 += TMath::Sin(3*phiV0) * fMultv0;
1328  sumMa += fMultv0;
1329  }
1330  }//----- channel loop ----------
1331 
1332 
1333  if(gPsiN==3){
1334  QxanCor = Qxan3/sumMa; //3rd order event plane
1335  QyanCor = Qyan3/sumMa;
1336  QxcnCor = Qxcn3/sumMc;
1337  QycnCor = Qycn3/sumMc;
1338 
1339  if(fHAvgerageQnV0C && fHAvgerageQnV0A && icentV0Qn < 91){
1340  QxanCor -= fHAvgerageQnV0A->GetBinContent(icentV0Qn,3); //x = Cos
1341  QxcnCor -= fHAvgerageQnV0C->GetBinContent(icentV0Qn,3); //x = Cos
1342  QyanCor -= fHAvgerageQnV0A->GetBinContent(icentV0Qn,4); //y = Sin
1343  QycnCor -= fHAvgerageQnV0C->GetBinContent(icentV0Qn,4); //y = Sin
1344  }
1345  //printf("\n .... I am using my own V0 gain correction for Psi3...\n");
1346  }
1347  else{
1348  QxanCor = Qxan2/sumMa; //2nd order Event plane
1349  QyanCor = Qyan2/sumMa;
1350  QxcnCor = Qxcn2/sumMc;
1351  QycnCor = Qycn2/sumMc;
1352 
1353  if(fHAvgerageQnV0C && fHAvgerageQnV0A && icentV0Qn < 91){
1354  QxanCor -= fHAvgerageQnV0A->GetBinContent(icentV0Qn,1); //x = Cos
1355  QxcnCor -= fHAvgerageQnV0C->GetBinContent(icentV0Qn,1); //x = Cos
1356  QyanCor -= fHAvgerageQnV0A->GetBinContent(icentV0Qn,2); //y = Sin
1357  QycnCor -= fHAvgerageQnV0C->GetBinContent(icentV0Qn,2); //y = Sin
1358  }
1359  //printf("\n .... I am using my own V0 gain correction for Psi2...\n");
1360  }
1361 
1362 
1363  //------------ For V0-Qn Recenter and Event plane ----------
1364  fV0CQ2xVsCentRun->Fill(centrCL1,Qxcn2/sumMc);
1365  fV0CQ2yVsCentRun->Fill(centrCL1,Qycn2/sumMc);
1366  fV0AQ2xVsCentRun->Fill(centrCL1,Qxan2/sumMa);
1367  fV0AQ2yVsCentRun->Fill(centrCL1,Qyan2/sumMa);
1368 
1369  fV0CQ3xVsCentRun->Fill(centrCL1,Qxcn3/sumMc);
1370  fV0CQ3yVsCentRun->Fill(centrCL1,Qycn3/sumMc);
1371  fV0AQ3xVsCentRun->Fill(centrCL1,Qxan3/sumMa);
1372  fV0AQ3yVsCentRun->Fill(centrCL1,Qyan3/sumMa);
1373 
1374  PsiNV0C = 1.0/gPsiN*( TMath::ATan2(QycnCor,QxcnCor) + TMath::Pi() );
1375  //if(PsiNV0C<0.) PsiNV0C += 2*TMath::Pi()/gPsiN;
1376  PsiNV0A = 1.0/gPsiN*( TMath::ATan2(QyanCor,QxanCor) + TMath::Pi() );
1377  //if(PsiNV0A<0.) PsiNV0A += 2*TMath::Pi()/gPsiN;
1378 
1379  fHV0CEventPlaneVsCent->Fill(EvtCent,PsiNV0C);
1380  fHV0AEventPlaneVsCent->Fill(EvtCent,PsiNV0A);
1381  //-------------------------------------------------------------
1382 
1383 
1384 
1385 
1386 
1387 
1388 
1389 
1390 
1391 
1392 
1393 
1394  //----- Psi_N from TPC sub & full event plane -------
1395  /*
1396  PsiNTPCA = (1.0/gPsiN)*( TMath::ATan2(sumTPCQn2y[0],sumTPCQn2x[0] )) ; // negetive eta
1397  if(PsiNTPCA<0.) PsiNTPCA += 2*TMath::Pi()/gPsiN;
1398  PsiNTPCC = (1.0/gPsiN)*( TMath::ATan2(sumTPCQn2y[1],sumTPCQn2x[1] )) ; // positive eta
1399  if(PsiNTPCC<0.) PsiNTPCC += 2*TMath::Pi()/gPsiN; */
1400 
1401  // Enable periodicity PsiN directly:
1402  PsiNTPCA = (1.0/gPsiN)*( TMath::ATan2(sumTPCQn2y[0],sumTPCQn2x[0]) + TMath::Pi() ) ; // negetive eta
1403  //if(PsiNTPCA<0.) PsiNTPCA += 2*TMath::Pi()/gPsiN;
1404  PsiNTPCC = (1.0/gPsiN)*( TMath::ATan2(sumTPCQn2y[1],sumTPCQn2x[1]) + TMath::Pi() ) ; // positive eta
1405  //if(PsiNTPCC<0.) PsiNTPCC += 2*TMath::Pi()/gPsiN;
1406 
1407 
1408  //fHTPCAEventPlaneVsCent->Fill(EvtCent,TMath::Cos(PsiNTPCA-PsiNTPCC));
1409  //fHTPCCEventPlaneVsCent->Fill(EvtCent,TMath::Cos(PsiNTPCC-PsiNV0C));
1410 
1411  //fHTPCAEventPlaneVsCent->Fill(EvtCent,PsiNTPCA);
1412  //fHTPCCEventPlaneVsCent->Fill(EvtCent,PsiNTPCC);
1413 
1414  fEventCount++;
1415 
1416 
1417 
1418  PsiNTPCF = (1.0/gPsiN)*( TMath::ATan2(sumTPCQn2y[3],sumTPCQn2x[3] )) ; // FUll TPC eta
1419  if(PsiNTPCF<0.) PsiNTPCF += 2*TMath::Pi()/gPsiN;
1420 
1421  //Psi3TPCA = (1.0/3)*( TMath::ATan2(sumTPCQn3y[0],sumTPCQn3x[0] )) ; // negetive eta ; Psi3 not needed now.
1422  //if(Psi3TPCA<0.) Psi3TPCA += 2*TMath::Pi()/3;
1423  //Psi3TPCC = (1.0/3)*( TMath::ATan2(sumTPCQn3y[1],sumTPCQn3x[1] )) ; // positive eta
1424  //if(Psi3TPCC<0.) Psi3TPCC += 2*TMath::Pi()/3;
1425  //----------------------------------------------------------------
1426 
1427 
1428 
1429 
1430 
1431 
1432 
1433 
1434 
1435 
1436 
1437  //---- Copies of TPC-Q vectors to remove track -by- track auto-correlation -----
1438 
1439  Double_t sumQxTPCneg;
1440  Double_t sumQyTPCneg;
1441  Double_t sumQxTPCpos;
1442  Double_t sumQyTPCpos;
1443  Double_t sumQxTPCneg2;
1444  Double_t sumQyTPCneg2;
1445  Double_t sumQxTPCpos2;
1446  Double_t sumQyTPCpos2;
1447 
1448 
1449 // Fill Centrality for run-by-run Event weight correction.
1450  fCentDistAfter->Fill(centrality);
1451 
1452 
1453 
1454 
1455 
1456 
1457  //--------- Track variable for PID/Charge studies ----------------
1458  Double_t PDGmassPion = 0.13957;
1459  Double_t PDGmassKaon = 0.49368;
1460  Double_t PDGmassProton = 0.93827;
1461 
1462  PDGmassProton *= PDGmassProton;
1463  PDGmassPion *= PDGmassPion;
1464  PDGmassKaon *= PDGmassKaon;
1465 
1466  Double_t mass, mom, beta, dEdx1;
1467  Double_t length, tofTime, dEdx2;
1468  Double_t dPhi1,dPhi2,dPt1,dPt2,dEta1,dEta2;
1469  Double_t ptw1, ptw2, w1NUA, w2NUA;
1470  Double_t nSigTOFpion, nSigTPCpion;
1471  Double_t nSigTOFkaon, nSigTPCkaon;
1472  Double_t nSigTOFproton,nSigTPCproton;
1473  Double_t nSigTOFpion2, nSigTPCpion2;
1474  Double_t nSigTOFkaon2, nSigTPCkaon2;
1475  Double_t nSigTOFproton2,nSigTPCproton2;
1476 
1477  Double_t c = TMath::C()*1.E-9;// velocity of light m/ns
1478  //Double_t dcaXY, dcaZ ;
1479  Double_t probMis, WgtEP;
1480  Int_t TOFmatch=0;
1481  Int_t ptBin,gCharge1,gCharge2;
1482 
1483 
1484 
1485  //----------- Set the desired Harmonic ------------
1486  Int_t n = gN;
1487  Int_t m = gM;
1488  Int_t p =n+m;
1489  //------------------------------------------------
1490 
1491 
1492  Int_t skipPairHBT = 0;
1493 
1494 
1495 
1496 
1497 
1498 
1499  Double_t ptwPion1,ptwKaon1,ptwProton1;
1500  Double_t ptwPion2,ptwKaon2,ptwProton2;
1501  Double_t wNUAPion1,wNUAKaon1,wNUAProton1;
1502  Double_t wNUAPion2,wNUAKaon2,wNUAProton2;
1503 
1504  Double_t WgtEPPion = 1.0;
1505  Double_t WgtEPKaon = 1.0;
1506  Double_t WgtEPProton = 1.0;
1507 
1508 
1509  Bool_t isPion1 = kFALSE;
1510  Bool_t isKaon1 = kFALSE;
1511  Bool_t isProton1 = kFALSE;
1512  Bool_t isPion2 = kFALSE;
1513  Bool_t isKaon2 = kFALSE;
1514  Bool_t isProton2 = kFALSE;
1515 
1516  Int_t multPOI1st = 0;
1517  Int_t multPOI2nd = 0;
1518 
1519 
1520  const Int_t maxTrack = 40000;
1521  Float_t nSigPionTPC[maxTrack] = {0.,};
1522  Float_t nSigKaonTPC[maxTrack] = {0.,};
1523  Float_t nSigProtonTPC[maxTrack] = {0.,};
1524  Float_t nSigPionTOF[maxTrack] = {0.,};
1525  Float_t nSigKaonTOF[maxTrack] = {0.,};
1526  Float_t nSigProtonTOF[maxTrack] = {0.,};
1527 
1528 //Dont break segment for higher tracks:
1529  if(ntracks>maxTrack) return;
1530 
1531  //vector<float> testVar;
1532  //testVar.reserve(20000);
1533  //testVar.reserve(ntracks+1);
1534 
1535 /*vector<float> nSigPionTPC;
1536  vector<float> nSigKaonTPC;
1537  vector<float> nSigProtonTPC;
1538  vector<float> nSigPionTOF;
1539  vector<float> nSigKaonTOF;
1540  vector<float> nSigProtonTOF;
1541 
1542  nSigPionTPC.reserve(ntracks+1);
1543  nSigKaonTPC.reserve(ntracks+1);
1544  nSigProtonTPC.reserve(ntracks+1);
1545  nSigPionTOF.reserve(ntracks+1);
1546  nSigKaonTOF.reserve(ntracks+1);
1547  nSigProtonTOF.reserve(ntracks+1);
1548  */
1549 
1550 
1551  //Calling fPIDResponse in nested loop is CPU expensive.
1552  //Store nSigma values in a array:
1553 
1554  for(Int_t itrack = 0; itrack < ntracks; itrack++) {
1555 
1556  AliAODTrack *trackForPID=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(itrack));
1557 
1558  // Array method:
1559  if(trackForPID){
1560  nSigPionTPC[itrack] = fPIDResponse->NumberOfSigmasTPC(trackForPID, AliPID::kPion);
1561  nSigKaonTPC[itrack] = fPIDResponse->NumberOfSigmasTPC(trackForPID, AliPID::kKaon);
1562  nSigProtonTPC[itrack] = fPIDResponse->NumberOfSigmasTPC(trackForPID, AliPID::kProton);
1563  nSigPionTOF[itrack] = fPIDResponse->NumberOfSigmasTOF(trackForPID, AliPID::kPion);
1564  nSigKaonTOF[itrack] = fPIDResponse->NumberOfSigmasTOF(trackForPID, AliPID::kKaon);
1565  nSigProtonTOF[itrack] = fPIDResponse->NumberOfSigmasTOF(trackForPID, AliPID::kProton);
1566  }
1567  else{
1568  nSigPionTPC[itrack] = -99;
1569  nSigKaonTPC[itrack] = -99;
1570  nSigProtonTPC[itrack] = -99;
1571  nSigPionTOF[itrack] = -99;
1572  nSigKaonTOF[itrack] = -99;
1573  nSigProtonTOF[itrack] = -99;
1574  }
1575 
1576  /*
1577  // Vector method:
1578  if(trackForPID){
1579  nSigPionTPC.push_back(fPIDResponse->NumberOfSigmasTPC(trackForPID, AliPID::kPion));
1580  nSigKaonTPC.push_back(fPIDResponse->NumberOfSigmasTPC(trackForPID, AliPID::kKaon));
1581  nSigProtonTPC.push_back(fPIDResponse->NumberOfSigmasTPC(trackForPID, AliPID::kProton));
1582 
1583  nSigPionTOF.push_back(fPIDResponse->NumberOfSigmasTOF(trackForPID, AliPID::kPion));
1584  nSigKaonTOF.push_back(fPIDResponse->NumberOfSigmasTOF(trackForPID, AliPID::kKaon));
1585  nSigProtonTOF.push_back(fPIDResponse->NumberOfSigmasTOF(trackForPID, AliPID::kProton));
1586  }
1587  else{
1588  nSigPionTPC.push_back(-100.);
1589  nSigKaonTPC.push_back(-100.);
1590  nSigProtonTPC.push_back(-100.);
1591 
1592  nSigPionTOF.push_back(-100.);
1593  nSigKaonTOF.push_back(-100.);
1594  nSigProtonTOF.push_back(-100.);
1595  }
1596  */
1597 
1598  }//1st track loop for PID storing.
1599 
1600 
1601 
1602 
1603 
1604 
1605 
1606 
1607 
1608  for(Int_t itrack = 0; itrack < ntracks; itrack++) {
1609 
1610  AliAODTrack *track=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(itrack));
1611  if(!track) continue;
1612 
1613  if(!track->TestFilterBit(fFilterBit)) continue;
1614 
1615 
1616 
1617  mom = track->P();
1618  dPt1 = track->Pt();
1619  dPhi1 = track->Phi();
1620  dEta1 = track->Eta();
1621  gCharge1 = track->Charge();
1622  dEdx1 = track->GetDetPid()->GetTPCsignal();
1623  //dcaXY = track->DCA();
1624  //dcaZ = track->ZAtDCA();
1625 
1626 
1627  //-------------- Check TOF status ------------------
1628  AliPIDResponse::EDetPidStatus status;
1629  status = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,track);
1630  TOFmatch = 0;
1631  if(status==AliPIDResponse::kDetPidOk){
1632  TOFmatch++;
1633  }
1634  probMis = fPIDResponse->GetTOFMismatchProbability(track);
1635  fHistTOFMatchCount->Fill(TOFmatch,probMis);
1636 
1637  mass = -9.9;
1638  beta = -0.9;
1639 
1640  if(TOFmatch>0 && probMis < 0.01) {
1641  //This conditions are called when detector status is checked above :
1642  //if((track->IsOn(AliAODTrack::kTOFin)) && (track->IsOn(AliAODTrack::kTIME)) && (track->IsOn(AliAODTrack::kTOFout))) {
1643  //if((track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid))) { //Naghmeh used it
1644  tofTime = track->GetTOFsignal(); // in pico seconds
1645  length = track->GetIntegratedLength();
1646  tofTime = tofTime*1.E-3; // ns
1647  if (tofTime <= 0) tofTime = 9999;
1648  length = length*0.01; // in meters
1649  tofTime = tofTime*c;
1650  beta = length/tofTime;
1651  mass = mom*mom*(1./(beta*beta) - 1);
1652  }//------------ TOF signal -------------------------
1653 
1654 
1655 
1656 
1657  //QA histograms before applying track cuts:
1658  fHistEtaPtBefore->Fill(dEta1,dPt1);
1659  fHistTPCdEdxvsPBefore->Fill(mom*gCharge1,dEdx1);
1660  fHistTOFBetavsPBefore->Fill(mom*gCharge1,beta);
1661  fHistTOFMassvsPtBefore->Fill(dPt1*gCharge1,mass);
1662 
1663  //-------- Apply Default track cuts for analysis: ---------
1664  if((dPt1 > fMaxPtCut) || (dPt1 < fMinPtCut) || (dEta1 > fMaxEtaCut) || (dEta1 < fMinEtaCut) || (dEdx1 < fdEdxMin) || (track->GetTPCNcls() < fTPCclustMin) || (track->Chi2perNDF() < fTrkChi2Min) || (track->DCA() > fDCAxyMax) || (track->ZAtDCA() > fDCAzMax) || !(TMath::Abs(gCharge1)))
1665  continue;
1666 
1667 
1668 
1669  /*//tested above in one if(...)
1670  if(track->GetDetPid()->GetTPCsignal() < 10.0) continue;
1671  if(track->GetTPCNcls() < 70) continue;
1672  if(track->Chi2perNDF() < 0.2) continue;
1673  if(dPt1 < fMinPtCut || dPt1 > fMaxPtCut) continue;
1674  if(eta< fMinEtaCut || eta > fMaxEtaCut) continue;
1675  if(!TMath::Abs(gCharge1)) continue;
1676  //if(!(track->TestFilterBit(fFilterBit))) continue; */
1677  //-----------------------------------------------------
1678 
1679 
1680  //============= charged hadron analysis: ============
1681  sumQxTPCneg = sumTPCQn2x[0]; // first copy to remove 1st track fron Q vector (AutoCorrelation)
1682  sumQyTPCneg = sumTPCQn2y[0];
1683  sumQxTPCpos = sumTPCQn2x[1];
1684  sumQyTPCpos = sumTPCQn2y[1];
1685 
1686  //--------------------- PID signals 1st track-------------------------
1687  /*nSigTOFpion = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
1688  nSigTOFkaon = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
1689  nSigTOFproton = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
1690  nSigTPCpion = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
1691  nSigTPCkaon = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
1692  nSigTPCproton = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);*/
1693 
1694  //Vector/Array both works same way:
1695  nSigTPCpion = nSigPionTPC[itrack];
1696  nSigTPCkaon = nSigKaonTPC[itrack];
1697  nSigTPCproton = nSigProtonTPC[itrack];
1698  nSigTOFpion = nSigPionTOF[itrack];
1699  nSigTOFkaon = nSigKaonTOF[itrack];
1700  nSigTOFproton = nSigProtonTOF[itrack];
1701 
1702  multPOI1st++;
1703 
1704  //if(itrack%100==0)
1705  //cout<<"Trk "<<itrack<<" pt1 = "<<dPt1<<"\tnSigPion = "<<nSigTPCpion<<"\tnSigKaon = "<<nSigTPCkaon <<"\tnSigprot = "<<nSigTPCproton<<endl;
1706 
1707  isPion1 = kFALSE;
1708  isKaon1 = kFALSE;
1709  isProton1 = kFALSE;
1710 
1711  //------> Pion
1712  if(dPt1<0.6 && TMath::Abs(nSigTPCpion)<=2.5){
1713  isPion1 = kTRUE;
1714  }
1715  else if(dPt1>=0.6 && dPt1<=2.0 && TMath::Abs(nSigTPCpion)<=2.5 && TMath::Abs(nSigTOFpion)<=2.0 ){
1716  isPion1 = kTRUE;
1717  }
1718  //------> Kaon
1719  if(dPt1<0.6 && TMath::Abs(nSigTPCkaon)<=2.5){
1720  isKaon1 = kTRUE;
1721  }
1722  else if(dPt1>=0.6 && dPt1<=2.0 && TMath::Abs(nSigTPCkaon)<=2.5 && TMath::Abs(nSigTOFkaon)<=2.0){
1723  isKaon1 = kTRUE;
1724  }
1725  //------> Proton
1726  if(dPt1<0.8 && TMath::Abs(nSigTPCproton)<=2.5){
1727  isProton1 = kTRUE;
1728  }
1729  else if(dPt1>=0.8 && dPt1<=3.5 && TMath::Abs(nSigTPCproton)<=2.5 && TMath::Abs(nSigTOFproton)<=2.5){
1730  isProton1 = kTRUE;
1731  }
1732 
1733  //-----------------------------------------------------------------
1734 
1735 
1736 
1737 
1738 
1739  //================= MC wgt and NUA wgt for PID =================
1740  ptwPion1 = 1.0;
1741  ptwKaon1 = 1.0;
1742  ptwProton1 = 1.0;
1743  wNUAPion1 = 1.0;
1744  wNUAKaon1 = 1.0;
1745  wNUAProton1 = 1.0;
1746 
1747  //------ get MC weight and NUA for Pion track1 --------------
1748  if(isPion1){
1749  if(fFB_Efficiency_Cent[cent10bin]){ // <-------------------- !!!! WARNING: use Pion Efficiency file when available.
1750  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt1);
1751  ptwPion1 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
1752  }
1753  else{ ptwPion1 = 1.0; }
1754 
1755  if(gCharge1>0){
1756  if(fHCorrectNUAposPion[cForNUA]){
1757  iBinNUA = fHCorrectNUAposPion[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1758  wNUAPion1 = fHCorrectNUAposPion[cForNUA]->GetBinContent(iBinNUA);
1759  }
1760  else{ wNUAPion1 = 1.0; }
1761  }
1762  else{
1763  if(fHCorrectNUAnegPion[cForNUA]){
1764  iBinNUA = fHCorrectNUAnegPion[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1765  wNUAPion1 = fHCorrectNUAnegPion[cForNUA]->GetBinContent(iBinNUA);
1766  }
1767  else{ wNUAPion1 = 1.0; }
1768  }
1769  }
1770 
1771 
1772  //------ get MC weight and NUA for Kaon track1 --------------
1773  if(isKaon1){
1774  if(fFB_Efficiency_Cent[cent10bin]){ // <-------------------- !!!! WARNING: use Kaon Efficiency file when available.
1775  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt1);
1776  ptwKaon1 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
1777  }
1778  else{ ptwKaon1 = 1.0; }
1779 
1780  if(gCharge1>0){
1781  if(fHCorrectNUAposKaon[cForNUA]){
1782  iBinNUA = fHCorrectNUAposKaon[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1783  wNUAKaon1 = fHCorrectNUAposKaon[cForNUA]->GetBinContent(iBinNUA);
1784  }
1785  else{ wNUAKaon1 = 1.0; }
1786  }
1787  else{
1788  if(fHCorrectNUAnegKaon[cForNUA]){
1789  iBinNUA = fHCorrectNUAnegKaon[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1790  wNUAKaon1 = fHCorrectNUAnegKaon[cForNUA]->GetBinContent(iBinNUA);
1791  }
1792  else{ wNUAKaon1 = 1.0; }
1793  }
1794  }
1795  //------ get MC weight and NUA for Proton track1 --------------
1796  if(isProton1){
1797  if(fFB_Efficiency_Cent[cent10bin]){ // <-------------------- !!!! WARNING: use Proton Efficiency file when available.
1798  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt1);
1799  ptwProton1 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
1800  }
1801  else{ ptwProton1 = 1.0; }
1802 
1803  if(gCharge1>0){
1804  if(fHCorrectNUAposProton[cForNUA]){
1805  iBinNUA = fHCorrectNUAposProton[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1806  wNUAProton1 = fHCorrectNUAposProton[cForNUA]->GetBinContent(iBinNUA);
1807  }
1808  else{ wNUAProton1 = 1.0; }
1809  }
1810  else{
1811  if(fHCorrectNUAnegProton[cForNUA]){
1812  iBinNUA = fHCorrectNUAnegProton[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1813  wNUAProton1 = fHCorrectNUAnegProton[cForNUA]->GetBinContent(iBinNUA);
1814  }
1815  else{ wNUAProton1 = 1.0; }
1816  }
1817  }
1818  //=========================== X ===============================
1819 
1820 
1821 
1822  //------ get MC weight and NUA for Charged track 1--------------
1823  ptw1 = 1.0;
1824 
1825  if(fFB_Efficiency_Cent[cent10bin]){
1826  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt1);
1827  ptw1 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
1828  }
1829  else{ ptw1 = 1.0; }
1830 
1831 
1832  if(gCharge1>0){
1833  if(fHCorrectNUApos[cForNUA]){
1834  iBinNUA = fHCorrectNUApos[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1835  w1NUA = fHCorrectNUApos[cForNUA]->GetBinContent(iBinNUA);
1836  }
1837  else{ w1NUA = 1.0; }
1838  }
1839  else{
1840  if(fHCorrectNUAneg[cForNUA]){
1841  iBinNUA = fHCorrectNUAneg[cForNUA]->FindBin(pVtxZ,dPhi1,dEta1);
1842  w1NUA = fHCorrectNUAneg[cForNUA]->GetBinContent(iBinNUA);
1843  }
1844  else{ w1NUA = 1.0; }
1845  }
1846 
1847 
1848 
1849  //-------- Remove track 1 from EP calculation ----------
1850  if(dEta1 < -0.05){
1851  sumQxTPCneg -= w1NUA*TMath::Cos(gPsiN*dPhi1);
1852  sumQyTPCneg -= w1NUA*TMath::Sin(gPsiN*dPhi1); // [0] = eta <-0.05
1853  }
1854  else if(dEta1 > 0.05){
1855  sumQxTPCpos -= w1NUA*TMath::Cos(gPsiN*dPhi1);
1856  sumQyTPCpos -= w1NUA*TMath::Sin(gPsiN*dPhi1); // [1] = eta > 0.05
1857  }
1858  //-----------------------------------------------------
1859 
1860 
1861 
1862 
1863 
1864 
1865 
1866  multPOI2nd = 0;
1867 
1868  //---2nd track loop (nested)---
1869  for(Int_t jtrack = 0; jtrack < ntracks; jtrack++) {
1870 
1871  if(jtrack==itrack) continue;
1872 
1873  AliAODTrack *track2=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(jtrack));
1874  if(!track2) continue;
1875 
1876  if(!(track2->TestFilterBit(fFilterBit))) continue;
1877 
1878  dPt2 = track2->Pt();
1879  dPhi2 = track2->Phi();
1880  dEta2 = track2->Eta();
1881  gCharge2= track2->Charge();
1882  dEdx2 = track2->GetDetPid()->GetTPCsignal();
1883 
1884  if((dPt2 > fMaxPtCut) || (dPt2 < fMinPtCut) || (dEta2 > fMaxEtaCut) || (dEta2 < fMinEtaCut) || (dEdx2 < fdEdxMin) || (track2->GetTPCNcls() < fTPCclustMin) || (track2->Chi2perNDF() < fTrkChi2Min) || (track2->DCA() > fDCAxyMax) || (track2->ZAtDCA() > fDCAzMax) || !(TMath::Abs(gCharge2)))
1885  continue;
1886 
1887 
1888 
1889 
1890  /*
1891  if(track2->GetDetPid()->GetTPCsignal() < 10.0) continue;
1892  if(track2->GetTPCNcls() < 70) continue;
1893  if(track2->Chi2perNDF() < 0.2) continue;
1894  if(dPt2 < fMinPtCut || dPt2 > fMaxPtCut) continue;
1895  if(dEta2< fMinEtaCut || dEta2 > fMaxEtaCut) continue;
1896  if(!TMath::Abs(gCharge2)) continue;
1897  //if(!(track2->TestFilterBit(fFilterBit))) continue; */
1898  //-----------------------------------------------------
1899 
1900  //calling fPIDResponse is too costly for CPU
1901  //--------------------- PID signals 2nd track-------------------------
1902  /*nSigTOFpion2 = fPIDResponse->NumberOfSigmasTOF(track2, AliPID::kPion);
1903  nSigTOFkaon2 = fPIDResponse->NumberOfSigmasTOF(track2, AliPID::kKaon);
1904  nSigTOFproton2 = fPIDResponse->NumberOfSigmasTOF(track2, AliPID::kProton);
1905  nSigTPCpion2 = fPIDResponse->NumberOfSigmasTPC(track2, AliPID::kPion);
1906  nSigTPCkaon2 = fPIDResponse->NumberOfSigmasTPC(track2, AliPID::kKaon);
1907  nSigTPCproton2 = fPIDResponse->NumberOfSigmasTPC(track2, AliPID::kProton); */
1908 
1909 
1910 
1911  //Vector/Array both works same way:
1912  nSigTPCpion2 = nSigPionTPC[jtrack];
1913  nSigTPCkaon2 = nSigKaonTPC[jtrack];
1914  nSigTPCproton2 = nSigProtonTPC[jtrack];
1915  nSigTOFpion2 = nSigPionTOF[jtrack];
1916  nSigTOFkaon2 = nSigKaonTOF[jtrack];
1917  nSigTOFproton2 = nSigProtonTOF[jtrack];
1918 
1919 
1920  multPOI2nd++;
1921 
1922 
1923  isPion2 = kFALSE;
1924  isKaon2 = kFALSE;
1925  isProton2 = kFALSE;
1926 
1927  //------> Pion 2
1928  if(dPt2<0.6 && TMath::Abs(nSigTPCpion2)<=2.5){
1929  isPion2 = kTRUE;
1930  }
1931  else if(dPt2>=0.6 && dPt2<=2.0 && TMath::Abs(nSigTPCpion2)<=2.5 && TMath::Abs(nSigTOFpion2)<=2.0 ){
1932  isPion2 = kTRUE;
1933  }
1934  //------> Kaon 2
1935  if(dPt2<0.6 && TMath::Abs(nSigTPCkaon2)<=2.5){
1936  isKaon2 = kTRUE;
1937  }
1938  else if(dPt2>=0.6 && dPt2<=2.0 && TMath::Abs(nSigTPCkaon2)<=2.5 && TMath::Abs(nSigTOFkaon2)<=2.0){
1939  isKaon2 = kTRUE;
1940  }
1941  //------> Proton 2
1942  if(dPt2<0.8 && TMath::Abs(nSigTPCproton2)<=2.5){
1943  isProton2 = kTRUE;
1944  }
1945  else if(dPt2>=0.8 && dPt2<=3.5 && TMath::Abs(nSigTPCproton2)<=2.5 && TMath::Abs(nSigTOFproton2)<=2.5){
1946  isProton2 = kTRUE;
1947  }
1948  //----------------------------------------------------------------
1949 
1950 
1951 
1952 
1953 
1954  //================= MC wgt and NUA wgt for PID =================
1955  WgtEPPion = 1.0;
1956  WgtEPKaon = 1.0;
1957  WgtEPProton = 1.0;
1958 
1959  ptwPion2 = 1.0;
1960  ptwKaon2 = 1.0;
1961  ptwProton2 = 1.0;
1962  wNUAPion2 = 1.0;
1963  wNUAKaon2 = 1.0;
1964  wNUAProton2 = 1.0;
1965 
1966  //------ get MC weight and NUA for Pion track2 --------------
1967  if(isPion2){
1968  if(fFB_Efficiency_Cent[cent10bin]){ // <-------------------- !!!! WARNING: use Pion Efficiency file when available.
1969  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt2);
1970  ptwPion2 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
1971  }
1972  else{ ptwPion2 = 1.0; }
1973 
1974  if(gCharge2>0){
1975  if(fHCorrectNUAposPion[cForNUA]){
1976  iBinNUA = fHCorrectNUAposPion[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
1977  wNUAPion2 = fHCorrectNUAposPion[cForNUA]->GetBinContent(iBinNUA);
1978  }
1979  else{ wNUAPion2 = 1.0; }
1980  }
1981  else{
1982  if(fHCorrectNUAnegPion[cForNUA]){
1983  iBinNUA = fHCorrectNUAnegPion[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
1984  wNUAPion2 = fHCorrectNUAnegPion[cForNUA]->GetBinContent(iBinNUA);
1985  }
1986  else{ wNUAPion2 = 1.0; }
1987  }
1988 
1989  WgtEPPion = ptwPion1*ptwPion2*wNUAPion1*wNUAPion2;
1990  }
1991 
1992  //------ get MC weight and NUA for Kaon track2 --------------
1993  if(isKaon2){
1994  if(fFB_Efficiency_Cent[cent10bin]){ // <-------------------- !!!! WARNING: use Kaon Efficiency file when available.
1995  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt2);
1996  ptwKaon2 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
1997  }
1998  else{ ptwKaon2 = 1.0; }
1999 
2000  if(gCharge2>0){
2001  if(fHCorrectNUAposKaon[cForNUA]){
2002  iBinNUA = fHCorrectNUAposKaon[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
2003  wNUAKaon2 = fHCorrectNUAposKaon[cForNUA]->GetBinContent(iBinNUA);
2004  }
2005  else{ wNUAKaon2 = 1.0; }
2006  }
2007  else{
2008  if(fHCorrectNUAnegKaon[cForNUA]){
2009  iBinNUA = fHCorrectNUAnegKaon[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
2010  wNUAKaon2 = fHCorrectNUAnegKaon[cForNUA]->GetBinContent(iBinNUA);
2011  }
2012  else{ wNUAKaon2 = 1.0; }
2013  }
2014 
2015  WgtEPKaon = ptwKaon1*ptwKaon2*wNUAKaon1*wNUAKaon2;
2016  }
2017  //------ get MC weight and NUA for Proton track2 --------------
2018  if(isProton2){
2019  if(fFB_Efficiency_Cent[cent10bin]){ // <-------------------- !!!! WARNING: use Proton Efficiency file when available.
2020  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt2);
2021  ptwProton2 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
2022  }
2023  else{ ptwProton2 = 1.0; }
2024 
2025  if(gCharge2>0){
2026  if(fHCorrectNUAposProton[cForNUA]){
2027  iBinNUA = fHCorrectNUAposProton[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
2028  wNUAProton2 = fHCorrectNUAposProton[cForNUA]->GetBinContent(iBinNUA);
2029  }
2030  else{ wNUAProton2 = 1.0; }
2031  }
2032  else{
2033  if(fHCorrectNUAnegProton[cForNUA]){
2034  iBinNUA = fHCorrectNUAnegProton[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
2035  wNUAProton2 = fHCorrectNUAnegProton[cForNUA]->GetBinContent(iBinNUA);
2036  }
2037  else{ wNUAProton2 = 1.0; }
2038  }
2039 
2040  WgtEPProton = ptwProton1*ptwProton2*wNUAProton1*wNUAProton2;
2041  }
2042  //========================== X ================================
2043 
2044 
2045 
2046 
2047 
2048  //------ get MC weight and NUA (Charged) for track 2--------------
2049  WgtEP = 1.0;
2050  ptw2 = 1.0;
2051 
2052  if(fFB_Efficiency_Cent[cent10bin]){
2053  ptBin = fFB_Efficiency_Cent[cent10bin]->FindBin(dPt2);
2054  ptw2 = 1.0/fFB_Efficiency_Cent[cent10bin]->GetBinContent(ptBin);
2055  }
2056  else{ ptw2 = 1.0; }
2057 
2058  if(gCharge2>0){
2059  if(fHCorrectNUApos[cForNUA]){
2060  iBinNUA = fHCorrectNUApos[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
2061  w2NUA = fHCorrectNUApos[cForNUA]->GetBinContent(iBinNUA);
2062  }
2063  else{ w2NUA = 1.0; }
2064  }
2065  else{
2066  if(fHCorrectNUAneg[cForNUA]){
2067  iBinNUA = fHCorrectNUAneg[cForNUA]->FindBin(pVtxZ,dPhi2,dEta2);
2068  w2NUA = fHCorrectNUAneg[cForNUA]->GetBinContent(iBinNUA);
2069  }
2070  else{ w2NUA = 1.0; }
2071  }
2072 
2073  //---------- Remove track2 from EP calculation ---------
2074  sumQxTPCneg2 = sumQxTPCneg; //second copy to remove 2nd track from EP
2075  sumQyTPCneg2 = sumQyTPCneg;
2076  sumQxTPCpos2 = sumQxTPCpos;
2077  sumQyTPCpos2 = sumQyTPCpos;
2078  //------------------------------------------------------
2079 
2080 
2081  if(dEta2 < -0.05){
2082  sumQxTPCneg2 -= w2NUA*TMath::Cos(gPsiN*dPhi2);
2083  sumQyTPCneg2 -= w2NUA*TMath::Sin(gPsiN*dPhi2); // [0] = eta <-0.05
2084  }
2085  else if(dEta2 > 0.05){
2086  sumQxTPCpos2 -= w2NUA*TMath::Cos(gPsiN*dPhi2);
2087  sumQyTPCpos2 -= w2NUA*TMath::Sin(gPsiN*dPhi2); // [1] = eta > 0.05
2088  }
2089 
2090  /*
2091  // track by track EP:
2092  PsiNTPCA = (1.0/gPsiN)*( TMath::ATan2(sumQyTPCneg2,sumQxTPCneg2) );
2093  if(PsiNTPCA<0.) PsiNTPCA += 2*TMath::Pi()/gPsiN;
2094  PsiNTPCC = (1.0/gPsiN)*( TMath::ATan2(sumQyTPCpos2,sumQxTPCpos2) );
2095  if(PsiNTPCC<0.) PsiNTPCC += 2*TMath::Pi()/gPsiN; */
2096 
2097  PsiNTPCA = (1.0/gPsiN)*( TMath::ATan2(sumQyTPCneg2,sumQxTPCneg2) + TMath::Pi() );
2098  //if(PsiNTPCA<0.) PsiNTPCA += 2*TMath::Pi()/gPsiN;
2099  PsiNTPCC = (1.0/gPsiN)*( TMath::ATan2(sumQyTPCpos2,sumQxTPCpos2) + TMath::Pi() );
2100  //if(PsiNTPCC<0.) PsiNTPCC += 2*TMath::Pi()/gPsiN;
2101 
2102 
2103  fHTPCAEventPlaneVsCent->Fill(EvtCent,PsiNTPCA);
2104  fHTPCCEventPlaneVsCent->Fill(EvtCent,PsiNTPCC);
2105 
2106 
2107  // combined weight for EP:
2108  WgtEP = ptw1*ptw2*w1NUA*w2NUA;
2109  //-----------------------------------------------------------
2110 
2111 
2112  //cout<<", passes 2 ";
2113 
2114  if(gCharge1!=gCharge2) {
2115  fHist_Corr3p_EP_Norm_PN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEP);
2116  fHist_Corr3p_EP_Norm_PN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEP);
2117  fHist_Corr3p_EP_Norm_PN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEP);
2118  fHist_Corr3p_EP_Norm_PN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEP);
2119 
2120  //-------------> PID CME ---------------
2121  //Pion:
2122  if(isPion1 && isPion2){
2123  fHist_Corr3p_Pion_EP_Norm_PN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPPion);
2124  fHist_Corr3p_Pion_EP_Norm_PN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPPion);
2125  fHist_Corr3p_Pion_EP_Norm_PN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPPion);
2126  fHist_Corr3p_Pion_EP_Norm_PN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPPion);
2127  }
2128  //Kaon:
2129  if(isKaon1 && isKaon2){
2130  fHist_Corr3p_Kaon_EP_Norm_PN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPKaon);
2131  fHist_Corr3p_Kaon_EP_Norm_PN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPKaon);
2132  fHist_Corr3p_Kaon_EP_Norm_PN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPKaon);
2133  fHist_Corr3p_Kaon_EP_Norm_PN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPKaon);
2134  }
2135  //Proton:
2136  if(isProton1 && isProton2){
2137  //cout<<"#pair pt1 = "<<dPt1<<"\tpt2 = "<<dPt2<<"\tnSigp1 = "<<nSigTPCproton<<"\tnSigp2 = "<<nSigTPCproton2<<endl;
2138  fHist_Corr3p_Proton_EP_Norm_PN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPProton);
2139  fHist_Corr3p_Proton_EP_Norm_PN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPProton);
2140  fHist_Corr3p_Proton_EP_Norm_PN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPProton);
2141  fHist_Corr3p_Proton_EP_Norm_PN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPProton);
2142  }
2143  //------------------------------------
2144  }
2145  else if(gCharge1>0 && gCharge2>0 && skipPairHBT==0) {
2146  fHist_Corr3p_EP_Norm_PP[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEP);
2147  fHist_Corr3p_EP_Norm_PP[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEP);
2148  fHist_Corr3p_EP_Norm_PP[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEP);
2149  fHist_Corr3p_EP_Norm_PP[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEP);
2150 
2151  //-------------> PID CME ---------------
2152  //Pion:
2153  if(isPion1 && isPion2){
2154  fHist_Corr3p_Pion_EP_Norm_PP[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPPion);
2155  fHist_Corr3p_Pion_EP_Norm_PP[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPPion);
2156  fHist_Corr3p_Pion_EP_Norm_PP[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPPion);
2157  fHist_Corr3p_Pion_EP_Norm_PP[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPPion);
2158  }
2159  //Kaon:
2160  if(isKaon1 && isKaon2){
2161  fHist_Corr3p_Kaon_EP_Norm_PP[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPKaon);
2162  fHist_Corr3p_Kaon_EP_Norm_PP[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPKaon);
2163  fHist_Corr3p_Kaon_EP_Norm_PP[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPKaon);
2164  fHist_Corr3p_Kaon_EP_Norm_PP[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPKaon);
2165  }
2166  //Proton:
2167  if(isProton1 && isProton2){
2168  fHist_Corr3p_Proton_EP_Norm_PP[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPProton);
2169  fHist_Corr3p_Proton_EP_Norm_PP[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPProton);
2170  fHist_Corr3p_Proton_EP_Norm_PP[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPProton);
2171  fHist_Corr3p_Proton_EP_Norm_PP[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPProton);
2172  }
2173  //------------------------------------
2174  }
2175 
2176  else if(gCharge1<0 && gCharge2<0 && skipPairHBT==0){
2177  fHist_Corr3p_EP_Norm_NN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEP);
2178  fHist_Corr3p_EP_Norm_NN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEP);
2179  fHist_Corr3p_EP_Norm_NN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEP);
2180  fHist_Corr3p_EP_Norm_NN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEP);
2181 
2182  //-------------> PID CME ---------------
2183  //Pion:
2184  if(isPion1 && isPion2){
2185  fHist_Corr3p_Pion_EP_Norm_NN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPPion);
2186  fHist_Corr3p_Pion_EP_Norm_NN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPPion);
2187  fHist_Corr3p_Pion_EP_Norm_NN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPPion);
2188  fHist_Corr3p_Pion_EP_Norm_NN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPPion);
2189  }
2190  //Kaon:
2191  if(isKaon1 && isKaon2){
2192  fHist_Corr3p_Kaon_EP_Norm_NN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPKaon);
2193  fHist_Corr3p_Kaon_EP_Norm_NN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPKaon);
2194  fHist_Corr3p_Kaon_EP_Norm_NN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPKaon);
2195  fHist_Corr3p_Kaon_EP_Norm_NN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPKaon);
2196  }
2197  //Proton:
2198  if(isProton1 && isProton2){
2199  fHist_Corr3p_Proton_EP_Norm_NN[QAindex][0]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0A), WgtEPProton);
2200  fHist_Corr3p_Proton_EP_Norm_NN[QAindex][1]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNV0C), WgtEPProton);
2201  fHist_Corr3p_Proton_EP_Norm_NN[QAindex][2]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCA),WgtEPProton);
2202  fHist_Corr3p_Proton_EP_Norm_NN[QAindex][3]->Fill(EvtCent, TMath::Cos(n*dPhi1 + m*dPhi2 - p*PsiNTPCC),WgtEPProton);
2203  }
2204  //------------------------------------
2205  }
2206 
2207  //cout<<", passes 3 WgtEP = "<<WgtEP<<endl;
2208  }//-------- nested track loop ends ------------------
2209 
2210 
2211  //cout<<" It doesn't break here 6 \n";
2212 
2213 
2214 
2215  //============ PID business starts here =============
2216 
2217  fHistEtaPtAfter->Fill(dEta1,dPt1);
2218 
2219  fHistTPCdEdxvsPAfter->Fill(mom*gCharge1,dEdx1);
2220 
2221  if(TOFmatch>0 && probMis < 0.01){
2222  fHistTOFBetavsPAfter->Fill(mom*gCharge1,beta);
2223  }
2224 
2225  //nSigTOFpion=fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
2226  //nSigTOFkaon=fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
2227  //nSigTOFproton=fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
2228 
2229  fHistTOFMatchCount->Fill(TOFmatch+2,nSigTOFpion);
2230  fHistTOFMatchCount->Fill(TOFmatch+4,nSigTOFkaon);
2231  fHistTOFMatchCount->Fill(TOFmatch+6,nSigTOFproton);
2232 
2233  if(!TOFmatch || probMis > 0.01){ // I dont want mismatched track in my signal distribution
2234  nSigTOFpion = -99;
2235  nSigTOFkaon = -99;
2236  nSigTOFproton = -99;
2237  }
2238 
2239  //nSigTPCpion = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
2240  //nSigTPCkaon = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
2241  //nSigTPCproton = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2242 
2243  //0=pi, 1=K, 2=Proton
2244  fHistTPCTOFnSigmavsPtAfter[0]->Fill(dPt1,nSigTPCpion,nSigTOFpion);
2245  fHistTPCTOFnSigmavsPtAfter[1]->Fill(dPt1,nSigTPCkaon,nSigTOFkaon);
2246  fHistTPCTOFnSigmavsPtAfter[2]->Fill(dPt1,nSigTPCproton,nSigTOFproton);
2247 
2248 
2249 
2250  //cout<<" It doesn't break here 7 \n";
2251 
2252 
2253 
2254 
2255  if(TMath::Abs(nSigTPCpion)<=fNSigmaCut){
2256  fHistTOFnSigmavsPtAfter[0]->Fill(dPt1*gCharge1,nSigTOFpion);
2257  fHistPtwithTPCNsigma[0]->Fill(dPt1*gCharge1);
2258  fHistTPCdEdxvsPtPIDAfter[0]->Fill(dPt1,dEdx1);
2259  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
2260  fHistPtwithTOFSignal[0]->Fill(dPt1*gCharge1);
2261  }
2262  }
2263  if(TMath::Abs(nSigTPCkaon)<=fNSigmaCut){
2264  fHistPtwithTPCNsigma[1]->Fill(dPt1*gCharge1);
2265  fHistTOFnSigmavsPtAfter[1]->Fill(dPt1*gCharge1,nSigTOFkaon);
2266  fHistTPCdEdxvsPtPIDAfter[1]->Fill(dPt1,dEdx1);
2267  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
2268  fHistPtwithTOFSignal[1]->Fill(dPt1*gCharge1);
2269  }
2270  }
2271  if(TMath::Abs(nSigTPCproton)<=fNSigmaCut){
2272  fHistPtwithTPCNsigma[2]->Fill(dPt1*gCharge1);
2273  fHistTOFnSigmavsPtAfter[2]->Fill(dPt1*gCharge1,nSigTOFproton);
2274  fHistTPCdEdxvsPtPIDAfter[2]->Fill(dPt1,dEdx1);
2275  if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
2276  fHistPtwithTOFSignal[2]->Fill(dPt1*gCharge1);
2277  }
2278  }
2279 
2280 
2281 
2282 
2283  //========> nSigmaTOF distribution for circular cut <===============
2284  //if(TMath::Sqrt(nSigTPCpion*nSigTPCpion+nSigTOFpion*nSigTOFpion)<=fNSigmaCut){
2285  //fHistTOFnSigmavsPtAfter[0]->Fill(dPt1*gCharge1,nSigTOFpion);
2286  //if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
2287  //fHistPtwithTOFSignal[0]->Fill(dPt1*gCharge1);
2288  //}
2289  //}
2290  //if(TMath::Sqrt(nSigTPCkaon*nSigTPCkaon+nSigTOFkaon*nSigTOFkaon)<=fNSigmaCut){
2291  //fHistTOFnSigmavsPtAfter[1]->Fill(dPt1*gCharge1,nSigTOFkaon);
2292  //if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
2293  //fHistPtwithTOFSignal[1]->Fill(dPt1*gCharge1);
2294  //}
2295  //}
2296  //if(TMath::Sqrt(nSigTPCproton*nSigTPCproton+nSigTOFproton*nSigTOFproton)<=fNSigmaCut){
2297  //fHistTOFnSigmavsPtAfter[2]->Fill(dPt1*gCharge1,nSigTOFproton);
2298  //if(TOFmatch>0 && probMis < 0.01 && beta>0.2){
2299  //fHistPtwithTOFSignal[2]->Fill(dPt1*gCharge1);
2300  //}
2301  //}
2302 
2303 
2304 
2305 
2306  if(TMath::Abs(nSigTPCpion)<=fNSigmaCut){
2307  fHistPtwithTPCNsigma[0]->Fill(dPt1*gCharge1);
2308  fHistTPCdEdxvsPtPIDAfter[0]->Fill(dPt1,dEdx1);
2309  }
2310  if(TMath::Abs(nSigTPCkaon)<=fNSigmaCut){
2311  fHistPtwithTPCNsigma[1]->Fill(dPt1*gCharge1);
2312  fHistTPCdEdxvsPtPIDAfter[1]->Fill(dPt1,dEdx1);
2313  }
2314  if(TMath::Abs(nSigTPCproton)<=fNSigmaCut){
2315  fHistPtwithTPCNsigma[2]->Fill(dPt1*gCharge1);
2316  fHistTPCdEdxvsPtPIDAfter[2]->Fill(dPt1,dEdx1);
2317  }
2318 
2319 
2320 
2321  //-------------- Fill NUA for Charged tracks ----------------
2322  if(gCharge1>0){
2323  fHist3DEtaPhiVz_Pos_Run[3][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2324  }
2325  else if(gCharge1<0){
2326  fHist3DEtaPhiVz_Neg_Run[3][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2327  }
2328 
2329  if(bFillNUAHistPID){
2330  //============== Fill NUA Histograms for Pion ---------------------
2331  if(dPt1<0.6 && TMath::Abs(nSigTPCpion)<=2.5){
2332  if(gCharge1>0){
2333  fHist3DEtaPhiVz_Pos_Run[0][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2334  }
2335  else if(gCharge1<0){
2336  fHist3DEtaPhiVz_Neg_Run[0][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2337  }
2338  }
2339  else if(dPt1>=0.6 && dPt1<=2.0 && TMath::Abs(nSigTPCpion)<=2.5 && TMath::Abs(nSigTOFpion)<=2.0 ){
2340  if(gCharge1>0){
2341  fHist3DEtaPhiVz_Pos_Run[0][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2342  }
2343  else if(gCharge1<0){
2344  fHist3DEtaPhiVz_Neg_Run[0][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2345  }
2346  }
2347 
2348  //============== Fill NUA Histograms for Kaon ---------------------
2349  if(dPt1<0.6 && TMath::Abs(nSigTPCkaon)<=2.5){
2350  if(gCharge1>0){
2351  fHist3DEtaPhiVz_Pos_Run[1][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2352  }
2353  else if(gCharge1<0){
2354  fHist3DEtaPhiVz_Neg_Run[1][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2355  }
2356  }
2357  else if(dPt1>=0.6 && dPt1<=2.0 && TMath::Abs(nSigTPCkaon)<=2.5 && TMath::Abs(nSigTOFkaon)<=2.0){
2358  if(gCharge1>0){
2359  fHist3DEtaPhiVz_Pos_Run[1][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2360  }
2361  else if(gCharge1<0){
2362  fHist3DEtaPhiVz_Neg_Run[1][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2363  }
2364  }
2365 
2366  //============== Fill NUA Histograms for proton ---------------------
2367  if(dPt1<0.8 && TMath::Abs(nSigTPCproton)<=2.5){
2368  if(gCharge1>0){
2369  fHist3DEtaPhiVz_Pos_Run[2][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2370  }
2371  else if(gCharge1<0){
2372  fHist3DEtaPhiVz_Neg_Run[2][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2373  }
2374  }
2375  else if(dPt1>=0.8 && dPt1<=3.5 && TMath::Abs(nSigTPCproton)<=2.5 && TMath::Abs(nSigTOFproton)<=2.5){
2376  if(gCharge1>0){
2377  fHist3DEtaPhiVz_Pos_Run[2][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2378  }
2379  else if(gCharge1<0){
2380  fHist3DEtaPhiVz_Neg_Run[2][cForNUA]->Fill(pVtxZ,dPhi1,dEta1);
2381  }
2382  }
2383 
2384  }// Fill NUA for PID or not?
2385 
2386 
2387 
2388 
2389  //------------------------------------------------------------
2390 
2391  // TOF matching is not needed from data as done by Davide using MC.
2392  /*
2393  if(!TOFmatch || probMis > 0.01 || beta>0.2) continue;
2394 
2395  // nSigmaTPC distribution for Fixed nSigmaTOF cut
2396  if(TMath::Abs(nSigTOFpion)<=fNSigmaCut){
2397  fHistTPCnSigmavsPtAfter[0]->Fill(dPt1*gCharge1,nSigTPCpion);
2398  fHistPtwithTOFmasscut[0]->Fill(dPt1*gCharge1);
2399  }
2400  if(TMath::Abs(nSigTOFkaon)<=fNSigmaCut){
2401  fHistTPCnSigmavsPtAfter[1]->Fill(dPt1*gCharge1,nSigTPCkaon);
2402  fHistPtwithTOFmasscut[1]->Fill(dPt1*gCharge1);
2403  }
2404  if(TMath::Abs(nSigTOFproton)<=fNSigmaCut){
2405  fHistTPCnSigmavsPtAfter[2]->Fill(dPt1*gCharge1,nSigTPCproton);
2406  fHistPtwithTOFmasscut[2]->Fill(dPt1*gCharge1);
2407  } */
2408 
2409 
2410 
2411  }
2412 //===================== track loop ends ============================
2413 
2414 
2415 
2416  //---------- Do event by event business here -------
2417 
2418 
2419 
2420 
2421 
2422  // ----------------- Store Q-Vectors --------------
2423  //V0A-V0C
2424  fHist_Reso2n_EP_Norm_Det[QAindex][0]->Fill(EvtCent, TMath::Cos(gPsiN*(PsiNV0A-PsiNV0C)));
2425  //V0A-TPC
2426  fHist_Reso2n_EP_Norm_Det[QAindex][1]->Fill(EvtCent, TMath::Cos(gPsiN*(PsiNV0A-PsiNTPCF)));
2427  //V0C-TPC
2428  fHist_Reso2n_EP_Norm_Det[QAindex][2]->Fill(EvtCent, TMath::Cos(gPsiN*(PsiNV0C-PsiNTPCF)));
2429  //TPCa -TPCc
2430  fHist_Reso2n_EP_Norm_Det[QAindex][3]->Fill(EvtCent, TMath::Cos(gPsiN*(PsiNTPCA-PsiNTPCC)));
2431 
2432 
2433 
2434 
2435  //--------- Store TPC-Qn for Recenter ---------
2436  sumTPCQn2x[0] = sumTPCQn2x[0]/multEtaNeg;
2437  sumTPCQn2y[0] = sumTPCQn2y[0]/multEtaNeg;
2438  sumTPCQn3x[0] = sumTPCQn3x[0]/multEtaNeg;
2439  sumTPCQn3y[0] = sumTPCQn3y[0]/multEtaNeg;
2440  sumTPCQn4x[0] = sumTPCQn4x[0]/multEtaNeg;
2441  sumTPCQn4y[0] = sumTPCQn4y[0]/multEtaNeg;
2442 
2443  sumTPCQn2x[1] = sumTPCQn2x[1]/multEtaPos;
2444  sumTPCQn2y[1] = sumTPCQn2y[1]/multEtaPos;
2445  sumTPCQn3x[1] = sumTPCQn3x[1]/multEtaPos;
2446  sumTPCQn3y[1] = sumTPCQn3y[1]/multEtaPos;
2447  sumTPCQn4x[1] = sumTPCQn4x[1]/multEtaPos;
2448  sumTPCQn4y[1] = sumTPCQn4y[1]/multEtaPos;
2449 
2450  fTPCAQ2xVsCentRun->Fill(EvtCent,sumTPCQn2x[0]);
2451  fTPCAQ2yVsCentRun->Fill(EvtCent,sumTPCQn2y[0]);
2452  fTPCCQ2xVsCentRun->Fill(EvtCent,sumTPCQn2x[1]);
2453  fTPCCQ2yVsCentRun->Fill(EvtCent,sumTPCQn2y[1]);
2454 
2455  fTPCAQ3xVsCentRun->Fill(EvtCent,sumTPCQn3x[0]);
2456  fTPCAQ3yVsCentRun->Fill(EvtCent,sumTPCQn3y[0]);
2457  fTPCCQ3xVsCentRun->Fill(EvtCent,sumTPCQn3x[1]);
2458  fTPCCQ3yVsCentRun->Fill(EvtCent,sumTPCQn3y[1]);
2459 
2460  fTPCAQ4xVsCentRun->Fill(EvtCent,sumTPCQn4x[0]);
2461  fTPCAQ4yVsCentRun->Fill(EvtCent,sumTPCQn4y[0]);
2462  fTPCCQ4xVsCentRun->Fill(EvtCent,sumTPCQn4x[1]);
2463  fTPCCQ4yVsCentRun->Fill(EvtCent,sumTPCQn4y[1]);
2464 
2465 
2466 
2467 
2468 
2469 
2470  PostData(1,fListHist);
2471 
2472  fHistEventCount->Fill(14.5); //15th bin is last one
2473  stepCount++;
2474 
2475 
2476 
2477 
2478 
2479 
2480  //if(multPOI2nd!=(multPOI1st-1))
2481  //cout<<"mismatched "<< "\tPOIs1st = "<<multPOI1st<<"\tPOIs2nd = "<<multPOI2nd <<" for Event = "<<fEventCount<<endl;
2482 
2483  //if(fEventCount%10==0)
2484  //cout<<"Ev = "<<fEventCount<<"\tMult = "<<multEtaFull<<"\tPOIs1st = "<<multPOI1st<<"\t POI2nd = "<< multPOI2nd <<endl;
2485  //cout<<"Ev = "<<fEventCount<<"\tMult = "<<multEtaFull<<"\tPOIs = "<<multPOI1st<<"\tRealTime = "<< watch.RealTime() <<"\tCPUTime = "<< watch.CpuTime()<<endl;
2486 
2487  //watch.Stop();
2488 
2489 
2490 
2491 
2492 
2493 }//================ UserExec ==============
2494 
2495 
2496 
2497 
2498 
2499 
2500 
2501 
2502 
2503 
2504 
2505 
2506 
2507 
2508 
2509 
2510 
2511 
2512 
2513 
2514 
2515 
2517 
2518 
2519 
2520 
2521 double AliAnalysisTaskCMEV0PID::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
2522 {
2523  // calculate sqrt of weighted distance to other vertex
2524  if (!v0 || !v1) {
2525  AliDebug(2,"\n\n ::GetWDist => One of vertices is not valid\n\n");
2526  return 0;
2527  }
2528  static TMatrixDSym vVb(3);
2529  double dist = -1;
2530  double dx = v0->GetX()-v1->GetX();
2531  double dy = v0->GetY()-v1->GetY();
2532  double dz = v0->GetZ()-v1->GetZ();
2533  double cov0[6],cov1[6];
2534  v0->GetCovarianceMatrix(cov0);
2535  v1->GetCovarianceMatrix(cov1);
2536  vVb(0,0) = cov0[0]+cov1[0];
2537  vVb(1,1) = cov0[2]+cov1[2];
2538  vVb(2,2) = cov0[5]+cov1[5];
2539  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
2540  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
2541  vVb.InvertFast();
2542  if (!vVb.IsValid()) {
2543  AliDebug(2,"Singular Matrix\n");
2544  return dist;
2545  }
2546  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
2547  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
2548  return dist>0 ? TMath::Sqrt(dist) : -1;
2549 }
2550 
2551 
2553  { // check for multi-vertexer pile-up
2554  const int kMinPlpContrib = 5;
2555  const double kMaxPlpChi2 = 5.0;
2556  const double kMinWDist = 15;
2557 
2558  const AliVVertex* vtPrm = 0;
2559  const AliVVertex* vtPlp = 0;
2560 
2561  int nPlp = 0;
2562 
2563  if(!(nPlp=faod->GetNumberOfPileupVerticesTracks()))
2564  return kFALSE;
2565 
2566  vtPrm = faod->GetPrimaryVertex();
2567  if(vtPrm == faod->GetPrimaryVertexSPD())
2568  return kTRUE; // there are pile-up vertices but no primary
2569 
2570  //int bcPrim = vtPrm->GetBC();
2571 
2572  for(int ipl=0;ipl<nPlp;ipl++) {
2573  vtPlp = (const AliVVertex*)faod->GetPileupVertexTracks(ipl);
2574  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
2575  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
2576  //int bcPlp = vtPlp->GetBC();
2577  //if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2)
2578  // return kTRUE; // pile-up from other BC
2579 
2580  double wDst = GetWDist(vtPrm,vtPlp);
2581  if (wDst<kMinWDist) continue;
2582 
2583  return kTRUE; // pile-up: well separated vertices
2584  }
2585  return kFALSE;
2586 }
2587 
2588 
2589 
2591 
2592  if(bApplyMCcorr){
2593  if(!gGrid){
2594  TGrid::Connect("alien://");
2595  }
2596 
2597  if(!mfileFBHijing){
2598 
2599  mfileFBHijing = TFile::Open(sMCfilePath,"READ");
2600 
2601  fListFBHijing = dynamic_cast<TList*>(mfileFBHijing->FindObjectAny("fMcEffiHij"));
2602 
2603  if(!fListFBHijing){
2604  std::cout<<"\n\n !!!!**** Warning: FB Efficiency File/List not found *****\n\n"<<std::endl;
2605  //exit(1);
2606  }
2607  else if(fListFBHijing) {
2608  cout<<"\n =========> Info: FB Efficiency is Used from file = "<<sMCfilePath.Data()<<endl;
2609  for(int i=0;i<10;i++) {
2610  fFB_Efficiency_Cent[i] = (TH1D *) fListFBHijing->FindObject(Form("eff_unbiased_%d",i));
2611  //std::cout<<" input MC hist"<<i<<" = "<<fFB_Efficiency_Cent[i]->GetName()<<std::endl;
2612  }
2613  }
2614  }
2615  }
2616  else{ // if MC efficiency Not used/ file not found, then use weight = 1.
2617  for(int i=0;i<10;i++){
2618  fFB_Efficiency_Cent[i] = new TH1D(Form("eff_unbiased_%d",i),"",1,0,50.);
2619  fFB_Efficiency_Cent[i]->SetBinContent(1,1.0);
2620  }
2621  if(bApplyMCcorr){ printf("\n!!***** !!!! WARNING !!!! *****!!\n MC correction File not found, using Efficiency = 1.0 !!\n\n");}
2622  }
2623 }
2624 
2625 
2626 
2627 
2629 
2630  Int_t cIndex = 0;
2631 
2632  if(fCent<5.0) {
2633  cIndex = 0;
2634  }
2635  else if(fCent>=5.0 && fCent<10){
2636  cIndex = 1;
2637  }
2638  else if(fCent>=10.0) {
2639  cIndex = abs(fCent/10.0)+1;
2640  }
2641  return cIndex;
2642 }
2643 
2645  //std::cout<<" centrality outlier function called "<<std::endl;
2646  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.};
2647 
2648  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.};
2649 
2650  for(int i=0;i<90;i++) {
2651  hCentvsTPCmultCuts->SetBinContent(i+1,1,fMeanTPC[i]);
2652  hCentvsTPCmultCuts->SetBinContent(i+1,2,fSigmaTPC[i]);
2653  }
2654 }
2655 
2656 
2657 
2658 
2659 
2660 
2662 
2663  fHistTaskConfigParameters = new TH1F("fHistTaskConfigParameters","Task parameters",20,0,20);
2664  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(1,"FilterBit");
2665  fHistTaskConfigParameters->SetBinContent(1,fFilterBit);
2666  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(2,"n#sigmaTPC");
2667  fHistTaskConfigParameters->SetBinContent(2,fNSigmaCut);
2668  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(3,"MinPt");
2669  fHistTaskConfigParameters->SetBinContent(3,fMinPtCut);
2670  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(4,"MaxPt");
2671  fHistTaskConfigParameters->SetBinContent(4,fMaxPtCut);
2672  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(5,"MinEta");
2673  fHistTaskConfigParameters->SetBinContent(5,fMinEtaCut);
2674  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(6,"MaxEta");
2675  fHistTaskConfigParameters->SetBinContent(6,fMaxEtaCut);
2676 
2677  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(11,"CentralityMin");
2679  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(12,"CentralityMax");
2681 
2682  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(13,"VertexMin(cm)");
2683  fHistTaskConfigParameters->GetXaxis()->SetBinLabel(14,"VertexMax(cm)");
2684 
2686 
2687 
2688 
2689  fHistPileUpCount = new TH1F("fHistPileUpCount", "fHistPileUpCount", 15, 0., 15.);
2690  fHistPileUpCount->GetXaxis()->SetBinLabel(1,"plpMV");
2691  fHistPileUpCount->GetXaxis()->SetBinLabel(2,"fromSPD");
2692  fHistPileUpCount->GetXaxis()->SetBinLabel(3,"RefMultComb08");
2693  fHistPileUpCount->GetXaxis()->SetBinLabel(4,"IncompleteDAQ");
2694  fHistPileUpCount->GetXaxis()->SetBinLabel(5,"abs(V0M-CL1)>5.0");
2695  fHistPileUpCount->GetXaxis()->SetBinLabel(6,"missingVtx");
2696  fHistPileUpCount->GetXaxis()->SetBinLabel(7,"inconsistentVtx");
2697  Int_t puConst = fPileUpConstParm;
2698  fHistPileUpCount->GetXaxis()->SetBinLabel(8,Form("multESDTPCDif>%d",puConst));
2699  fHistPileUpCount->GetXaxis()->SetBinLabel(9,Form("multGlobTPCDif>%d",puConst));
2700  fHistPileUpCount->GetXaxis()->SetBinLabel(10,"PileUpMultSelTask");
2702 
2703 
2704  fHistMultSelPUCount = new TH1F("fHistMultSelPileUpCount", "no of PU Event from MultSelTask", 10, 0., 10);
2705  fHistMultSelPUCount->GetXaxis()->SetBinLabel(1,"PileUp");
2706  fHistMultSelPUCount->GetXaxis()->SetBinLabel(2,"PileUpMV");
2707  fHistMultSelPUCount->GetXaxis()->SetBinLabel(3,"PileUpMultBins");
2708  fHistMultSelPUCount->GetXaxis()->SetBinLabel(4,"InconsistentVtx");
2709  fHistMultSelPUCount->GetXaxis()->SetBinLabel(5,"TrackletVsCluster");
2710  fHistMultSelPUCount->GetXaxis()->SetBinLabel(6,"IncompleteDAQ");
2711  fHistMultSelPUCount->GetXaxis()->SetBinLabel(7,"NotGoodVertex2016");
2713 
2714 
2715  fHistEventCount = new TH1F("fHistEventCount","Event counts",15,0,15);
2716  fHistEventCount->GetXaxis()->SetBinLabel(1,"Called UserExec()");
2717  fHistEventCount->GetXaxis()->SetBinLabel(2,"Called Exec()");
2718  fHistEventCount->GetXaxis()->SetBinLabel(3,"AOD Exist");
2719  fHistEventCount->GetXaxis()->SetBinLabel(4,"Vz < 10");
2720  fHistEventCount->GetXaxis()->SetBinLabel(5,Form("%2.0f<Cent<%2.0f",fCentralityPercentMin,fCentralityPercentMax));
2721  fHistEventCount->GetXaxis()->SetBinLabel(6,"noAODtrack > 2 ");
2722  fHistEventCount->GetXaxis()->SetBinLabel(7,"TPC vs Global");
2723  fHistEventCount->GetXaxis()->SetBinLabel(8,"TPC128 vs ESD");
2724  fHistEventCount->GetXaxis()->SetBinLabel(9,"Cent vs TPC");
2725  fHistEventCount->GetXaxis()->SetBinLabel(10,"mult eta+/- > 2");
2726  fHistEventCount->GetXaxis()->SetBinLabel(11,"centCL1 < 90");
2727  fHistEventCount->GetXaxis()->SetBinLabel(15,"Survived Events");
2728  fListHist->Add(fHistEventCount);
2729 
2730  //fHistEventCount->Fill(1);
2731 
2732 }
2733 
2734 
2735 
2736 
2737 
2739 {
2740  if(fListV0MCorr){
2741  fHCorrectV0M = (TH1D *) fListV0MCorr->FindObject(Form("fHistV0Gain_Run%d",run));
2742  fHAvgerageQnV0A = (TH2D *) fListV0MCorr->FindObject(Form("fHistAvgQnV0A_Run%d",run));
2743  fHAvgerageQnV0C = (TH2D *) fListV0MCorr->FindObject(Form("fHistAvgQnV0C_Run%d",run));
2744 
2745  if(fHCorrectV0M){
2746  cout<<"\n =========== Info:: Setting up V0 gain correction for run = "<<run<<"============"<<endl;
2747  }
2748  else{
2749  cout<<"\n =========== Info:: No V0 gain correction..!!! for run = "<<run<<"============"<<endl;
2750  }
2751  }
2752  else{
2753  cout<<"\n ======== Error:: List of V0 gain correction not found for run "<<run<<"============"<<endl;
2754  fHCorrectV0M = new TH1D("fHCorrectV0M","",64,0,64);
2755  for(int i=1;i<=64;i++){
2756  fHCorrectV0M->SetBinContent(i,1.0);
2757  }
2758  fHAvgerageQnV0A = new TH2D("fHAvgerageQnV0A_empty","<Cos2>,<Sin2>,<Cos3>,<Sin3> V0A",90,0,90,8,0,8);
2759  fHAvgerageQnV0C = new TH2D("fHAvgerageQnV0C_empty","<Cos2>,<Sin2>,<Cos3>,<Sin3> V0A",90,0,90,8,0,8);
2760  for(int i=1;i<=90;i++){
2761  for(int j=1;j<=8;j++){
2762  fHAvgerageQnV0A->SetBinContent(i,j,0.0);
2763  fHAvgerageQnV0C->SetBinContent(i,j,0.0);
2764  }
2765  }
2766  }
2767 }
2768 
2769 
2770 
2771 
2772 
2773 
2774 
2776 {
2777 
2778  if(fListNUACorr){
2779  for(int i=0;i<5;i++){
2780  fHCorrectNUApos[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Pos_Cent%d_Run%d",i,run));
2781  fHCorrectNUAneg[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Neg_Cent%d_Run%d",i,run));
2782  }
2783  if(fHCorrectNUApos[0] && fHCorrectNUApos[0]){
2784  cout<<"\n=========== Info:: Setting up NUA corrections for run "<<run<<"============"<<endl;
2785  }
2786  }
2787  else {
2788  printf("\n ******** Warning: No NUA Correction for Charge Particle in run %d, Use Wgt = 1.0 ********* \n",run);
2789  for(int i=0;i<5;i++){
2790  fHCorrectNUApos[i] = new TH3D(Form("fHCorrectNUApos_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2791  fHCorrectNUAneg[i] = new TH3D(Form("fHCorrectNUAneg_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2792  fHCorrectNUApos[i]->SetBinContent(1,1,1,1.0);
2793  fHCorrectNUAneg[i]->SetBinContent(1,1,1,1.0);
2794  //exit(1);
2795  }
2796  }
2797 
2798  //=================== PID: ==========================
2799  if(fListNUACorr){
2800  for(int i=0;i<5;i++){
2801  fHCorrectNUAposPion[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Pion_Pos_Cent%d_Run%d",i,run)); //
2802  fHCorrectNUAnegPion[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Pion_Neg_Cent%d_Run%d",i,run));
2803 
2804  fHCorrectNUAposKaon[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Kaon_Pos_Cent%d_Run%d",i,run)); //
2805  fHCorrectNUAnegKaon[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Kaon_Neg_Cent%d_Run%d",i,run));
2806 
2807  fHCorrectNUAposProton[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Proton_Pos_Cent%d_Run%d",i,run));
2808  fHCorrectNUAnegProton[i] = (TH3D *) fListNUACorr->FindObject(Form("fHist_NUA_VzPhiEta_Proton_Neg_Cent%d_Run%d",i,run));
2809  }
2811  cout<<"\n=========== Info:: Setting up --> PID NUA corrections for run = "<<run<<"============"<<endl;
2812  }
2813  else{
2814  cout<<"\n=========== WARNING :: PID NUA corrections NOT found for run = "<<run<<"============"<<endl;
2815  }
2816  }
2817  else {
2818  printf("\n ******** Warning: No NUA Correction for PID for run %d, Use Wgt = 1.0 ********* \n",run);
2819  for(int i=0;i<5;i++){
2820  fHCorrectNUAposPion[i] = new TH3D(Form("fHCorrectNUAPionpos_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2821  fHCorrectNUAposPion[i] = new TH3D(Form("fHCorrectNUAPionneg_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2822  fHCorrectNUAposPion[i]->SetBinContent(1,1,1,1.0);
2823  fHCorrectNUAposPion[i]->SetBinContent(1,1,1,1.0);
2824 
2825  fHCorrectNUAposKaon[i] = new TH3D(Form("fHCorrectNUAKaonpos_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2826  fHCorrectNUAposKaon[i] = new TH3D(Form("fHCorrectNUAKaonneg_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2827  fHCorrectNUAposKaon[i]->SetBinContent(1,1,1,1.0);
2828  fHCorrectNUAposKaon[i]->SetBinContent(1,1,1,1.0);
2829 
2830  fHCorrectNUAposProton[i] = new TH3D(Form("fHCorrectNUAProtonpos_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2831  fHCorrectNUAposProton[i] = new TH3D(Form("fHCorrectNUAProtonneg_cent%d",i),"",1,-10,10,1,0,6.284,1,-0.9,0.9);
2832  fHCorrectNUAposProton[i]->SetBinContent(1,1,1,1.0);
2833  fHCorrectNUAposProton[i]->SetBinContent(1,1,1,1.0);
2834  //exit(1);
2835  }
2836  }
2837 
2838 }
2839 
Int_t GetCentralityScaled0to10(Double_t fCent)
double Double_t
Definition: External.C:58
Definition: External.C:260
TProfile * fV0AQ2xVsCentRun
V0C Average <Qn>, n=2,3.
TH2D * fHAvgerageQnV0A
To read Gain Correction file.
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
TProfile * fHist_Corr3p_Kaon_EP_Norm_PP[2][4]
AliAnalysisUtils * fAnalysisUtil
TProfile * fHist_Corr3p_Proton_EP_Norm_NN[2][4]
Double_t mass
AliMultSelection * fMultSelection
PID response Handler.
centrality
TProfile * fHist_Corr3p_Kaon_EP_Norm_PN[2][4]
char Char_t
Definition: External.C:18
TH1F * fHistPtwithTPCNsigma[3]
last in the list
TProfile * fHist_Reso2n_EP_Norm_Det[2][4]
TCanvas * c
Definition: TestFitELoss.C:172
TH2F * fHistEtaPtAfter
Eta-Pt acceptance.
TProfile * fHist_Corr3p_EP_Norm_NN[2][4]
TH3D * fHCorrectNUAnegProton[5]
5 centrality bin, read NUA from file
TH3D * fHCorrectNUAposPion[5]
5 centrality bin, read NUA from file
TH1F * fCentDistAfter
without PileUp cut
TProfile * fHist_Corr3p_Pion_EP_Norm_PP[2][4]
TH3D * fHCorrectNUAnegPion[5]
5 centrality bin, read NUA from file
AliPIDResponse * fPIDResponse
aod
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Definition: External.C:252
TH3F * fHist3DEtaPhiVz_Neg_Run[4][5]
4 particle 5 centrality bin
TProfile * fHist_Corr3p_Proton_EP_Norm_PP[2][4]
Definition: External.C:228
Definition: External.C:212
TH1D * fHCorrectV0M
with PileUp cut
TH2F * fHistTPCvsGlobalMultBefore
Eta-Pt acceptance.
TH3D * fHCorrectNUAposKaon[5]
5 centrality bin, read NUA from file
TProfile * fHist_Corr3p_Pion_EP_Norm_PN[2][4]
TProfile * fHist_Corr3p_Kaon_EP_Norm_NN[2][4]
TProfile * fHist_Corr3p_Pion_EP_Norm_NN[2][4]
TProfile * fHist_Corr3p_Proton_EP_Norm_PN[2][4]
TH1F * fCentDistBefore
To fill VOM multiplicity.
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 * fHCorrectNUAnegKaon[5]
5 centrality bin, read NUA from file
TH3D * fHCorrectNUAneg[5]
5 centrality bin, read NUA from file
double GetWDist(const AliVVertex *v0, const AliVVertex *v1)
TH1D * fFB_Efficiency_Cent[10]
4 particle 5 centrality bin
TProfile * fHist_Corr3p_EP_Norm_PN[2][4]
5 centrality bin, read NUA from file
const char Option_t
Definition: External.C:48
TH3D * fHCorrectNUAposProton[5]
5 centrality bin, read NUA from file
bool Bool_t
Definition: External.C:53
void SetupMCcorrectionMap(TString sMCfilePath)
TH2D * fHAvgerageQnV0C
V0A Average <Qn>, n=2,3.
TProfile * fHist_Corr3p_EP_Norm_PP[2][4]