update from Marco
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskV2AllChAOD.cxx
CommitLineData
10a99a07 1/**************************************************************************
2 * Copyright(c) 1998-2009, 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//-----------------------------------------------------------------
17// AliAnalysisTaskV2AllChAOD class
18//-----------------------------------------------------------------
19
20#include "TChain.h"
21#include "TTree.h"
22#include "TLegend.h"
23#include "TH1F.h"
24#include "TH2F.h"
25#include "THnSparse.h"
26#include "TProfile.h"
27#include "TCanvas.h"
28#include "AliAnalysisTask.h"
29#include "AliAODTrack.h"
30#include "AliAODMCParticle.h"
31#include "AliVParticle.h"
32#include "AliAODEvent.h"
33#include "AliAODInputHandler.h"
34#include "AliAnalysisTaskV2AllChAOD.h"
35#include "AliAnalysisTaskESDfilter.h"
36#include "AliAnalysisDataContainer.h"
37#include "AliSpectraAODTrackCuts.h"
38#include "AliSpectraAODEventCuts.h"
39#include "AliPIDCombined.h"
40#include "AliCentrality.h"
41#include "TProof.h"
42#include "AliVEvent.h"
43#include "AliStack.h"
44#include <TMCProcess.h>
45
46#include <iostream>
47
48using namespace std;
49
50ClassImp(AliAnalysisTaskV2AllChAOD)
51
52//________________________________________________________________________
53AliAnalysisTaskV2AllChAOD::AliAnalysisTaskV2AllChAOD(const char *name) : AliAnalysisTaskSE(name),
54 fAOD(0x0),
55 fTrackCuts(0x0),
56 fEventCuts(0x0),
57 fIsMC(0),
58 fCharge(0),
59 fVZEROside(0),
60 fOutput(0x0),
61 fOutput_lq(0x0),
62 fOutput_sq(0x0),
63 fnCentBins(20),
64 fnQvecBins(100),
65 fIsQvecCalibMode(0),
66 fQvecUpperLim(100),
67 fCutLargeQperc(90.),
68 fCutSmallQperc(10.),
69 fEtaGapMin(-0.5),
70 fEtaGapMax(0.5),
d64e71aa 71 fTrkBit(272),
72 fEtaCut(0.8),
73 fMinPt(0),
74 fMaxPt(20.0),
75 fMinTPCNcls(70),
30e20fad 76 fFillTHn(kFALSE),
10a99a07 77 fResSP(0),
af960b8a 78 fResSP_vs_Cent(0),
79 f2partCumQA_vs_Cent(0),
80 f2partCumQB_vs_Cent(0),
bc76879c 81 fEta_vs_Phi_bef(0),
2ba0e068 82 fEta_vs_PhiA(0),
83 fEta_vs_PhiB(0),
10a99a07 84 fResSP_lq(0),
af960b8a 85 fResSP_vs_Cent_lq(0),
86 f2partCumQA_vs_Cent_lq(0),
87 f2partCumQB_vs_Cent_lq(0),
88 fResSP_sq(0),
af960b8a 89 fResSP_vs_Cent_sq(0),
90 f2partCumQA_vs_Cent_sq(0),
91 f2partCumQB_vs_Cent_sq(0)
10a99a07 92{
93
94 for (Int_t i = 0; i< 9; i++){
30e20fad 95
96 fResSP_vs_Qvec[i] = 0;
97
10a99a07 98 fv2SPGap1A[i] = 0;
10a99a07 99 fv2SPGap1B[i] = 0;
bc76879c 100
101 fSinGap1Aq[i] = 0;
102 fCosGap1Aq[i] = 0;
103 fSinGap1Bq[i] = 0;
104 fCosGap1Bq[i] = 0;
10a99a07 105
106 fSinGap1A[i] = 0;
107 fCosGap1A[i] = 0;
108 fSinGap1B[i] = 0;
109 fCosGap1B[i] = 0;
110
111 //large q
112 fv2SPGap1A_lq[i] = 0;
113 fv2SPGap1B_lq[i] = 0;
bc76879c 114
115 fSinGap1Aq_lq[i] = 0;
116 fCosGap1Aq_lq[i] = 0;
117 fSinGap1Bq_lq[i] = 0;
118 fCosGap1Bq_lq[i] = 0;
119
10a99a07 120 fSinGap1A_lq[i] = 0;
121 fCosGap1A_lq[i] = 0;
122 fSinGap1B_lq[i] = 0;
123 fCosGap1B_lq[i] = 0;
124
125 //small q
126 fv2SPGap1A_sq[i] = 0;
127 fv2SPGap1B_sq[i] = 0;
bc76879c 128
129 fSinGap1Aq_sq[i] = 0;
130 fCosGap1Aq_sq[i] = 0;
131 fSinGap1Bq_sq[i] = 0;
132 fCosGap1Bq_sq[i] = 0;
133
10a99a07 134 fSinGap1A_sq[i] = 0;
135 fCosGap1A_sq[i] = 0;
136 fSinGap1B_sq[i] = 0;
137 fCosGap1B_sq[i] = 0;
138
139 }
140
141 // Default constructor
142 DefineInput(0, TChain::Class());
143 DefineOutput(1, TList::Class());
144 DefineOutput(2, AliSpectraAODEventCuts::Class());
145 DefineOutput(3, AliSpectraAODTrackCuts::Class());
146 DefineOutput(4, TList::Class());
147 DefineOutput(5, TList::Class());
148}
149
150//________________________________________________________________________
151void AliAnalysisTaskV2AllChAOD::UserCreateOutputObjects()
152{
153 // create output objects
154 fOutput=new TList();
155 fOutput->SetOwner();
156 fOutput->SetName("fOutput");
157
158 fOutput_lq=new TList();
159 fOutput_lq->SetOwner();
160 fOutput_lq->SetName("fOutput_lq");
161
162 fOutput_sq=new TList();
163 fOutput_sq->SetOwner();
164 fOutput_sq->SetName("fOutput_sq");
165
166 if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
167 if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
168
962f3d11 169 if( fFillTHn ){
170 //dimensions of THnSparse for Q vector checks
af960b8a 171 const Int_t nvarev=5;
172 // cent Q vec Qa Qb res
173 Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, 100, 100, 100};
174 Double_t xminHistRealEv[nvarev] = { 0., 0., 0., 0., 0.};
175 Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 10., 10., 1.};
962f3d11 176 THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
177 NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
178 NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
179 NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
180 NSparseHistEv->GetAxis(1)->SetName("Q_vec");
181 NSparseHistEv->GetAxis(2)->SetTitle("Q_vec (A)");
182 NSparseHistEv->GetAxis(2)->SetName("Qvec_A");
183 NSparseHistEv->GetAxis(3)->SetTitle("Q_vec (B)");
184 NSparseHistEv->GetAxis(3)->SetName("Qvec_B");
af960b8a 185 NSparseHistEv->GetAxis(4)->SetTitle("resolution");
186 NSparseHistEv->GetAxis(4)->SetName("res");
962f3d11 187 fOutput->Add(NSparseHistEv);
188 }
189
10a99a07 190 // binning common to all the THn
191 //change it according to your needs + move it to global variables -> setter/getter
192// 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};
193// const Int_t nptBins = 31;
194 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};
195 const Int_t nptBins = 33;
196
10a99a07 197 fResSP = new TProfile("fResSP", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
198 fOutput->Add(fResSP);
199
30e20fad 200 fResSP_vs_Cent = new TProfile("fResSP_vs_Cent", "Resolution; centrality; Resolution", 20., 0., 100.);
af960b8a 201 fOutput->Add(fResSP_vs_Cent);
202
30e20fad 203 f2partCumQA_vs_Cent = new TProfile("f2partCumQA_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
af960b8a 204 fOutput->Add(f2partCumQA_vs_Cent);
205
30e20fad 206 f2partCumQB_vs_Cent = new TProfile("f2partCumQB_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
af960b8a 207 fOutput->Add(f2partCumQB_vs_Cent);
208
30e20fad 209 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.);
bc76879c 210 fOutput->Add(fEta_vs_Phi_bef);
d64e71aa 211
30e20fad 212 fEta_vs_PhiA = new TH2D("fEta_vs_PhiA","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
2ba0e068 213 fOutput->Add(fEta_vs_PhiA);
214
30e20fad 215 fEta_vs_PhiB = new TH2D("fEta_vs_PhiB","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
2ba0e068 216 fOutput->Add(fEta_vs_PhiB);
d64e71aa 217
10a99a07 218 //large q resolution
219 fResSP_lq = new TProfile("fResSP_lq", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
220 fOutput_lq->Add(fResSP_lq);
221
30e20fad 222 fResSP_vs_Cent_lq = new TProfile("fResSP_vs_Cent_lq", "Resolution; centrality; Resolution", 20., 0., 100.);
af960b8a 223 fOutput_lq->Add(fResSP_vs_Cent_lq);
224
30e20fad 225 f2partCumQA_vs_Cent_lq = new TProfile("f2partCumQA_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
af960b8a 226 fOutput_lq->Add(f2partCumQA_vs_Cent_lq);
227
30e20fad 228 f2partCumQB_vs_Cent_lq = new TProfile("f2partCumQB_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
af960b8a 229 fOutput_lq->Add(f2partCumQB_vs_Cent_lq);
230
10a99a07 231 //small q resolution
232 fResSP_sq = new TProfile("fResSP_sq", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
233 fOutput_sq->Add(fResSP_sq);
af960b8a 234
30e20fad 235 fResSP_vs_Cent_sq = new TProfile("fResSP_vs_Cent_sq", "Resolution; centrality; Resolution", 20., 0., 100.);
af960b8a 236 fOutput_sq->Add(fResSP_vs_Cent_sq);
237
30e20fad 238 f2partCumQA_vs_Cent_sq = new TProfile("f2partCumQA_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
af960b8a 239 fOutput_sq->Add(f2partCumQA_vs_Cent_sq);
240
30e20fad 241 f2partCumQB_vs_Cent_sq = new TProfile("f2partCumQB_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
af960b8a 242 fOutput_sq->Add(f2partCumQB_vs_Cent_sq);
d64e71aa 243
10a99a07 244 for (Int_t iC = 0; iC < 9; iC++){
30e20fad 245
246 fResSP_vs_Qvec[iC] = new TProfile(Form("fResSP_vs_Qvec_%d", iC), "Resolution; Qvec (V0A); Resolution", 10., 0., 100.);
247 fOutput->Add(fResSP_vs_Qvec[iC]);
10a99a07 248
249 fv2SPGap1A[iC] = new TProfile(Form("fv2SPGap1A_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
250 fOutput->Add(fv2SPGap1A[iC]);
251
10a99a07 252 fv2SPGap1B[iC] = new TProfile(Form("fv2SPGap1B_%d", iC), "v_{2}{2} vs p_{T}; p_{T} (GeV/c); v_{2}{2}", nptBins, ptBins);
253 fOutput->Add(fv2SPGap1B[iC]);
254
bc76879c 255 fSinGap1Aq[iC] = new TProfile(Form("fSinGap1Aq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
256 fOutput->Add(fSinGap1Aq[iC]);
10a99a07 257
bc76879c 258 fCosGap1Aq[iC] = new TProfile(Form("fCosGap1Aq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
259 fOutput->Add(fCosGap1Aq[iC]);
260
261 fSinGap1Bq[iC] = new TProfile(Form("fSinGap1Bq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
262 fOutput->Add(fSinGap1Bq[iC]);
263
264 fCosGap1Bq[iC] = new TProfile(Form("fCosGap1Bq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
265 fOutput->Add(fCosGap1Bq[iC]);
266
10a99a07 267 fSinGap1A[iC] = new TProfile(Form("fSinGap1A_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
268 fOutput->Add(fSinGap1A[iC]);
269
270 fCosGap1A[iC] = new TProfile(Form("fCosGap1A_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
271 fOutput->Add(fCosGap1A[iC]);
272
273 fSinGap1B[iC] = new TProfile(Form("fSinGap1B_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
274 fOutput->Add(fSinGap1B[iC]);
275
276 fCosGap1B[iC] = new TProfile(Form("fCosGap1B_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
277 fOutput->Add(fCosGap1B[iC]);
278
279 //large q
280 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);
281 fOutput_lq->Add(fv2SPGap1A_lq[iC]);
282
283 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);
284 fOutput_lq->Add(fv2SPGap1B_lq[iC]);
285
bc76879c 286 fSinGap1Aq_lq[iC] = new TProfile(Form("fSinGap1Aq_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
287 fOutput_lq->Add(fSinGap1Aq_lq[iC]);
288
289 fCosGap1Aq_lq[iC] = new TProfile(Form("fCosGap1Aq_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
290 fOutput_lq->Add(fCosGap1Aq_lq[iC]);
291
292 fSinGap1Bq_lq[iC] = new TProfile(Form("fSinGap1Bq_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
293 fOutput_lq->Add(fSinGap1Bq_lq[iC]);
294
295 fCosGap1Bq_lq[iC] = new TProfile(Form("fCosGap1Bq_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
296 fOutput_lq->Add(fCosGap1Bq_lq[iC]);
297
10a99a07 298 fSinGap1A_lq[iC] = new TProfile(Form("fSinGap1A_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
299 fOutput_lq->Add(fSinGap1A_lq[iC]);
300
301 fCosGap1A_lq[iC] = new TProfile(Form("fCosGap1A_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
302 fOutput_lq->Add(fCosGap1A_lq[iC]);
303
304 fSinGap1B_lq[iC] = new TProfile(Form("fSinGap1B_lq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
305 fOutput_lq->Add(fSinGap1B_lq[iC]);
306
307 fCosGap1B_lq[iC] = new TProfile(Form("fCosGap1B_lq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
308 fOutput_lq->Add(fCosGap1B_lq[iC]);
309
310 //small q
311 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);
312 fOutput_sq->Add(fv2SPGap1A_sq[iC]);
313
314 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);
315 fOutput_sq->Add(fv2SPGap1B_sq[iC]);
316
bc76879c 317 fSinGap1Aq_sq[iC] = new TProfile(Form("fSinGap1Aq_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
318 fOutput_sq->Add(fSinGap1Aq_sq[iC]);
319
320 fCosGap1Aq_sq[iC] = new TProfile(Form("fCosGap1Aq_sq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
321 fOutput_sq->Add(fCosGap1Aq_sq[iC]);
322
323 fSinGap1Bq_sq[iC] = new TProfile(Form("fSinGap1Bq_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
324 fOutput_sq->Add(fSinGap1Bq_sq[iC]);
325
326 fCosGap1Bq_sq[iC] = new TProfile(Form("fCosGap1Bq_sq_%d", iC), "p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
327 fOutput_sq->Add(fCosGap1Bq_sq[iC]);
328
10a99a07 329 fSinGap1A_sq[iC] = new TProfile(Form("fSinGap1A_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
330 fOutput_sq->Add(fSinGap1A_sq[iC]);
331
332 fCosGap1A_sq[iC] = new TProfile(Form("fCosGap1A_sq_%d", iC), ";p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
333 fOutput_sq->Add(fCosGap1A_sq[iC]);
334
335 fSinGap1B_sq[iC] = new TProfile(Form("fSinGap1B_sq_%d", iC), ";p_{T} (GeV/c);#LT sin(2*#phi) #GT", nptBins, ptBins);
336 fOutput_sq->Add(fSinGap1B_sq[iC]);
337
338 fCosGap1B_sq[iC] = new TProfile(Form("fCosGap1B_sq_%d", iC), "p_{T} (GeV/c);#LT cos(2*#phi) #GT", nptBins, ptBins);
339 fOutput_sq->Add(fCosGap1B_sq[iC]);
340 };
341
342 PostData(1, fOutput );
343 PostData(2, fEventCuts);
344 PostData(3, fTrackCuts);
345 PostData(4, fOutput_lq );
346 PostData(5, fOutput_sq );
347}
348
349//________________________________________________________________________
350
351void AliAnalysisTaskV2AllChAOD::UserExec(Option_t *)
352{
353 //Printf("An event");
354 // main event loop
355 fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
356 if (!fAOD) {
357 AliWarning("ERROR: AliAODEvent not available \n");
358 return;
359 }
360
361 if (strcmp(fAOD->ClassName(), "AliAODEvent"))
362 {
363 AliFatal("Not processing AODs");
364 }
365
366 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
367
368 Double_t Qvec=0.;//in case of MC we save space in the memory
369 if(!fIsMC){
370 if(fIsQvecCalibMode){
371 if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
372 else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
373 }
374 else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
375 }
376
377 Double_t Cent=fEventCuts->GetCent();
378
30e20fad 379 Int_t centV0 = -1;
10a99a07 380 if ((Cent > 0) && (Cent <= 5.0))
381 centV0 = 0;
382 else if ((Cent > 5.0) && (Cent <= 10.0))
383 centV0 = 1;
384 else if ((Cent > 10.0) && (Cent <= 20.0))
385 centV0 = 2;
386 else if ((Cent > 20.0) && (Cent <= 30.0))
387 centV0 = 3;
388 else if ((Cent > 30.0) && (Cent <= 40.0))
389 centV0 = 4;
390 else if ((Cent > 40.0) && (Cent <= 50.0))
391 centV0 = 5;
392 else if ((Cent > 50.0) && (Cent <= 60.0))
393 centV0 = 6;
394 else if ((Cent > 60.0) && (Cent <= 70.0))
395 centV0 = 7;
396 else if ((Cent > 70.0) && (Cent <= 80.0))
397 centV0 = 8;
398
399
400 Double_t QxGap1A = 0., QyGap1A = 0.;
401 Double_t QxGap1B = 0., QyGap1B = 0.;
402 Int_t multGap1A = 0, multGap1B = 0;
d57dd83f 403
404 for (Int_t loop = 0; loop < 2; loop++){
10a99a07 405
d57dd83f 406 //main loop on tracks
407 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
408 AliAODTrack* track = fAOD->GetTrack(iTracks);
409 if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
410 if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
10a99a07 411
bc76879c 412 fEta_vs_Phi_bef->Fill( track->Eta(), track->Phi() );
d57dd83f 413
414 if (loop == 0) {
415
416 if (track->Eta() > fEtaGapMax){
417 QxGap1A += TMath::Cos(2.*track->Phi());
418 QyGap1A += TMath::Sin(2.*track->Phi());
419 multGap1A++;
10a99a07 420
bc76879c 421 fSinGap1Aq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
422 fCosGap1Aq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
423
2ba0e068 424 fEta_vs_PhiA->Fill( track->Eta(), track->Phi() );
10a99a07 425
962f3d11 426 if (Qvec > fCutLargeQperc && Qvec < 100.){
bc76879c 427 fSinGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
428 fCosGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
d57dd83f 429 }
10a99a07 430
962f3d11 431 if (Qvec > 0. && Qvec < fCutSmallQperc){
bc76879c 432 fSinGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
433 fCosGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
d57dd83f 434 }
435 }
436
437 if (track->Eta() < fEtaGapMin){
438 QxGap1B += TMath::Cos(2.*track->Phi());
439 QyGap1B += TMath::Sin(2.*track->Phi());
440 multGap1B++;
10a99a07 441
bc76879c 442 fCosGap1Bq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
443 fSinGap1Bq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
444
2ba0e068 445 fEta_vs_PhiB->Fill( track->Eta(), track->Phi() );
10a99a07 446
962f3d11 447 if (Qvec > fCutLargeQperc && Qvec < 100.){
bc76879c 448 fSinGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
449 fCosGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
d57dd83f 450 }
10a99a07 451
962f3d11 452 if (Qvec > 0. && Qvec < fCutSmallQperc){
bc76879c 453 fSinGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
454 fCosGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
d57dd83f 455 }
10a99a07 456 }
10a99a07 457
d57dd83f 458 } else {
10a99a07 459
d57dd83f 460 //eval v2 scalar product
461 if (track->Eta() < fEtaGapMin && multGap1A > 0){
462 Double_t v2SPGap1A = (TMath::Cos(2.*track->Phi())*QxGap1A + TMath::Sin(2.*track->Phi())*QyGap1A)/(Double_t)multGap1A;
463 fv2SPGap1A[centV0]->Fill(track->Pt(), v2SPGap1A);
bc76879c 464
465 fSinGap1A[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
466 fCosGap1A[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
10a99a07 467
962f3d11 468 if (Qvec > fCutLargeQperc && Qvec < 100.){
d57dd83f 469 fv2SPGap1A_lq[centV0]->Fill(track->Pt(), v2SPGap1A);
bc76879c 470 fSinGap1A_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
471 fCosGap1A_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
472 }
10a99a07 473
962f3d11 474 if (Qvec > 0. && Qvec < fCutSmallQperc){
d57dd83f 475 fv2SPGap1A_sq[centV0]->Fill(track->Pt(), v2SPGap1A);
bc76879c 476 fSinGap1A_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
477 fCosGap1A_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
478 }
d57dd83f 479 }
10a99a07 480
d57dd83f 481 if (track->Eta() > fEtaGapMax && multGap1B > 0){
482 Double_t v2SPGap1B = (TMath::Cos(2.*track->Phi())*QxGap1B + TMath::Sin(2.*track->Phi())*QyGap1B)/(Double_t)multGap1B;
483 fv2SPGap1B[centV0]->Fill(track->Pt(), v2SPGap1B);
bc76879c 484
485 fCosGap1B[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
486 fSinGap1B[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
10a99a07 487
962f3d11 488 if (Qvec > fCutLargeQperc && Qvec < 100.){
d57dd83f 489 fv2SPGap1B_lq[centV0]->Fill(track->Pt(), v2SPGap1B);
bc76879c 490 fSinGap1B_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
491 fCosGap1B_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
492 }
d57dd83f 493
962f3d11 494 if (Qvec > 0. && Qvec < fCutSmallQperc){
d57dd83f 495 fv2SPGap1B_sq[centV0]->Fill(track->Pt(), v2SPGap1B);
bc76879c 496 fSinGap1B_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
497 fCosGap1B_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
498 }
d57dd83f 499 }
500 }// end else
501 } // end loop on tracks
502 } // end loop
503
10a99a07 504
505 if (multGap1A > 0 && multGap1B > 0){
506 Double_t res = (QxGap1A*QxGap1B + QyGap1A*QyGap1B)/(Double_t)multGap1A/(Double_t)multGap1B;
507 fResSP->Fill((Double_t)centV0, res);
30e20fad 508 fResSP_vs_Cent->Fill(Cent, res);
509 fResSP_vs_Qvec[centV0]->Fill(Qvec,res);
af960b8a 510
511 Double_t f2partCumQA = -999.;
512 if(multGap1A>1)
513 f2partCumQA = ( ( (QxGap1A*QxGap1A + QyGap1A*QyGap1A) - (Double_t)multGap1A ) / ((Double_t)multGap1A*((Double_t)multGap1A-1)) );
30e20fad 514 if(f2partCumQA>0)f2partCumQA_vs_Cent->Fill((Double_t)Cent,f2partCumQA);
af960b8a 515
516 Double_t f2partCumQB = -999.;
517 if(multGap1B>1)
518 f2partCumQB = ( ( (QxGap1B*QxGap1B + QyGap1B*QyGap1B) - (Double_t)multGap1B ) / ((Double_t)multGap1B*((Double_t)multGap1B-1)) );
30e20fad 519 if(f2partCumQB>0)f2partCumQB_vs_Cent->Fill((Double_t)Cent,f2partCumQB);
10a99a07 520
962f3d11 521 if (Qvec > fCutLargeQperc && Qvec < 100.){
10a99a07 522 fResSP_lq->Fill((Double_t)centV0, res);
30e20fad 523 fResSP_vs_Cent_lq->Fill(Cent, res);
524 if(f2partCumQA>0)f2partCumQA_vs_Cent_lq->Fill((Double_t)Cent,f2partCumQA);
525 if(f2partCumQB>0)f2partCumQB_vs_Cent_lq->Fill((Double_t)Cent,f2partCumQB);
962f3d11 526 }
527
528 if (Qvec > 0. && Qvec < fCutSmallQperc){
10a99a07 529 fResSP_sq->Fill((Double_t)centV0, res);
30e20fad 530 fResSP_vs_Cent_sq->Fill(Cent, res);
531 if(f2partCumQA>0)f2partCumQA_vs_Cent_sq->Fill((Double_t)Cent,f2partCumQA);
532 if(f2partCumQB>0)f2partCumQB_vs_Cent_sq->Fill((Double_t)Cent,f2partCumQB);
962f3d11 533 }
534
535 if( fFillTHn ){
536 Double_t QA = TMath::Sqrt( (QxGap1A*QxGap1A + QyGap1A*QyGap1A)/multGap1A );
537 Double_t QB = TMath::Sqrt( (QxGap1B*QxGap1B + QyGap1B*QyGap1B)/multGap1B );
538
af960b8a 539 Double_t varEv[5];
540 varEv[0]=Cent;
962f3d11 541 varEv[1]=Qvec;
542 varEv[2]=(Double_t)QA;
af960b8a 543 varEv[2]=(Double_t)QB;
544 varEv[4]=(Double_t)res;
962f3d11 545 ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
546 }
10a99a07 547 }
548
10a99a07 549 PostData(1, fOutput );
550 PostData(2, fEventCuts);
551 PostData(3, fTrackCuts);
552 PostData(4, fOutput_lq );
553 PostData(5, fOutput_sq );
554}
555
556//_________________________________________________________________
557Bool_t AliAnalysisTaskV2AllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
558
559 //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
560 //FIXME should update EventCuts?
561 //FIXME add track->GetXYZ(p) method
562
563 double xyz[3],cov[3];
564
565 if (!trk->GetXYZ(xyz)) { // dca is not stored
566 AliExternalTrackParam etp;
567 etp.CopyFromVTrack(trk);
568 AliVEvent* ev = (AliVEvent*)trk->GetEvent();
569 if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
570 if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
571 }
572 p[0] = xyz[0];
573 p[1] = xyz[1];
574 return kTRUE;
575
576}
577
578//_________________________________________________________________
579void AliAnalysisTaskV2AllChAOD::Terminate(Option_t *)
580{
581 // Terminate
582}