update from Marco
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskV2AllChAOD.cxx
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
48 using namespace std;
49
50 ClassImp(AliAnalysisTaskV2AllChAOD)
51
52 //________________________________________________________________________
53 AliAnalysisTaskV2AllChAOD::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),
71   fTrkBit(272),
72   fEtaCut(0.8),
73   fMinPt(0),
74   fMaxPt(20.0),
75   fMinTPCNcls(70),
76   fFillTHn(kFALSE),
77   fResSP(0),
78   fResSP_vs_Cent(0),
79   f2partCumQA_vs_Cent(0),
80   f2partCumQB_vs_Cent(0),
81   fEta_vs_Phi_bef(0),
82   fEta_vs_PhiA(0),
83   fEta_vs_PhiB(0),
84   fResSP_lq(0),
85   fResSP_vs_Cent_lq(0),
86   f2partCumQA_vs_Cent_lq(0),
87   f2partCumQB_vs_Cent_lq(0),
88   fResSP_sq(0),
89   fResSP_vs_Cent_sq(0),
90   f2partCumQA_vs_Cent_sq(0),
91   f2partCumQB_vs_Cent_sq(0)
92 {
93   
94   for (Int_t i = 0; i< 9; i++){
95     
96     fResSP_vs_Qvec[i] = 0;
97     
98     fv2SPGap1A[i] = 0;
99     fv2SPGap1B[i] = 0;
100
101     fSinGap1Aq[i] = 0;
102     fCosGap1Aq[i] = 0;
103     fSinGap1Bq[i] = 0;
104     fCosGap1Bq[i] = 0;
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;
114     
115     fSinGap1Aq_lq[i] = 0;
116     fCosGap1Aq_lq[i] = 0;
117     fSinGap1Bq_lq[i] = 0;
118     fCosGap1Bq_lq[i] = 0;
119     
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;
128     
129     fSinGap1Aq_sq[i] = 0;
130     fCosGap1Aq_sq[i] = 0;
131     fSinGap1Bq_sq[i] = 0;
132     fCosGap1Bq_sq[i] = 0;
133     
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 //________________________________________________________________________
151 void 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   
169   if( fFillTHn ){ 
170     //dimensions of THnSparse for Q vector checks
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.};
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");
185     NSparseHistEv->GetAxis(4)->SetTitle("resolution");
186     NSparseHistEv->GetAxis(4)->SetName("res");
187     fOutput->Add(NSparseHistEv);
188   }
189   
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   
197   fResSP = new TProfile("fResSP", "Resolution; centrality; Resolution", 9, -0.5, 8.5);
198   fOutput->Add(fResSP);
199   
200   fResSP_vs_Cent = new TProfile("fResSP_vs_Cent", "Resolution; centrality; Resolution", 20., 0., 100.);
201   fOutput->Add(fResSP_vs_Cent);
202   
203   f2partCumQA_vs_Cent = new TProfile("f2partCumQA_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
204   fOutput->Add(f2partCumQA_vs_Cent);
205   
206   f2partCumQB_vs_Cent = new TProfile("f2partCumQB_vs_Cent", "Resolution; centrality; Resolution", 100., 0., 100.);
207   fOutput->Add(f2partCumQB_vs_Cent);
208
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.);
210   fOutput->Add(fEta_vs_Phi_bef);
211   
212   fEta_vs_PhiA = new TH2D("fEta_vs_PhiA","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
213   fOutput->Add(fEta_vs_PhiA);
214   
215   fEta_vs_PhiB = new TH2D("fEta_vs_PhiB","eta vs phi distribution;#eta;#phi",200.,-1.,1.,175.,0.,7.);
216   fOutput->Add(fEta_vs_PhiB);
217   
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   
222   fResSP_vs_Cent_lq = new TProfile("fResSP_vs_Cent_lq", "Resolution; centrality; Resolution", 20., 0., 100.);
223   fOutput_lq->Add(fResSP_vs_Cent_lq);
224   
225   f2partCumQA_vs_Cent_lq = new TProfile("f2partCumQA_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
226   fOutput_lq->Add(f2partCumQA_vs_Cent_lq);
227   
228   f2partCumQB_vs_Cent_lq = new TProfile("f2partCumQB_vs_Cent_lq", "Resolution; centrality; Resolution", 100., 0., 100.);
229   fOutput_lq->Add(f2partCumQB_vs_Cent_lq);
230   
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);
234
235   fResSP_vs_Cent_sq = new TProfile("fResSP_vs_Cent_sq", "Resolution; centrality; Resolution", 20., 0., 100.);
236   fOutput_sq->Add(fResSP_vs_Cent_sq);
237   
238   f2partCumQA_vs_Cent_sq = new TProfile("f2partCumQA_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
239   fOutput_sq->Add(f2partCumQA_vs_Cent_sq);
240   
241   f2partCumQB_vs_Cent_sq = new TProfile("f2partCumQB_vs_Cent_sq", "Resolution; centrality; Resolution", 100., 0., 100.);
242   fOutput_sq->Add(f2partCumQB_vs_Cent_sq);
243   
244   for (Int_t iC = 0; iC < 9; iC++){
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]);
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
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
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]);
257       
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
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
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
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
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
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
351 void 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   
379   Int_t centV0 = -1;
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;
403   
404   for (Int_t loop = 0; loop < 2; loop++){
405
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)
411     
412       fEta_vs_Phi_bef->Fill( track->Eta(), track->Phi() );
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++;
420
421           fSinGap1Aq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
422           fCosGap1Aq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
423           
424           fEta_vs_PhiA->Fill( track->Eta(), track->Phi() );
425                     
426           if (Qvec > fCutLargeQperc && Qvec < 100.){
427             fSinGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
428             fCosGap1Aq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
429           }
430       
431           if (Qvec > 0. && Qvec < fCutSmallQperc){
432             fSinGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
433             fCosGap1Aq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
434           }
435         }
436     
437       if (track->Eta() < fEtaGapMin){
438         QxGap1B += TMath::Cos(2.*track->Phi());
439         QyGap1B += TMath::Sin(2.*track->Phi());
440         multGap1B++;
441                     
442         fCosGap1Bq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
443         fSinGap1Bq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
444         
445         fEta_vs_PhiB->Fill( track->Eta(), track->Phi() );
446                     
447         if (Qvec > fCutLargeQperc && Qvec < 100.){
448           fSinGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
449           fCosGap1Bq_lq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
450         }
451       
452         if (Qvec > 0. && Qvec < fCutSmallQperc){
453           fSinGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
454           fCosGap1Bq_sq[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
455         }
456       }
457   
458     } else {
459       
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);
464
465           fSinGap1A[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
466           fCosGap1A[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
467       
468           if (Qvec > fCutLargeQperc && Qvec < 100.){
469             fv2SPGap1A_lq[centV0]->Fill(track->Pt(), v2SPGap1A);
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           }
473       
474           if (Qvec > 0. && Qvec < fCutSmallQperc){
475             fv2SPGap1A_sq[centV0]->Fill(track->Pt(), v2SPGap1A);
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           }
479         }
480       
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);
484           
485           fCosGap1B[centV0]->Fill(track->Pt(), TMath::Cos(2.*track->Phi()));
486           fSinGap1B[centV0]->Fill(track->Pt(), TMath::Sin(2.*track->Phi()));
487       
488           if (Qvec > fCutLargeQperc && Qvec < 100.){
489             fv2SPGap1B_lq[centV0]->Fill(track->Pt(), v2SPGap1B);
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           }
493       
494           if (Qvec > 0. && Qvec < fCutSmallQperc){
495             fv2SPGap1B_sq[centV0]->Fill(track->Pt(), v2SPGap1B);
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           }
499         }
500       }// end else 
501     } // end loop on tracks
502   } // end loop
503   
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);
508     fResSP_vs_Cent->Fill(Cent, res);
509     fResSP_vs_Qvec[centV0]->Fill(Qvec,res);
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)) );
514     if(f2partCumQA>0)f2partCumQA_vs_Cent->Fill((Double_t)Cent,f2partCumQA);
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)) );
519     if(f2partCumQB>0)f2partCumQB_vs_Cent->Fill((Double_t)Cent,f2partCumQB);
520         
521     if (Qvec > fCutLargeQperc && Qvec < 100.){
522       fResSP_lq->Fill((Double_t)centV0, res);
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);
526     }
527     
528     if (Qvec > 0. && Qvec < fCutSmallQperc){
529       fResSP_sq->Fill((Double_t)centV0, res);
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);
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   
539       Double_t varEv[5];
540       varEv[0]=Cent;
541       varEv[1]=Qvec;
542       varEv[2]=(Double_t)QA;
543       varEv[2]=(Double_t)QB;
544       varEv[4]=(Double_t)res;
545       ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
546     }
547   }
548   
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 //_________________________________________________________________
557 Bool_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 //_________________________________________________________________
579 void   AliAnalysisTaskV2AllChAOD::Terminate(Option_t *)
580 {
581   // Terminate
582 }