Fix a typo
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Hypernuclei / AliAnalysisTaskNucleiv2SP.cxx
1 /**************************************************************************
2  * Contributors are not mentioned at all.                                 *
3  *                                                                        *
4  * Permission to use, copy, modify and distribute this software and its   *
5  * documentation strictly for non-commercial purposes is hereby granted   *
6  * without fee, provided that the above copyright noticxse appears in all *
7  * copies and that both the copyright notice and this permission notice   *
8  * appear in the supporting documentation. The authors make no claims     *
9  * about the suitability of this software for any purpose. It is          *
10  * provided "as is" without express or implied warranty.                  *
11  **************************************************************************/
12 //-----------------------------------------------------------------
13 //                 AliAnalysisTaskNucleiv2SP class
14 //-----------------------------------------------------------------
15
16 class TTree;
17 class TParticle;
18 class TVector3;
19
20 #include "AliAnalysisManager.h"
21 #include <AliMCEventHandler.h>
22 #include <AliMCEvent.h>
23 #include <AliStack.h>
24
25 class AliESDVertex;
26 class AliAODVertex;
27 class AliESDv0;
28 class AliAODv0; 
29
30 #include <iostream>
31
32 #include "TList.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "TH3.h"
36 #include "TNtuple.h"
37 #include "TGraph.h"
38 #include "TF1.h"
39 #include "TCanvas.h"
40 #include "TMath.h"
41 #include "TChain.h"
42 #include "Riostream.h"
43 #include "AliLog.h"
44 #include "AliCascadeVertexer.h"
45 #include "AliESDEvent.h"
46 #include "AliESDtrack.h"
47 #include "AliExternalTrackParam.h"
48 #include "AliAODEvent.h"
49 #include "AliInputEventHandler.h"
50 #include "AliESDcascade.h"
51 #include "AliAODcascade.h"
52 #include "AliAnalysisTaskNucleiv2SP.h"
53 #include "AliESDtrackCuts.h"
54 #include "AliCentrality.h"
55
56
57 ClassImp(AliAnalysisTaskNucleiv2SP)
58
59 using std::cout;
60 using std::endl;
61     
62 //________________________________________________________________________
63 AliAnalysisTaskNucleiv2SP::AliAnalysisTaskNucleiv2SP() 
64 : AliAnalysisTaskSE(), 
65   fisPrimCut(kFALSE),
66   fptc(1),     
67   fListHist(0), 
68   fHistEventMultiplicity(0), 
69   fHistTrackMultiplicity(0),
70   fHistTrackMultiplicityCentral(0),    
71   fHistTrackMultiplicitySemiCentral(0),
72   fHistTrackMultiplicityMB(0),
73   fhBB(0),
74   fhBBDeu(0),
75   fhTOF(0),
76   fhMassTOF(0),
77   EPVzAvsCentrality(0), 
78   EPVzCvsCentrality(0), 
79   EPTPCvsCentrality(0), 
80   EPVzvsCentrality(0), 
81   EPTPCpvsCentrality(0), 
82   EPTPCnvsCentrality(0), 
83   hEvPlaneTPCvsEvPVz05(0),                      
84   hEvPlaneTPCvsEvPVz075(0), 
85   hEvPlaneTPCvsEvPVz1530(0),
86   hEvPlaneTPCvsEvPVz3050(0),                      
87   hEvPlaneTPCvsEvPVz2040(0),                      
88   hEvPlaneTPCvsEvPVz4060(0),       
89   hCos2DeltaTPCVzAvsCentrality(0),
90   hCos2DeltaTPCVzCvsCentrality(0),
91   hCos2DeltaVzAVzCvsCentrality(0),
92   hCos2DeltaVzMVzAvsCentrality(0),
93   hCos2DeltaVzMVzCvsCentrality(0),
94   hCos2DeltaVzATPCvsCentrality(0),
95   hCos2DeltaVzCTPCvsCentrality(0),
96   hCos2DeltaVzCVzAvsCentrality(0),
97   hCos2DeltaVzMTPCpvsCentrality(0),
98   hCos2DeltaVzMTPCnvsCentrality(0),
99   hCos2DeltaTPCpTPCnvsCentrality(0),
100   hQVzAQVzCvsCentrality(0),
101   hQxVzAvsCentrality(0),
102   hQyVzAvsCentrality(0),
103   hQxVzCvsCentrality(0),
104   hQyVzCvsCentrality(0),
105   hQxVzMvsCentrality(0),
106   hQyVzMvsCentrality(0),
107   ftree(0),           
108   tCentrality(0),     
109   tType(0),  
110   tHasTOF(0),    
111   tpT(0),  
112   tMassTOF(0),
113   tuqV0A(0),
114   tuqV0C(0),
115   tCharge(0),
116   tCosdeltaphiTPC(0),
117   tCosdeltaphiV0M(0),
118   tCosdeltaphiV0A(0),
119   tCosdeltaphiV0C(0),
120   timpactXY(0),
121   timpactZ(0),
122   tpull(0),
123   tphi(0),
124   fESDtrackCuts(0),
125   fESDtrackCutsEP(0),
126   fPIDResponse(0)
127 {
128   // Dummy Constructor 
129   fESDtrackCuts   = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
130   fESDtrackCutsEP = new AliESDtrackCuts("AliESDtrackCutsEP","AliESDtrackCutsEP");
131   Initialize();
132   cout<<"Dummy constructor"<<endl;
133 }
134
135 //________________________________________________________________________
136 AliAnalysisTaskNucleiv2SP::AliAnalysisTaskNucleiv2SP(const char *name) 
137 : AliAnalysisTaskSE(name), 
138   fisPrimCut(kFALSE),
139   fptc(1),     
140   fListHist(0), 
141   fHistEventMultiplicity(0), 
142   fHistTrackMultiplicity(0),
143   fHistTrackMultiplicityCentral(0),    
144   fHistTrackMultiplicitySemiCentral(0),
145   fHistTrackMultiplicityMB(0),
146   fhBB(0),
147   fhBBDeu(0),
148   fhTOF(0),
149   fhMassTOF(0),
150   EPVzAvsCentrality(0), 
151   EPVzCvsCentrality(0), 
152   EPTPCvsCentrality(0), 
153   EPVzvsCentrality(0), 
154   EPTPCpvsCentrality(0), 
155   EPTPCnvsCentrality(0), 
156   hEvPlaneTPCvsEvPVz05(0),                      
157   hEvPlaneTPCvsEvPVz075(0), 
158   hEvPlaneTPCvsEvPVz1530(0),
159   hEvPlaneTPCvsEvPVz3050(0),                      
160   hEvPlaneTPCvsEvPVz2040(0),                      
161   hEvPlaneTPCvsEvPVz4060(0),       
162   hCos2DeltaTPCVzAvsCentrality(0),
163   hCos2DeltaTPCVzCvsCentrality(0),
164   hCos2DeltaVzAVzCvsCentrality(0),
165   hCos2DeltaVzMVzAvsCentrality(0),
166   hCos2DeltaVzMVzCvsCentrality(0),
167   hCos2DeltaVzATPCvsCentrality(0),
168   hCos2DeltaVzCTPCvsCentrality(0),
169   hCos2DeltaVzCVzAvsCentrality(0),
170   hCos2DeltaVzMTPCpvsCentrality(0),
171   hCos2DeltaVzMTPCnvsCentrality(0),
172   hCos2DeltaTPCpTPCnvsCentrality(0),
173   hQVzAQVzCvsCentrality(0),
174   hQxVzAvsCentrality(0),
175   hQyVzAvsCentrality(0),
176   hQxVzCvsCentrality(0),
177   hQyVzCvsCentrality(0),
178   hQxVzMvsCentrality(0),
179   hQyVzMvsCentrality(0),
180   ftree(0),           
181   tCentrality(0),     
182   tType(0),  
183   tHasTOF(0),    
184   tpT(0),  
185   tMassTOF(0),
186   tuqV0A(0),
187   tuqV0C(0),
188   tCharge(0),
189   tCosdeltaphiTPC(0),
190   tCosdeltaphiV0M(0),
191   tCosdeltaphiV0A(0),
192   tCosdeltaphiV0C(0),
193   timpactXY(0),
194   timpactZ(0),
195   tpull(0),
196   tphi(0),
197   fESDtrackCuts(0),
198   fESDtrackCutsEP(0),
199   fPIDResponse(0)
200 {
201   // Define input and output slots here
202   // Input slot #0 works with a TChain
203   //DefineInput(0, TChain::Class());
204   // Output slot #0 writes into a TList container ()
205
206   //
207   // create track cuts
208   //
209   fESDtrackCuts   = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
210   fESDtrackCutsEP = new AliESDtrackCuts("AliESDtrackCutsEP","AliESDtrackCutsEP");
211   //
212   cout<<"Real constructor"<<endl;
213   Initialize();
214
215   DefineInput(0, TChain::Class());
216   DefineOutput(1, TList::Class());
217   DefineOutput(2, TTree::Class());
218   
219 }
220
221 void AliAnalysisTaskNucleiv2SP::Initialize()
222 {
223   //
224   // updating parameters in case of changes
225   //
226   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(fisPrimCut,kTRUE);
227   fESDtrackCuts->SetMaxDCAToVertexXY(3);
228   fESDtrackCuts->SetMaxDCAToVertexZ(2);
229   fESDtrackCuts->SetEtaRange(-0.8,0.8);
230   
231   fESDtrackCutsEP = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); 
232
233 }
234
235 //________________________________________________________________________
236 Float_t AliAnalysisTaskNucleiv2SP::GetEventPlaneForCandidate(AliESDtrack* track0, const TVector2* q,AliEventplane *pl){
237   
238   // remove autocorrelations 
239   
240   TArrayF* qx = 0x0;
241   TArrayF* qy = 0x0;
242   TVector2 qcopy; 
243   // if(!fEtaGap){
244   qx = pl->GetQContributionXArray();
245   qy = pl->GetQContributionYArray();
246   qcopy = *q;
247   
248   TVector2 q0;
249   if((track0->GetID()) < qx->fN){
250     q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));
251   }
252   
253   qcopy = qcopy - q0;
254   
255   return qcopy.Phi()/2.;
256   
257 }
258 //________________________________________________________________________
259 Float_t AliAnalysisTaskNucleiv2SP::GetPhi0Pi(Float_t phi){
260   // Sets the phi angle in the range 0-pi
261   Float_t result=phi;
262   while(result<0){
263     result=result+TMath::Pi();
264   }
265   while(result>TMath::Pi()){
266     result=result-TMath::Pi();
267   }
268    return result;
269 }
270
271 //==================DEFINITION OF OUTPUT OBJECTS==============================
272
273 void AliAnalysisTaskNucleiv2SP::UserCreateOutputObjects()
274 {
275   //-------------------------------------------------------
276   fListHist = new TList();
277   fListHist->SetOwner();  // IMPORTANT!
278   
279   if(! fHistEventMultiplicity ){
280
281     fHistEventMultiplicity   = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 12 , -0.5,11.5);
282
283     fHistEventMultiplicity->GetXaxis()->SetBinLabel(1,"All Events");
284     fHistEventMultiplicity->GetXaxis()->SetBinLabel(2,"Events w/PV");
285     fHistEventMultiplicity->GetXaxis()->SetBinLabel(3,"Events w/|Vz|<10cm");
286     fHistEventMultiplicity->GetXaxis()->SetBinLabel(4,"Central Events");
287     fHistEventMultiplicity->GetXaxis()->SetBinLabel(5,"Semi-Central Events");
288     fHistEventMultiplicity->GetXaxis()->SetBinLabel(6,"MB Events");
289     //from HF
290     fHistEventMultiplicity->GetXaxis()->SetBinLabel(7,"nEventsAnal");
291     fHistEventMultiplicity->GetXaxis()->SetBinLabel(8,"nEvSelected");
292     fHistEventMultiplicity->GetXaxis()->SetBinLabel(9,"nCandidatesSelected");
293     fHistEventMultiplicity->GetXaxis()->SetBinLabel(10,"out of pt bounds");
294     fHistEventMultiplicity->GetXaxis()->SetBinLabel(11,"mismatch lab");
295     fHistEventMultiplicity->GetXaxis()->SetBinLabel(12,"non valid TPC EP");
296     fListHist->Add(fHistEventMultiplicity);
297   }
298
299   if(! fHistTrackMultiplicity ){
300     fHistTrackMultiplicity  = new TH2F( "fHistTrackMultiplicity", "Nb of Tracks MB Events |Vz| < 10", 250,0, 25000,105,-0.5,104.5);
301     fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");
302     fHistTrackMultiplicity->GetYaxis()->SetTitle("Percentile");
303     fListHist->Add(fHistTrackMultiplicity);
304   } 
305
306   if(! fHistTrackMultiplicityCentral ){
307     fHistTrackMultiplicityCentral  = new TH2F( "fHistTrackMultiplicityCentral", "Nb of Tracks MB Events |Vz| < 10", 250,0, 25000,105,-0.5,104.5);
308     fHistTrackMultiplicityCentral->GetXaxis()->SetTitle("Number of tracks");
309     fHistTrackMultiplicityCentral->GetYaxis()->SetTitle("Percentile");
310     fListHist->Add(fHistTrackMultiplicityCentral);
311   } 
312   if(! fHistTrackMultiplicitySemiCentral ){
313     fHistTrackMultiplicitySemiCentral  = new TH2F( "fHistTrackMultiplicitySemiCentral", "Nb of Tracks MB Events |Vz| < 10", 250,0, 25000,105,-0.5,104.5);
314     fHistTrackMultiplicitySemiCentral->GetXaxis()->SetTitle("Number of tracks");
315     fHistTrackMultiplicitySemiCentral->GetYaxis()->SetTitle("Percentile");
316     fListHist->Add(fHistTrackMultiplicitySemiCentral);
317   } 
318   if(! fHistTrackMultiplicityMB ){
319     fHistTrackMultiplicityMB  = new TH2F( "fHistTrackMultiplicityMB", "Nb of Tracks MB Events |Vz| < 10", 250,0, 25000,105,-0.5,104.5);
320     fHistTrackMultiplicityMB->GetXaxis()->SetTitle("Number of tracks");
321     fHistTrackMultiplicityMB->GetYaxis()->SetTitle("Percentile");
322     fListHist->Add(fHistTrackMultiplicityMB);
323   } 
324  
325   if(! fhBB ){
326     fhBB = new TH2F( "fhBB" , "BetheBlochTPC" , 240,-6,6,250,0,1000);
327     fListHist->Add(fhBB);
328   }
329   
330   if(! fhBBDeu ){
331     fhBBDeu = new TH2F( "fhBBDeu" , "BetheBlochTPC - Deuteron" , 240,-6,6,250,0,1000);
332     fListHist->Add(fhBBDeu);
333   }
334  
335   if(! fhTOF ){
336     fhTOF = new TH2F( "fhTOF" , "Scatter Plot TOF" , 240,-6,6,500,0,1.2);
337     fListHist->Add(fhTOF);
338   }
339   if(! fhMassTOF){
340     fhMassTOF=new TH1F ("fhMassTOF","Particle Mass - TOF", 100,0 ,10);
341     fListHist->Add(fhMassTOF);
342   }
343   
344   EPVzAvsCentrality  = new TH2D("EPVzAvsCentrality" , "EPVzAvsCentrality" , 80, -2, 2,105,-0.5,105.5);
345   EPVzCvsCentrality  = new TH2D("EPVzCvsCentrality" , "EPVzCvsCentrality" , 80, -2, 2,105,-0.5,105.5);
346   EPTPCvsCentrality  = new TH2D("EPTPCvsCentrality" , "EPTPCvsCentrality" , 80, -2, 2,105,-0.5,105.5);
347   EPVzvsCentrality   = new TH2D("EPVzvsCentrality"  , "EPVzvsCentrality"  , 80, -2, 2,105,-0.5,105.5);
348   EPTPCpvsCentrality = new TH2D("EPTPCpvsCentrality", "EPTPCpvsCentrality", 80, -2, 2,105,-0.5,105.5);
349   EPTPCnvsCentrality = new TH2D("EPTPCnvsCentrality", "EPTPCnvsCentrality", 80, -2, 2,105,-0.5,105.5);
350
351   fListHist->Add(EPVzAvsCentrality);
352   fListHist->Add(EPVzCvsCentrality);
353   fListHist->Add(EPTPCvsCentrality);
354   fListHist->Add(EPVzvsCentrality);
355   fListHist->Add(EPTPCpvsCentrality);
356   fListHist->Add(EPTPCnvsCentrality);
357   
358   hEvPlaneTPCvsEvPVz05   = new TH2F("hEvPlaneTPCvsEvPVz05"  ,"hEvPlaneTPCvsEvPVz05"  ,100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());                      
359   hEvPlaneTPCvsEvPVz075  = new TH2F("hEvPlaneTPCvsEvPVz075" ,"hEvPlaneTPCvsEvPVz075" ,100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi()); 
360   hEvPlaneTPCvsEvPVz1530 = new TH2F("hEvPlaneTPCvsEvPVz1530","hEvPlaneTPCvsEvPVz1530",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
361   hEvPlaneTPCvsEvPVz3050 = new TH2F("hEvPlaneTPCvsEvPVz3050","hEvPlaneTPCvsEvPVz3050",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());                      
362   hEvPlaneTPCvsEvPVz2040 = new TH2F("hEvPlaneTPCvsEvPVz2040","hEvPlaneTPCvsEvPVz2040",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());                      
363   hEvPlaneTPCvsEvPVz4060 = new TH2F("hEvPlaneTPCvsEvPVz4060","hEvPlaneTPCvsEvPVz4060",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());   
364   
365   fListHist->Add(hEvPlaneTPCvsEvPVz05);                      
366   fListHist->Add(hEvPlaneTPCvsEvPVz075); 
367   fListHist->Add(hEvPlaneTPCvsEvPVz1530);
368   fListHist->Add(hEvPlaneTPCvsEvPVz3050);                      
369   fListHist->Add(hEvPlaneTPCvsEvPVz2040);                      
370   fListHist->Add(hEvPlaneTPCvsEvPVz4060);   
371
372   hCos2DeltaTPCVzAvsCentrality   = new TH2F("hCos2DeltaTPCVzAvsCentrality"  ,"hCos2DeltaTPCVzAvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
373   hCos2DeltaTPCVzCvsCentrality   = new TH2F("hCos2DeltaTPCVzCvsCentrality"  ,"hCos2DeltaTPCVzCvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
374   hCos2DeltaVzAVzCvsCentrality   = new TH2F("hCos2DeltaVzAVzCvsCentrality"  ,"hCos2DeltaVzAVzCvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
375   hCos2DeltaVzMVzAvsCentrality   = new TH2F("hCos2DeltaVzMVzAvsCentrality"  ,"hCos2DeltaVzMVzAvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
376   hCos2DeltaVzMVzCvsCentrality   = new TH2F("hCos2DeltaVzMVzCvsCentrality"  ,"hCos2DeltaVzMVzCvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
377   hCos2DeltaVzATPCvsCentrality   = new TH2F("hCos2DeltaVzATPCvsCentrality"  ,"hCos2DeltaVzATPCvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
378   hCos2DeltaVzCTPCvsCentrality   = new TH2F("hCos2DeltaVzCTPCvsCentrality"  ,"hCos2DeltaVzCTPCvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
379   hCos2DeltaVzCVzAvsCentrality   = new TH2F("hCos2DeltaVzCVzAvsCentrality"  ,"hCos2DeltaVzCVzAvsCentrality"  ,100,-1.1,1.1,105,-0.5,105.5);
380   hCos2DeltaVzMTPCpvsCentrality  = new TH2F("hCos2DeltaVzMTPCpvsCentrality" ,"hCos2DeltaVzMTPCpvsCentrality" ,100,-1.1,1.1,105,-0.5,105.5);
381   hCos2DeltaVzMTPCnvsCentrality  = new TH2F("hCos2DeltaVzMTPCnvsCentrality" ,"hCos2DeltaVzMTPCnvsCentrality" ,100,-1.1,1.1,105,-0.5,105.5);
382   hCos2DeltaTPCpTPCnvsCentrality = new TH2F("hCos2DeltaTPCpTPCnvsCentrality","hCos2DeltaTPCpTPCnvsCentrality",100,-1.1,1.1,105,-0.5,105.5);
383
384   fListHist->Add(hCos2DeltaTPCVzAvsCentrality);
385   fListHist->Add(hCos2DeltaTPCVzCvsCentrality);
386   fListHist->Add(hCos2DeltaVzAVzCvsCentrality);
387   fListHist->Add(hCos2DeltaVzMVzAvsCentrality);
388   fListHist->Add(hCos2DeltaVzMVzCvsCentrality);
389   fListHist->Add(hCos2DeltaVzATPCvsCentrality);
390   fListHist->Add(hCos2DeltaVzCTPCvsCentrality);
391   fListHist->Add(hCos2DeltaVzCVzAvsCentrality);
392   fListHist->Add(hCos2DeltaVzMTPCpvsCentrality);  
393   fListHist->Add(hCos2DeltaVzMTPCnvsCentrality); 
394   fListHist->Add(hCos2DeltaTPCpTPCnvsCentrality);
395
396  
397   hQVzAQVzCvsCentrality = new TH2F("hQVzAQVzCvsCentrality","hQVzAQVzCvsCentrality",1000,-5,5,105,-0.5,105.5);
398   fListHist->Add(hQVzAQVzCvsCentrality);
399  
400   hQxVzAvsCentrality = new TH2F("hQxVzAvsCentrality","hQxVzAvsCentrality",100,-5,5,105,-0.5,105.5);
401   hQyVzAvsCentrality = new TH2F("hQyVzAvsCentrality","hQyVzAvsCentrality",100,-5,5,105,-0.5,105.5);
402   hQxVzCvsCentrality = new TH2F("hQxVzCvsCentrality","hQxVzCvsCentrality",100,-5,5,105,-0.5,105.5);
403   hQyVzCvsCentrality = new TH2F("hQyVzCvsCentrality","hQyVzCvsCentrality",100,-5,5,105,-0.5,105.5);
404   hQxVzMvsCentrality = new TH2F("hQxVzCvsCentrality","hQxVzCvsCentrality",100,-5,5,105,-0.5,105.5);
405   hQyVzMvsCentrality = new TH2F("hQyVzCvsCentrality","hQyVzCvsCentrality",100,-5,5,105,-0.5,105.5);
406
407   fListHist->Add(hQxVzAvsCentrality);
408   fListHist->Add(hQyVzAvsCentrality);
409   fListHist->Add(hQxVzCvsCentrality);
410   fListHist->Add(hQyVzCvsCentrality);
411   fListHist->Add(hQxVzMvsCentrality);
412   fListHist->Add(hQyVzMvsCentrality);
413
414   if(!ftree){
415    
416     ftree = new TTree("ftree","ftree");
417  
418     ftree->Branch("tCentrality"      ,&tCentrality      ,"tCentrality/D"    );
419     ftree->Branch("tType"            ,&tType            ,"tType/D"          );
420     ftree->Branch("tHasTOF"          ,&tHasTOF          ,"tHasTOF/D"        );
421     ftree->Branch("tpT"              ,&tpT              ,"tpT/D"            );
422     ftree->Branch("tMassTOF"         ,&tMassTOF         ,"tMassTOF/D"       );
423     ftree->Branch("tuqV0A"           ,&tuqV0A           ,"tuqV0A/D"         );
424     ftree->Branch("tuqV0C"           ,&tuqV0C           ,"tuqV0C/D"         );
425     ftree->Branch("tCharge"          ,&tCharge          ,"tCharge/D"        );
426     ftree->Branch("tCosdeltaphiTPC"  ,&tCosdeltaphiTPC  ,"tCosdeltaphiTPC/D");
427     ftree->Branch("tCosdeltaphiV0M"  ,&tCosdeltaphiV0M  ,"tCosdeltaphiV0M/D");
428     ftree->Branch("tCosdeltaphiV0A"  ,&tCosdeltaphiV0A  ,"tCosdeltaphiV0A/D");
429     ftree->Branch("tCosdeltaphiV0C"  ,&tCosdeltaphiV0C  ,"tCosdeltaphiV0C/D");
430     ftree->Branch("timpactXY"        ,&timpactXY        ,"timpactXY/D"      );
431     ftree->Branch("timpactZ"         ,&timpactZ         ,"timpactZ/D"       );
432     ftree->Branch("tpull"            ,&tpull            ,"tpull/D"          );
433     ftree->Branch("tphi"             ,&tphi             ,"tphi/D"          );
434
435   }
436
437   PostData(1,  fListHist);
438   PostData(2,  ftree);
439 }// end UserCreateOutputObjects
440
441
442 //====================== USER EXEC ========================
443
444 void AliAnalysisTaskNucleiv2SP::UserExec(Option_t *) 
445 {
446   // Main loop
447   // Called for EACH event
448   //  cout<<"AliAnalysisTaskNucleiv2SP Starting UserExec"<<endl;
449
450   Info("AliAnalysisTaskNucleiv2SP","Starting UserExec");  
451   
452   AliVEvent *event = InputEvent();
453   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
454   
455  
456   AliESDEvent* lESDevent = dynamic_cast<AliESDEvent*>(event);
457   if (!lESDevent) {
458     AliError("Cannot get the ESD event");
459     return;
460   }  
461   
462   fHistEventMultiplicity->Fill(1);
463   fHistEventMultiplicity->Fill(7);
464   
465   //_____________________________________________________
466   //   Centrality  
467   
468   AliCentrality *centrality = lESDevent->GetCentrality();
469   Float_t percentile=centrality->GetCentralityPercentile("V0M");
470   
471   Int_t TrackNumber = lESDevent->GetNumberOfTracks();
472   fHistTrackMultiplicity->Fill(TrackNumber,percentile); //tracce per evento
473   
474   //______________________________________________________
475   // PID
476   
477   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
478   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
479   fPIDResponse=inputHandler->GetPIDResponse(); 
480   
481   //=================================================================
482   
483   Float_t  impactXY=-999., impactZ=-999.;
484   Double_t TPCSignal=0.;
485   
486   ULong_t  status=0;
487  
488   Double_t pmax  = 10;
489   Double_t ptmax = 7;
490   // Primary vertex cut
491   
492   const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();
493     
494   if(vtx->GetNContributors()<1) {
495     
496     // SPD vertex cut
497     vtx = lESDevent->GetPrimaryVertexSPD();
498     
499     if(vtx->GetNContributors()<1) {
500       Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
501       return; // NO GOOD VERTEX, SKIP EVENT 
502     }
503   }
504   
505   fHistEventMultiplicity->Fill(2); // analyzed events with PV
506   
507   if(TMath::Abs(vtx->GetZ())>10) return;
508   fHistEventMultiplicity->Fill(3);
509
510   Bool_t isSelectedCentral     = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
511   Bool_t isSelectedSemiCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
512   Bool_t isSelectedMB          = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
513     
514   fHistTrackMultiplicity->Fill(TrackNumber,percentile); 
515     
516   Int_t eventtype = -999;
517   
518   //  cout<<"ET 1: "<<eventtype<<endl;
519   
520   if(isSelectedCentral){
521     fHistEventMultiplicity->Fill(4);
522     fHistTrackMultiplicityCentral->Fill(TrackNumber,percentile); 
523     eventtype =1;
524   }
525   
526   if(isSelectedSemiCentral){
527     fHistEventMultiplicity->Fill(5);
528     fHistTrackMultiplicitySemiCentral->Fill(TrackNumber,percentile); 
529     eventtype =2;
530   }
531   
532   if(isSelectedMB){
533     if(percentile<0)return;
534     if(percentile>=80)return;
535     fHistEventMultiplicity->Fill(6);
536     fHistTrackMultiplicityMB->Fill(TrackNumber,percentile); 
537     eventtype =3;
538   }
539   
540   //    cout<<"ET 2: "<<eventtype<<endl;
541   
542   if(eventtype!=1 && eventtype!=2 && eventtype!=3 )return;
543   
544   AliEventplane *pl=lESDevent->GetEventplane();
545   
546   
547   if(!pl ){
548     AliError("AliAnalysisTaskSENucleiv2SP::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
549     fHistEventMultiplicity->Fill(12);
550   }
551     
552   //Event plane from FLOW
553   
554   Double_t qxEPa = 0, qyEPa = 0;
555   Double_t qxEPc = 0, qyEPc = 0;
556   Double_t qxEP =  0 , qyEP = 0;
557   
558   Double_t evPlAngV0A = pl->CalculateVZEROEventPlane(lESDevent, 8, 2, qxEPa, qyEPa);
559   Double_t evPlAngV0C = pl->CalculateVZEROEventPlane(lESDevent, 9, 2, qxEPc, qyEPc);
560   Double_t evPlAngV0  = pl->CalculateVZEROEventPlane(lESDevent,10, 2, qxEP,  qyEP);
561   
562   Double_t Qx2  = 0, Qy2  = 0;
563   Double_t Qx2p = 0, Qy2p = 0;
564   Double_t Qx2n = 0, Qy2n = 0;
565  
566   for (Int_t iT = 0; iT < TrackNumber; iT++){
567     
568     AliESDtrack* track = lESDevent->GetTrack(iT);
569     
570     if (!track)
571       continue;
572     
573     if ((TMath::Abs(track->Eta()) > 0.8) || (track->Pt() < 0.2) || (track->GetTPCNcls() < 70) || (track->Pt() >= 20.0))
574       continue;
575     if(!fESDtrackCutsEP->AcceptTrack(track))
576       continue;
577     if(track->Eta()>0 && track->Eta()<0.8){
578       
579       Qx2p += TMath::Cos(2*track->Phi());
580       Qy2p += TMath::Sin(2*track->Phi());
581     }
582     if(track->Eta()<0 && track->Eta()> -0.8){
583       
584       Qx2n += TMath::Cos(2*track->Phi());
585       Qy2n += TMath::Sin(2*track->Phi());
586     }
587
588     if(track->Eta()>0 && track->Eta()<0.8){ //half TPC
589       Qx2 += TMath::Cos(2*track->Phi());
590       Qy2 += TMath::Sin(2*track->Phi());
591     }
592   }
593   
594   Double_t evPlAngTPC  = TMath::ATan2(Qy2,  Qx2) /2.;
595   Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
596   Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
597
598   EPVzAvsCentrality  ->Fill(evPlAngV0A  , percentile); 
599   EPVzCvsCentrality  ->Fill(evPlAngV0C  , percentile); 
600   EPVzvsCentrality   ->Fill(evPlAngV0   , percentile); 
601   EPTPCvsCentrality  ->Fill(evPlAngTPC  , percentile); 
602   EPTPCpvsCentrality ->Fill(evPlAngTPCp , percentile); 
603   EPTPCnvsCentrality ->Fill(evPlAngTPCn , percentile); 
604
605   if(percentile>=0 && percentile<=5)
606     hEvPlaneTPCvsEvPVz05  ->Fill(evPlAngTPC,evPlAngV0); 
607   if(percentile>=0 && percentile<=7.5)
608     hEvPlaneTPCvsEvPVz075 ->Fill(evPlAngTPC,evPlAngV0); 
609   if(percentile>=15 && percentile<=30)
610     hEvPlaneTPCvsEvPVz1530->Fill(evPlAngTPC,evPlAngV0);
611   if(percentile>=30 && percentile<50)
612     hEvPlaneTPCvsEvPVz3050->Fill(evPlAngTPC,evPlAngV0);  
613   if(percentile>=20 && percentile<=40)                    
614     hEvPlaneTPCvsEvPVz2040->Fill(evPlAngTPC,evPlAngV0);   
615   if(percentile>=40 && percentile<=60)                   
616     hEvPlaneTPCvsEvPVz4060->Fill(evPlAngTPC,evPlAngV0);           
617
618   // For TPC, V0M, V0c and V0A resolution 
619
620   hCos2DeltaTPCVzAvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngTPC - evPlAngV0A)) , percentile);
621   hCos2DeltaTPCVzCvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngTPC - evPlAngV0C)) , percentile);
622   hCos2DeltaVzAVzCvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngV0A - evPlAngV0C)) , percentile);
623   hCos2DeltaVzMVzAvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngV0  - evPlAngV0A)) , percentile);
624   hCos2DeltaVzMVzCvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngV0  - evPlAngV0C)) , percentile);
625   hCos2DeltaVzATPCvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngV0A - evPlAngTPC)) , percentile);
626   hCos2DeltaVzCTPCvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngV0C - evPlAngTPC)) , percentile);
627   hCos2DeltaVzCVzAvsCentrality  ->Fill(TMath::Cos(2.*(evPlAngV0C - evPlAngV0A)) , percentile);
628   hCos2DeltaVzMTPCpvsCentrality ->Fill(TMath::Cos(2.*(evPlAngV0  - evPlAngTPCp)), percentile);
629   hCos2DeltaVzMTPCnvsCentrality ->Fill(TMath::Cos(2.*(evPlAngV0  - evPlAngTPCn)), percentile);
630   hCos2DeltaTPCpTPCnvsCentrality->Fill(TMath::Cos(2.*(evPlAngTPCp- evPlAngTPCn)), percentile);
631
632   //Scalar Product
633   
634   Double_t  QV0AQV0C = qxEPa * qxEPc + qyEPa*qyEPc;
635   hQVzAQVzCvsCentrality->Fill(QV0AQV0C,percentile);
636   
637   //NUA correction
638
639   hQxVzAvsCentrality->Fill(qxEPa,percentile);
640   hQyVzAvsCentrality->Fill(qyEPa,percentile);
641   hQxVzCvsCentrality->Fill(qxEPc,percentile);
642   hQyVzCvsCentrality->Fill(qyEPc,percentile);
643   hQxVzMvsCentrality->Fill(qxEP ,percentile);
644   hQyVzMvsCentrality->Fill(qyEP ,percentile);
645
646   //====================================================================================================================
647   
648   // To remove auto-correlation
649   TVector2 *q = 0x0;
650   q = pl->GetQVector();
651
652   Float_t ptcExp  = -999;
653   Double_t pullTPC = -999;
654   Float_t deltaphiTPC = -3;
655   Float_t deltaphiV0  = -3;
656   Float_t deltaphiV0A = -3;
657   Float_t deltaphiV0C = -3;
658
659   Float_t  uqV0A = -999;
660   Float_t  uqV0C = -999; 
661
662   for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks
663       
664     AliESDtrack *esdtrack=lESDevent->GetTrack(j);
665     if (!fESDtrackCuts->AcceptTrack(esdtrack)) continue;
666     
667     status  = (ULong_t)esdtrack->GetStatus();
668    
669     Bool_t hasTOFout  = status&AliESDtrack::kTOFout; 
670     Bool_t hasTOF     = kFALSE;
671     if (hasTOFout) hasTOF = kTRUE;
672     Float_t length = esdtrack->GetIntegratedLength(); 
673     if (length < 350.) hasTOF = kFALSE;
674     
675     TPCSignal=esdtrack->GetTPCsignal(); 
676       
677     if(TPCSignal<10)continue;
678     if(TPCSignal>1000)continue;
679     if(!esdtrack->GetInnerParam()) continue;
680     
681     Double_t ptot = esdtrack->GetInnerParam()->GetP(); // momentum for dEdx determination
682  
683     if(ptot<0.2)continue;
684     fhBB->Fill(ptot*esdtrack->GetSign(),TPCSignal);
685     esdtrack->GetImpactParameters(impactXY, impactZ);
686               
687     ptcExp = -999;
688     if(fptc==1)
689       ptcExp  = AliExternalTrackParam::BetheBlochAleph(ptot/(0.938*2),1.45802,27.4992,4.00313e-15,2.48485,8.31768);
690     if(fptc==2)
691       ptcExp  = AliExternalTrackParam::BetheBlochAleph(ptot/(0.938*3),1.45802,27.4992,4.00313e-15,2.48485,8.31768);
692     if(fptc==3)
693       ptcExp  = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),1.74962,27.4992,4.00313e-15,2.42485,8.31768);
694     
695     pullTPC  = (TPCSignal - ptcExp)/(0.07*ptcExp);
696
697     Double_t p    = esdtrack->P();
698     Double_t tof  = esdtrack->GetTOFsignal()-fPIDResponse->GetTOFResponse().GetStartTime(p);
699     Double_t tPhi = esdtrack->Phi();
700
701     Float_t  beta = 0;
702     Float_t  gamma = 0;
703     Float_t  mass  = -99;
704     Double_t pt  = esdtrack->Pt();
705
706     if(fptc==3)
707       pt = 2*pt;
708
709     if(TMath::Abs(ptot) < pmax && TMath::Abs(pullTPC) <= 3 && TMath::Abs(pt) < ptmax){
710
711       fhBBDeu->Fill(ptot*esdtrack->GetSign(),TPCSignal);
712       
713       //
714       // Process TOF information
715       //
716       if (hasTOF) {
717         beta = length / (2.99792457999999984e-02 * tof);
718         gamma = 1/TMath::Sqrt(1 - beta*beta);
719         mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
720         
721         fhTOF->Fill(ptot*esdtrack->GetSign(),beta);
722         if(fptc==1){
723           if(TMath::Abs(mass) > 2.7)continue;
724           if(TMath::Abs(mass) < 1. )continue;
725         }
726         if(fptc==2){
727           if(TMath::Abs(mass) > 5.0)continue;
728           if(TMath::Abs(mass) < 1.8 )continue;
729         }
730         if(fptc==3){
731           if(TMath::Abs(mass) > 5.0)continue;
732           if(TMath::Abs(mass) < 1.8 )continue;
733         }
734         fhMassTOF->Fill(mass);
735       }
736    
737       // Event Plane
738       // Remove AutoCorrelation
739       
740       evPlAngTPC = GetEventPlaneForCandidate(esdtrack,q,pl);
741       
742       deltaphiTPC=TMath::Cos(2*GetPhi0Pi(tPhi-evPlAngTPC));
743       deltaphiV0 =TMath::Cos(2*GetPhi0Pi(tPhi-evPlAngV0 ));
744       deltaphiV0A=TMath::Cos(2*GetPhi0Pi(tPhi-evPlAngV0A));
745       deltaphiV0C=TMath::Cos(2*GetPhi0Pi(tPhi-evPlAngV0C));
746       
747       // Scalar Product
748       
749       uqV0A = TMath::Cos(2*tPhi)*qxEPa+TMath::Sin(2*tPhi)*qyEPa;
750       uqV0C = TMath::Cos(2*tPhi)*qxEPc+TMath::Sin(2*tPhi)*qyEPc;
751        
752       tCentrality      = percentile;
753       tType            = eventtype;
754       tHasTOF          = hasTOF;
755       tpT              = pt;
756       tMassTOF         = mass;
757       tuqV0A           = uqV0A;
758       tuqV0C           = uqV0C;
759       tCharge          = esdtrack->GetSign();
760       tCosdeltaphiTPC  = deltaphiTPC;
761       tCosdeltaphiV0M  = deltaphiV0;
762       tCosdeltaphiV0A  = deltaphiV0A;
763       tCosdeltaphiV0C  = deltaphiV0C;
764       timpactXY        = impactXY;
765       timpactZ         = impactZ;
766       tpull            = pullTPC;
767       tphi             = tPhi;
768
769       ftree->Fill();
770     } 
771   }  //track
772   
773   PostData(1, fListHist);
774   PostData(2, ftree);
775 } //end userexec
776
777
778 //________________________________________________________________________
779
780 void AliAnalysisTaskNucleiv2SP::Terminate(Option_t *) 
781 {
782   // Draw result to the screen
783   // Called once at the end of the query
784 }
785