1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------
17 // AliAnalysisTaskV2AllChAOD class
18 //-----------------------------------------------------------------
26 #include "THnSparse.h"
29 #include "AliAnalysisTask.h"
30 #include "AliAODTrack.h"
31 #include "AliAODMCParticle.h"
32 #include "AliVParticle.h"
33 #include "AliAODEvent.h"
34 #include "AliAODInputHandler.h"
35 #include "AliAnalysisTaskV2AllChAOD.h"
36 #include "AliAnalysisTaskESDfilter.h"
37 #include "AliAnalysisDataContainer.h"
38 #include "AliSpectraAODTrackCuts.h"
39 #include "AliSpectraAODEventCuts.h"
40 #include "AliPIDCombined.h"
41 #include "AliCentrality.h"
43 #include "AliVEvent.h"
45 #include <TMCProcess.h>
52 ClassImp(AliAnalysisTaskV2AllChAOD)
54 //________________________________________________________________________
55 AliAnalysisTaskV2AllChAOD::AliAnalysisTaskV2AllChAOD(const char *name) : AliAnalysisTaskSE(name),
84 f2partCumQA_vs_Cent(0),
85 f2partCumQB_vs_Cent(0),
91 f2partCumQA_vs_Cent_lq(0),
92 f2partCumQB_vs_Cent_lq(0),
95 f2partCumQA_vs_Cent_sq(0),
96 f2partCumQB_vs_Cent_sq(0),
98 fv2SPGap1A_inclusive_mb(0),
99 fv2SPGap1B_inclusive_mb(0),
100 fv2SPGap1A_inclusive_lq(0),
101 fv2SPGap1B_inclusive_lq(0),
102 fv2SPGap1A_inclusive_sq(0),
103 fv2SPGap1B_inclusive_sq(0),
104 fResSPmc_inclusive(0),
105 fv2SPGap1Amc_inclusive_mb(0),
106 fv2SPGap1Bmc_inclusive_mb(0),
107 fv2SPGap1Amc_inclusive_lq(0),
108 fv2SPGap1Bmc_inclusive_lq(0),
109 fv2SPGap1Amc_inclusive_sq(0),
110 fv2SPGap1Bmc_inclusive_sq(0),
116 fDoCentrSystCentrality(0)
119 for (Int_t i = 0; i< 9; i++){
121 fResSP_vs_Qvec[i] = 0;
137 fv2SPGap1A_lq[i] = 0;
138 fv2SPGap1B_lq[i] = 0;
140 fSinGap1Aq_lq[i] = 0;
141 fCosGap1Aq_lq[i] = 0;
142 fSinGap1Bq_lq[i] = 0;
143 fCosGap1Bq_lq[i] = 0;
151 fv2SPGap1A_sq[i] = 0;
152 fv2SPGap1B_sq[i] = 0;
154 fSinGap1Aq_sq[i] = 0;
155 fCosGap1Aq_sq[i] = 0;
156 fSinGap1Bq_sq[i] = 0;
157 fCosGap1Bq_sq[i] = 0;
166 fRecoEffList=new TList();
167 fRecoEffList->SetOwner();
168 fRecoEffList->SetName("fRecoEffList");
170 // Default constructor
171 DefineInput(0, TChain::Class());
172 DefineOutput(1, TList::Class());
173 DefineOutput(2, AliSpectraAODEventCuts::Class());
174 DefineOutput(3, AliSpectraAODTrackCuts::Class());
175 DefineOutput(4, TList::Class());
176 DefineOutput(5, TList::Class());
179 //________________________________________________________________________
180 void AliAnalysisTaskV2AllChAOD::UserCreateOutputObjects()
182 // create output objects
185 fOutput->SetName("fOutput");
187 fOutput_lq=new TList();
188 fOutput_lq->SetOwner();
189 fOutput_lq->SetName("fOutput_lq");
191 fOutput_sq=new TList();
192 fOutput_sq->SetOwner();
193 fOutput_sq->SetName("fOutput_sq");
195 if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
196 if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
199 //dimensions of THnSparse for Q vector checks
200 const Int_t nvarev=6;
201 // cent q-rec_perc qvec-rec q-gen_tracks qvec-gen_vzero Nch
202 Int_t binsHistRealEv[nvarev] = { fnCentBins, 100, fnQvecBins, fnQvecBins, fnQvecBins, fnNchBins};
203 Double_t xminHistRealEv[nvarev] = { 0., 0., 0., 0., 0., 0.};
204 Double_t xmaxHistRealEv[nvarev] = { 100., 100., fQvecUpperLim, fQvecUpperLim, fQvecUpperLim, 2000.};
206 THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
207 NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
208 NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
210 NSparseHistEv->GetAxis(1)->SetTitle("q-vec rec percentile");
211 NSparseHistEv->GetAxis(1)->SetName("Qrec_perc");
213 NSparseHistEv->GetAxis(2)->SetTitle("q-vec rec");
214 NSparseHistEv->GetAxis(2)->SetName("Qrec");
216 NSparseHistEv->GetAxis(3)->SetTitle("q-vec gen tracks");
217 NSparseHistEv->GetAxis(3)->SetName("Qgen_tracks");
219 NSparseHistEv->GetAxis(4)->SetTitle("q-vec gen vzero");
220 NSparseHistEv->GetAxis(4)->SetName("Qgen_vzero");
222 NSparseHistEv->GetAxis(5)->SetTitle("Ncharged");
223 NSparseHistEv->GetAxis(5)->SetName("Nch");
224 fOutput->Add(NSparseHistEv);
227 fCentrality = new TH1D("fCentrality", "centrality distribution; centrality", 200, 0., 100);
228 fOutput->Add(fCentrality);
230 fQvector = new TH1D("fQvector", "q-vector distribution; q-vector", fnQvecBins, 0., fQvecUpperLim);
231 fOutput->Add(fQvector);
233 fQvector_lq = new TH1D("fQvector_lq", "q-vector distribution; q-vector", fnQvecBins, 0., fQvecUpperLim);
234 fOutput_lq->Add(fQvector_lq);
236 fQvector_sq = new TH1D("fQvector_sq", "q-vector distribution; q-vector", fnQvecBins, 0., fQvecUpperLim);
237 fOutput_sq->Add(fQvector_sq);
239 // binning common to all the THn
240 //change it according to your needs + move it to global variables -> setter/getter
241 // Double_t ptBins[] = {0., 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 20.0};
242 // const Int_t nptBins = 31;
243 Double_t ptBins[] = {0., 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 3.0, 3.4, 3.8, 4.2, 4.6, 5.0, 5.5, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 20.0};
244 const Int_t nptBins = 33;
246 fResSP = new TProfile("fResSP", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
247 fOutput->Add(fResSP);
249 fResSP_vs_Cent = new TProfile("fResSP_vs_Cent", "Resolution; centrality; Resolution", 20., 0., 100.);
250 fOutput->Add(fResSP_vs_Cent);
252 f2partCumQA_vs_Cent = new TProfile("f2partCumQA_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
253 fOutput->Add(f2partCumQA_vs_Cent);
255 f2partCumQB_vs_Cent = new TProfile("f2partCumQB_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
256 fOutput->Add(f2partCumQB_vs_Cent);
258 fEta_vs_Phi_bef = new TH2D("fEta_vs_Phi_bef","eta vs phi distribution before eta gap;#eta;#phi",200.,-1.,1.,175.,0.,7.);
259 fOutput->Add(fEta_vs_Phi_bef);
261 fEta_vs_PhiA = new TH2D("fEta_vs_PhiA","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
262 fOutput->Add(fEta_vs_PhiA);
264 fEta_vs_PhiB = new TH2D("fEta_vs_PhiB","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
265 fOutput->Add(fEta_vs_PhiB);
268 fResSP_inclusive = new TProfile("fResSP_inclusive", "Resolution; ese; Resolution", 3, 0., 3.);
269 fOutput->Add(fResSP_inclusive);
271 fv2SPGap1A_inclusive_mb = new TProfile("fv2SPGap1A_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
272 fOutput->Add(fv2SPGap1A_inclusive_mb);
274 fv2SPGap1B_inclusive_mb = new TProfile("fv2SPGap1B_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
275 fOutput->Add(fv2SPGap1B_inclusive_mb);
279 fResSP_lq = new TProfile("fResSP_lq", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
280 fOutput_lq->Add(fResSP_lq);
282 fResSP_vs_Cent_lq = new TProfile("fResSP_vs_Cent_lq", "Resolution; centrality; Resolution", 20., 0., 100.);
283 fOutput_lq->Add(fResSP_vs_Cent_lq);
285 f2partCumQA_vs_Cent_lq = new TProfile("f2partCumQA_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
286 fOutput_lq->Add(f2partCumQA_vs_Cent_lq);
288 f2partCumQB_vs_Cent_lq = new TProfile("f2partCumQB_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
289 fOutput_lq->Add(f2partCumQB_vs_Cent_lq);
292 fv2SPGap1A_inclusive_lq = new TProfile("fv2SPGap1A_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
293 fOutput_lq->Add(fv2SPGap1A_inclusive_lq);
295 fv2SPGap1B_inclusive_lq = new TProfile("fv2SPGap1B_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
296 fOutput_lq->Add(fv2SPGap1B_inclusive_lq);
299 fResSP_sq = new TProfile("fResSP_sq", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
300 fOutput_sq->Add(fResSP_sq);
302 fResSP_vs_Cent_sq = new TProfile("fResSP_vs_Cent_sq", "Resolution; centrality; Resolution", 20., 0., 100.);
303 fOutput_sq->Add(fResSP_vs_Cent_sq);
305 f2partCumQA_vs_Cent_sq = new TProfile("f2partCumQA_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
306 fOutput_sq->Add(f2partCumQA_vs_Cent_sq);
308 f2partCumQB_vs_Cent_sq = new TProfile("f2partCumQB_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
309 fOutput_sq->Add(f2partCumQB_vs_Cent_sq);
312 fv2SPGap1A_inclusive_sq = new TProfile("fv2SPGap1A_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
313 fOutput_sq->Add(fv2SPGap1A_inclusive_sq);
315 fv2SPGap1B_inclusive_sq = new TProfile("fv2SPGap1B_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
316 fOutput_sq->Add(fv2SPGap1B_inclusive_sq);
318 for (Int_t iC = 0; iC < 9; iC++){
320 fResSP_vs_Qvec[iC] = new TProfile(Form("fResSP_vs_Qvec_%d", iC), "Resolution; Qvec (V0A); Resolution", 10., 0., 100.);
321 fOutput->Add(fResSP_vs_Qvec[iC]);
323 fv2SPGap1A[iC] = new TProfile(Form("fv2SPGap1A_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
324 fOutput->Add(fv2SPGap1A[iC]);
326 fv2SPGap1B[iC] = new TProfile(Form("fv2SPGap1B_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
327 fOutput->Add(fv2SPGap1B[iC]);
329 fSinGap1Aq[iC] = new TProfile(Form("fSinGap1Aq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
330 fOutput->Add(fSinGap1Aq[iC]);
332 fCosGap1Aq[iC] = new TProfile(Form("fCosGap1Aq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
333 fOutput->Add(fCosGap1Aq[iC]);
335 fSinGap1Bq[iC] = new TProfile(Form("fSinGap1Bq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
336 fOutput->Add(fSinGap1Bq[iC]);
338 fCosGap1Bq[iC] = new TProfile(Form("fCosGap1Bq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
339 fOutput->Add(fCosGap1Bq[iC]);
341 fSinGap1A[iC] = new TProfile(Form("fSinGap1A_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
342 fOutput->Add(fSinGap1A[iC]);
344 fCosGap1A[iC] = new TProfile(Form("fCosGap1A_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
345 fOutput->Add(fCosGap1A[iC]);
347 fSinGap1B[iC] = new TProfile(Form("fSinGap1B_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
348 fOutput->Add(fSinGap1B[iC]);
350 fCosGap1B[iC] = new TProfile(Form("fCosGap1B_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
351 fOutput->Add(fCosGap1B[iC]);
354 fv2SPGap1A_lq[iC] = new TProfile(Form("fv2SPGap1A_lq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
355 fOutput_lq->Add(fv2SPGap1A_lq[iC]);
357 fv2SPGap1B_lq[iC] = new TProfile(Form("fv2SPGap1B_lq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
358 fOutput_lq->Add(fv2SPGap1B_lq[iC]);
360 fSinGap1Aq_lq[iC] = new TProfile(Form("fSinGap1Aq_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
361 fOutput_lq->Add(fSinGap1Aq_lq[iC]);
363 fCosGap1Aq_lq[iC] = new TProfile(Form("fCosGap1Aq_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
364 fOutput_lq->Add(fCosGap1Aq_lq[iC]);
366 fSinGap1Bq_lq[iC] = new TProfile(Form("fSinGap1Bq_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
367 fOutput_lq->Add(fSinGap1Bq_lq[iC]);
369 fCosGap1Bq_lq[iC] = new TProfile(Form("fCosGap1Bq_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
370 fOutput_lq->Add(fCosGap1Bq_lq[iC]);
372 fSinGap1A_lq[iC] = new TProfile(Form("fSinGap1A_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
373 fOutput_lq->Add(fSinGap1A_lq[iC]);
375 fCosGap1A_lq[iC] = new TProfile(Form("fCosGap1A_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
376 fOutput_lq->Add(fCosGap1A_lq[iC]);
378 fSinGap1B_lq[iC] = new TProfile(Form("fSinGap1B_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
379 fOutput_lq->Add(fSinGap1B_lq[iC]);
381 fCosGap1B_lq[iC] = new TProfile(Form("fCosGap1B_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
382 fOutput_lq->Add(fCosGap1B_lq[iC]);
385 fv2SPGap1A_sq[iC] = new TProfile(Form("fv2SPGap1A_sq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
386 fOutput_sq->Add(fv2SPGap1A_sq[iC]);
388 fv2SPGap1B_sq[iC] = new TProfile(Form("fv2SPGap1B_sq_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
389 fOutput_sq->Add(fv2SPGap1B_sq[iC]);
391 fSinGap1Aq_sq[iC] = new TProfile(Form("fSinGap1Aq_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
392 fOutput_sq->Add(fSinGap1Aq_sq[iC]);
394 fCosGap1Aq_sq[iC] = new TProfile(Form("fCosGap1Aq_sq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
395 fOutput_sq->Add(fCosGap1Aq_sq[iC]);
397 fSinGap1Bq_sq[iC] = new TProfile(Form("fSinGap1Bq_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
398 fOutput_sq->Add(fSinGap1Bq_sq[iC]);
400 fCosGap1Bq_sq[iC] = new TProfile(Form("fCosGap1Bq_sq_%d", iC), "p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
401 fOutput_sq->Add(fCosGap1Bq_sq[iC]);
403 fSinGap1A_sq[iC] = new TProfile(Form("fSinGap1A_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
404 fOutput_sq->Add(fSinGap1A_sq[iC]);
406 fCosGap1A_sq[iC] = new TProfile(Form("fCosGap1A_sq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
407 fOutput_sq->Add(fCosGap1A_sq[iC]);
409 fSinGap1B_sq[iC] = new TProfile(Form("fSinGap1B_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
410 fOutput_sq->Add(fSinGap1B_sq[iC]);
412 fCosGap1B_sq[iC] = new TProfile(Form("fCosGap1B_sq_%d", iC), "p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
413 fOutput_sq->Add(fCosGap1B_sq[iC]);
417 fResSPmc_inclusive = new TProfile("fResSPmc_inclusive", "Resolution; ese; Resolution", 3, 0., 3.);
418 fOutput->Add(fResSPmc_inclusive);
420 fv2SPGap1Amc_inclusive_mb = new TProfile("fv2SPGap1Amc_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
421 fOutput->Add(fv2SPGap1Amc_inclusive_mb);
423 fv2SPGap1Bmc_inclusive_mb = new TProfile("fv2SPGap1Bmc_inclusive_mb", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
424 fOutput->Add(fv2SPGap1Bmc_inclusive_mb);
427 fv2SPGap1Amc_inclusive_lq = new TProfile("fv2SPGap1Amc_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
428 fOutput_lq->Add(fv2SPGap1Amc_inclusive_lq);
430 fv2SPGap1Bmc_inclusive_lq = new TProfile("fv2SPGap1Bmc_inclusive_lq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
431 fOutput_lq->Add(fv2SPGap1Bmc_inclusive_lq);
434 fv2SPGap1Amc_inclusive_sq = new TProfile("fv2SPGap1Amc_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
435 fOutput_sq->Add(fv2SPGap1Amc_inclusive_sq);
437 fv2SPGap1Bmc_inclusive_sq = new TProfile("fv2SPGap1Bmc_inclusive_sq", "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
438 fOutput_sq->Add(fv2SPGap1Bmc_inclusive_sq);
442 PostData(1, fOutput );
443 PostData(2, fEventCuts);
444 PostData(3, fTrackCuts);
445 PostData(4, fOutput_lq );
446 PostData(5, fOutput_sq );
449 //________________________________________________________________________
451 void AliAnalysisTaskV2AllChAOD::UserExec(Option_t *)
453 //Printf("An event");
455 fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
457 AliWarning("ERROR: AliAODEvent not available \n");
461 if (strcmp(fAOD->ClassName(), "AliAODEvent"))
463 AliFatal("Not processing AODs");
466 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
468 //Get q-vector percentile.
470 if(fIsMC && fQvecGen) Qvec = fEventCuts->GetQvecPercentileMC(fVZEROside, fQgenType);
471 else Qvec = fEventCuts->GetQvecPercentile(fVZEROside);
473 fQvector->Fill(Qvec);
474 if (Qvec > fCutLargeQperc && Qvec < 100.) fQvector_lq->Fill(Qvec);
475 if (Qvec > 0. && Qvec < fCutSmallQperc) fQvector_sq->Fill(Qvec);
477 Double_t Cent=(fDoCentrSystCentrality)?1.01*fEventCuts->GetCent():fEventCuts->GetCent();
478 fCentrality->Fill(Cent);
481 if ((Cent > 0) && (Cent <= 5.0))
483 else if ((Cent > 5.0) && (Cent <= 10.0))
485 else if ((Cent > 10.0) && (Cent <= 20.0))
487 else if ((Cent > 20.0) && (Cent <= 30.0))
489 else if ((Cent > 30.0) && (Cent <= 40.0))
491 else if ((Cent > 40.0) && (Cent <= 50.0))
493 else if ((Cent > 50.0) && (Cent <= 60.0))
495 else if ((Cent > 60.0) && (Cent <= 70.0))
497 else if ((Cent > 70.0) && (Cent <= 80.0))
500 if(fIsMC) MCclosure(Qvec); // fill mc histograms for montecarlo closure
502 Double_t QxGap1A = 0., QyGap1A = 0.;
503 Double_t QxGap1B = 0., QyGap1B = 0.;
504 Int_t multGap1A = 0, multGap1B = 0;
506 for (Int_t loop = 0; loop < 2; loop++){
508 //main loop on tracks
509 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
510 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
511 if(!track) AliFatal("Not a standard AOD");
512 if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
513 if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
515 fEta_vs_Phi_bef->Fill( track->Eta(), track->Phi() );
519 // 2) reject randomly tracks at high pT until the reconstruction efficiency becomes flat (add the following before the loop == 0 part): (mail by Alexandru)
521 Double_t recoEff = GetRecoEff(track->Pt(), centV0);
523 cout<<"No reconstruction efficiency!"<<endl;
527 Double_t rndPt = gRandom->Rndm();
528 // cout<<"rndPt: "<<rndPt<<endl;
530 Double_t minRecPt = GetRecoEff(0.200001, centV0);
531 // cout<<"minRecPt: "<<minRecPt<<endl;
533 if (rndPt > minRecPt/recoEff){
534 // cout<<"Track rejected: "<<iTracks<<" from "<<fAOD->GetNumberOfTracks()<<endl;
542 if (track->Eta() > fEtaGapMax){
543 QxGap1A += TMath::Cos(2.*track->Phi());
544 QyGap1A += TMath::Sin(2.*track->Phi());
547 fSinGap1Aq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
548 fCosGap1Aq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
550 fEta_vs_PhiA->Fill( track->Eta(), track->Phi() );
552 if (Qvec > fCutLargeQperc && Qvec < 100.){
553 fSinGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
554 fCosGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
557 if (Qvec > 0. && Qvec < fCutSmallQperc){
558 fSinGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
559 fCosGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
564 if (track->Eta() < fEtaGapMin){
565 QxGap1B += TMath::Cos(2.*track->Phi());
566 QyGap1B += TMath::Sin(2.*track->Phi());
569 fCosGap1Bq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
570 fSinGap1Bq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
572 fEta_vs_PhiB->Fill( track->Eta(), track->Phi() );
574 if (Qvec > fCutLargeQperc && Qvec < 100.){
575 fSinGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
576 fCosGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
579 if (Qvec > 0. && Qvec < fCutSmallQperc){
580 fSinGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
581 fCosGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
588 //eval v2 scalar product
589 if (track->Eta() < fEtaGapMin && multGap1A > 0){
590 Double_t v2SPGap1A = (TMath::Cos(2.*track->Phi())*QxGap1A + TMath::Sin(2.*track->Phi())*QyGap1A)/(Double_t)multGap1A;
591 fv2SPGap1A[centV0]->Fill(track->Pt(), v2SPGap1A);
593 fv2SPGap1A_inclusive_mb->Fill(track->Pt(), v2SPGap1A); //mb v2 for mc closure
595 fSinGap1A[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
596 fCosGap1A[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
598 if (Qvec > fCutLargeQperc && Qvec < 100.){
599 fv2SPGap1A_lq[centV0]->Fill(track->Pt(), v2SPGap1A);
600 fSinGap1A_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
601 fCosGap1A_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
603 fv2SPGap1A_inclusive_lq->Fill(track->Pt(), v2SPGap1A); //lq v2 for mc closure
606 if (Qvec > 0. && Qvec < fCutSmallQperc){
607 fv2SPGap1A_sq[centV0]->Fill(track->Pt(), v2SPGap1A);
608 fSinGap1A_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
609 fCosGap1A_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
611 fv2SPGap1A_inclusive_sq->Fill(track->Pt(), v2SPGap1A); //sq v2 for mc closure
616 if (track->Eta() > fEtaGapMax && multGap1B > 0){
617 Double_t v2SPGap1B = (TMath::Cos(2.*track->Phi())*QxGap1B + TMath::Sin(2.*track->Phi())*QyGap1B)/(Double_t)multGap1B;
618 fv2SPGap1B[centV0]->Fill(track->Pt(), v2SPGap1B);
620 fv2SPGap1B_inclusive_mb->Fill(track->Pt(), v2SPGap1B); //mb v2 for mc closure
622 fCosGap1B[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
623 fSinGap1B[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
625 if (Qvec > fCutLargeQperc && Qvec < 100.){
626 fv2SPGap1B_lq[centV0]->Fill(track->Pt(), v2SPGap1B);
627 fSinGap1B_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
628 fCosGap1B_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
630 fv2SPGap1B_inclusive_lq->Fill(track->Pt(), v2SPGap1B); //lq v2 for mc closure
633 if (Qvec > 0. && Qvec < fCutSmallQperc){
634 fv2SPGap1B_sq[centV0]->Fill(track->Pt(), v2SPGap1B);
635 fSinGap1B_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
636 fCosGap1B_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
638 fv2SPGap1B_inclusive_sq->Fill(track->Pt(), v2SPGap1B); //sq v2 for mc closure
643 } // end loop on tracks
647 if (multGap1A > 0 && multGap1B > 0){
648 Double_t res = (QxGap1A*QxGap1B + QyGap1A*QyGap1B)/(Double_t)multGap1A/(Double_t)multGap1B;
649 fResSP->Fill((Double_t)centV0, res);
651 fResSP_inclusive->Fill(0., res); //mb v2 for mc closure
652 fResSP_vs_Cent->Fill(Cent, res);
653 fResSP_vs_Qvec[centV0]->Fill(Qvec,res);
655 Double_t f2partCumQA = -999.;
657 f2partCumQA = ( ( (QxGap1A*QxGap1A + QyGap1A*QyGap1A) - (Double_t)multGap1A ) / ((Double_t)multGap1A*((Double_t)multGap1A-1)) );
658 if(f2partCumQA>0)f2partCumQA_vs_Cent->Fill((Double_t)Cent,f2partCumQA);
660 Double_t f2partCumQB = -999.;
662 f2partCumQB = ( ( (QxGap1B*QxGap1B + QyGap1B*QyGap1B) - (Double_t)multGap1B ) / ((Double_t)multGap1B*((Double_t)multGap1B-1)) );
663 if(f2partCumQB>0)f2partCumQB_vs_Cent->Fill((Double_t)Cent,f2partCumQB);
665 if (Qvec > fCutLargeQperc && Qvec < 100.){
666 fResSP_lq->Fill((Double_t)centV0, res);
667 fResSP_vs_Cent_lq->Fill(Cent, res);
668 if(f2partCumQA>0)f2partCumQA_vs_Cent_lq->Fill((Double_t)Cent,f2partCumQA);
669 if(f2partCumQB>0)f2partCumQB_vs_Cent_lq->Fill((Double_t)Cent,f2partCumQB);
671 fResSP_inclusive->Fill(1., res); //lq v2 for mc closure
674 if (Qvec > 0. && Qvec < fCutSmallQperc){
675 fResSP_sq->Fill((Double_t)centV0, res);
676 fResSP_vs_Cent_sq->Fill(Cent, res);
677 if(f2partCumQA>0)f2partCumQA_vs_Cent_sq->Fill((Double_t)Cent,f2partCumQA);
678 if(f2partCumQB>0)f2partCumQB_vs_Cent_sq->Fill((Double_t)Cent,f2partCumQB);
680 fResSP_inclusive->Fill(2., res); //sq v2 for mc closure
682 }// end multiplicity if
689 varEv[1]=(Double_t)Qvec; // qvec_rec_perc
691 Double_t qvzero = 0.;
692 if(fVZEROside==0)qvzero=(Double_t)fEventCuts->GetqV0A();
693 else if (fVZEROside==1)qvzero=(Double_t)fEventCuts->GetqV0C(); // qvec_rec
694 varEv[2]=(Double_t)qvzero; // qvec from VZERO
696 Double_t qgen_tracks = (Double_t)fEventCuts->CalculateQVectorMC(fVZEROside, 0);
697 varEv[3]= (Double_t)qgen_tracks;
699 Double_t qgen_vzero = (Double_t)fEventCuts->CalculateQVectorMC(fVZEROside, 1);
700 varEv[4]= (Double_t)qgen_vzero;
702 varEv[5]=(Double_t)fEventCuts->GetNch(); // Nch
704 ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
709 PostData(1, fOutput );
710 PostData(2, fEventCuts);
711 PostData(3, fTrackCuts);
712 PostData(4, fOutput_lq );
713 PostData(5, fOutput_sq );
716 //_________________________________________________________________
717 Bool_t AliAnalysisTaskV2AllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
719 //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
720 //FIXME should update EventCuts?
721 //FIXME add track->GetXYZ(p) method
723 double xyz[3],cov[3];
725 if (!trk->GetXYZ(xyz)) { // dca is not stored
726 AliExternalTrackParam etp;
727 etp.CopyFromVTrack(trk);
728 AliVEvent* ev = (AliVEvent*)trk->GetEvent();
729 if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
730 if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
738 //_________________________________________________________________
739 void AliAnalysisTaskV2AllChAOD::MCclosure(Double_t qvec){
740 // First do MC to fill up the MC particle array
742 TClonesArray *arrayMC = 0;
745 arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
747 AliFatal("Error: MC particles branch not found!\n");
750 Double_t QxGap1Amc = 0., QyGap1Amc = 0.;
751 Double_t QxGap1Bmc = 0., QyGap1Bmc = 0.;
752 Int_t multGap1Amc = 0, multGap1Bmc = 0;
754 for (Int_t loop = 0; loop < 2; loop++){
756 Int_t nMC = arrayMC->GetEntries();
758 for (Int_t iMC = 0; iMC < nMC; iMC++)
760 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
761 if(!partMC->Charge()) continue;//Skip neutrals
762 if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
764 if (!(partMC->IsPhysicalPrimary()))
767 if(partMC->Eta()<fTrackCuts->GetEtaMin() || partMC->Eta()>fTrackCuts->GetEtaMax()) continue;
769 //Printf("a particle");
774 if (partMC->Eta() > fEtaGapMax){
775 QxGap1Amc += TMath::Cos(2.*partMC->Phi());
776 QyGap1Amc += TMath::Sin(2.*partMC->Phi());
780 if (partMC->Eta() < fEtaGapMin){
781 QxGap1Bmc += TMath::Cos(2.*partMC->Phi());
782 QyGap1Bmc += TMath::Sin(2.*partMC->Phi());
788 //eval v2 scalar product
789 if (partMC->Eta() < fEtaGapMin && multGap1Amc > 0){
790 Double_t v2SPGap1Amc = (TMath::Cos(2.*partMC->Phi())*QxGap1Amc + TMath::Sin(2.*partMC->Phi())*QyGap1Amc)/(Double_t)multGap1Amc;
791 fv2SPGap1Amc_inclusive_mb->Fill(partMC->Pt(), v2SPGap1Amc);
793 if (qvec > fCutLargeQperc && qvec < 100.){
794 fv2SPGap1Amc_inclusive_lq->Fill(partMC->Pt(), v2SPGap1Amc);
797 if (qvec > 0. && qvec < fCutSmallQperc){
798 fv2SPGap1Amc_inclusive_sq->Fill(partMC->Pt(), v2SPGap1Amc);
803 if (partMC->Eta() > fEtaGapMax && multGap1Bmc > 0){
804 Double_t v2SPGap1Bmc = (TMath::Cos(2.*partMC->Phi())*QxGap1Bmc + TMath::Sin(2.*partMC->Phi())*QyGap1Bmc)/(Double_t)multGap1Bmc;
805 fv2SPGap1Bmc_inclusive_mb->Fill(partMC->Pt(), v2SPGap1Bmc);
807 if (qvec > fCutLargeQperc && qvec < 100.){
808 fv2SPGap1Bmc_inclusive_lq->Fill(partMC->Pt(), v2SPGap1Bmc);
811 if (qvec > 0. && qvec < fCutSmallQperc){
812 fv2SPGap1Bmc_inclusive_sq->Fill(partMC->Pt(), v2SPGap1Bmc);
819 } // end loop on partMCs
822 if (multGap1Amc > 0 && multGap1Bmc > 0){
823 Double_t resmc = (QxGap1Amc*QxGap1Bmc + QyGap1Amc*QyGap1Bmc)/(Double_t)multGap1Amc/(Double_t)multGap1Bmc;
824 fResSPmc_inclusive->Fill(0.,resmc);
826 if (qvec > fCutLargeQperc && qvec < 100.){
827 fResSPmc_inclusive->Fill(1.,resmc);
830 if (qvec > 0. && qvec < fCutSmallQperc){
831 fResSPmc_inclusive->Fill(2.,resmc);
839 //_________________________________________________________________
840 Double_t AliAnalysisTaskV2AllChAOD::GetRecoEff(Double_t pt, Int_t iC){
844 if(pt<0.2 || pt>100.) return 1.;
846 // // // //spectra ese binning
847 // // // //const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.,20.,25.,30.,35.,40.,50.,75.,100.};
848 // // // //const Int_t nptBins=34;
850 const Double_t fEpsilon=0.000001;
852 TH1F *h = (TH1F*)fRecoEffList->At(iC);
854 Int_t bin = h->FindBin(pt);
856 Double_t lowlim = h->GetBinLowEdge(bin);
857 Double_t uplim = h->GetBinLowEdge(bin) + h->GetBinWidth(bin);
859 if( pt>lowlim && pt<uplim ) return h->GetBinContent(bin);
860 if( pt == lowlim ) return h->GetBinContent( h->FindBin(pt+fEpsilon) );
861 if( pt == uplim ) return h->GetBinContent( h->FindBin(pt-fEpsilon) );
866 //_________________________________________________________________
867 void AliAnalysisTaskV2AllChAOD::Terminate(Option_t *)