Transition PWG4/JetTaske -> PWGJE
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetServices.cxx
1
2 // **************************************
3 // Task used for the correction of determiantion of reconstructed jet spectra
4 // Compares input (gen) and output (rec) jets   
5 // *******************************************
6
7
8 /**************************************************************************
9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10  *                                                                        *
11  * Author: The ALICE Off-line Project.                                    *
12  * Contributors are mentioned in the code where appropriate.              *
13  *                                                                        *
14  * Permission to use, copy, modify and distribute this software and its   *
15  * documentation strictly for non-commercial purposes is hereby granted   *
16  * without fee, provided that the above copyright notice appears in all   *
17  * copies and that both the copyright notice and this permission notice   *
18  * appear in the supporting documentation. The authors make no claims     *
19  * about the suitability of this software for any purpose. It is          *
20  * provided "as is" without express or implied warranty.                  *
21  **************************************************************************/
22
23  
24 #include <TROOT.h>
25 #include <TRandom.h>
26 #include <TString.h>
27 #include <TSystem.h>
28 #include <TInterpreter.h>
29 #include <TChain.h>
30 #include <TFile.h>
31 #include <TKey.h>
32 #include <TH1F.h>
33 #include <TH2F.h>
34 #include <TH3F.h>
35 #include <TProfile.h>
36 #include <TList.h>
37 #include <TLorentzVector.h>
38 #include <TRandom3.h>
39 #include <TClonesArray.h>
40 #include "TDatabasePDG.h"
41
42 #include "AliAnalysisTaskJetServices.h"
43 #include "AliCentrality.h"
44 #include "AliAnalysisDataContainer.h"
45 #include "AliAnalysisDataSlot.h"
46 #include "AliAnalysisManager.h"
47 #include "AliJetFinder.h"
48 #include "AliJetHeader.h"
49 #include "AliJetReader.h"
50 #include "AliJetReaderHeader.h"
51 #include "AliUA1JetHeaderV1.h"
52 #include "AliESDEvent.h"
53 #include "AliAODEvent.h"
54 #include "AliAODHandler.h"
55 #include "AliAODTrack.h"
56 #include "AliAODVZERO.h"
57 #include "AliAODJet.h"
58 #include "AliAODMCParticle.h"
59 #include "AliMCEventHandler.h"
60 #include "AliMCEvent.h"
61 #include "AliStack.h"
62 #include "AliGenPythiaEventHeader.h"
63 #include "AliJetKineReaderHeader.h"
64 #include "AliGenCocktailEventHeader.h"
65 #include "AliInputEventHandler.h"
66 #include "AliPhysicsSelection.h"
67 #include "AliTriggerAnalysis.h"
68
69 #include "AliAnalysisHelperJetTasks.h"
70
71 ClassImp(AliAnalysisTaskJetServices)
72
73 AliAODVZERO*  AliAnalysisTaskJetServices::fgAODVZERO = NULL;
74 AliAODHeader*  AliAnalysisTaskJetServices::fgAODHeader = NULL;
75 TClonesArray*   AliAnalysisTaskJetServices::fgAODVertices = NULL;
76
77 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): 
78   AliAnalysisTaskSE(),
79   fUseAODInput(kFALSE),
80   fUsePhysicsSelection(kFALSE),
81   fMC(kFALSE),
82   fFilterAODCollisions(kFALSE),
83   fPhysicsSelectionFlag(AliVEvent::kMB),
84   fSelectionInfoESD(0),
85   fEventCutInfoESD(0),
86   fFilterMask(0),
87   fRPMethod(0),
88   fCollisionType(kPbPb),
89   fAvgTrials(1),
90   fVtxXMean(0),
91   fVtxYMean(0),
92   fVtxZMean(0),
93   fVtxRCut(1.),
94   fVtxZCut(8.),
95   fPtMinCosmic(5.),
96   fRIsolMinCosmic(3.),
97   fMaxCosmicAngle(0.01),
98   fCentrality(101),
99   fTrackRecEtaWindow(0.9),
100   fMinTrackPt(0.15),
101   fRPAngle(0),
102   fPsiVZEROA(0),
103   fPsiVZEROC(0),
104   fRandomizer(0),
105   fNonStdFile(""),
106   fh1Xsec(0x0),
107   fh1Trials(0x0),
108   fh1PtHard(0x0),
109   fh1PtHardTrials(0x0),
110   fh1SelectionInfoESD(0x0),
111   fh1EventCutInfoESD(0),
112   fh1CentralityESD(0),
113   fh1Centrality(0),
114   fh1RP(0),
115   fh2TriggerCount(0x0),
116   fh2ESDTriggerCount(0x0),
117   fh2TriggerVtx(0x0),
118   fh2ESDTriggerVtx(0x0),
119   fh2ESDTriggerRun(0x0),
120   fh2VtxXY(0x0),
121   fh1NCosmicsPerEvent(0x0),
122   fp1RPXA(0x0),
123   fp1RPYA(0x0),
124   fp1RPXC(0x0),
125   fp1RPYC(0x0),
126   fp1CalibRPXA(0x0),
127   fp1CalibRPYA(0x0),
128   fp1CalibRPXC(0x0),
129   fp1CalibRPYC(0x0),
130   fh2RPAC(0x0),
131   fh2RPAT(0x0),
132   fh2RPCT(0x0),
133   fh2XYA(0x0),
134   fh2XYC(0x0),
135   fh2RPCentrality(0x0),
136   fh2RPACentrality(0x0),
137   fh2RPCCentrality(0x0),
138   fTriggerAnalysis(0x0),
139   fHistList(0x0)  
140 {
141   fRunRange[0] = fRunRange[1] = 0; 
142 }
143
144 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
145   AliAnalysisTaskSE(name),
146   fUseAODInput(kFALSE),
147   fUsePhysicsSelection(kFALSE),
148   fMC(kFALSE),
149   fFilterAODCollisions(kFALSE),
150   fPhysicsSelectionFlag(AliVEvent::kMB),
151   fSelectionInfoESD(0),
152   fEventCutInfoESD(0),
153   fFilterMask(0),
154   fRPMethod(0),
155   fCollisionType(kPbPb),
156   fAvgTrials(1),
157   fVtxXMean(0),
158   fVtxYMean(0),
159   fVtxZMean(0),
160   fVtxRCut(1.),
161   fVtxZCut(8.),
162   fPtMinCosmic(5.),
163   fRIsolMinCosmic(3.),
164   fMaxCosmicAngle(0.01),
165   fCentrality(101),
166   fTrackRecEtaWindow(0.9),
167   fMinTrackPt(0.15),
168   fRPAngle(0),
169   fPsiVZEROA(0),
170   fPsiVZEROC(0),
171   fRandomizer(0),
172   fNonStdFile(""),
173   fh1Xsec(0x0),
174   fh1Trials(0x0),
175   fh1PtHard(0x0),
176   fh1PtHardTrials(0x0),
177   fh1SelectionInfoESD(0x0),
178   fh1EventCutInfoESD(0),
179   fh1CentralityESD(0),
180   fh1Centrality(0),
181   fh1RP(0),
182   fh2TriggerCount(0x0),
183   fh2ESDTriggerCount(0x0),
184   fh2TriggerVtx(0x0),
185   fh2ESDTriggerVtx(0x0),
186   fh2ESDTriggerRun(0x0),
187   fh2VtxXY(0x0),
188   fh1NCosmicsPerEvent(0x0),
189   fp1RPXA(0x0),
190   fp1RPYA(0x0),
191   fp1RPXC(0x0),
192   fp1RPYC(0x0),
193   fp1CalibRPXA(0x0),
194   fp1CalibRPYA(0x0),
195   fp1CalibRPXC(0x0),
196   fp1CalibRPYC(0x0),
197   fh2RPAC(0x0),
198   fh2RPAT(0x0),
199   fh2RPCT(0x0),
200   fh2XYA(0x0),
201   fh2XYC(0x0),
202   fh2RPCentrality(0x0),
203   fh2RPACentrality(0x0),
204   fh2RPCCentrality(0x0),
205   fTriggerAnalysis(0x0),
206   fHistList(0x0)  
207 {
208   fRunRange[0] = fRunRange[1] = 0; 
209   DefineOutput(1,TList::Class());
210 }
211
212
213
214 Bool_t AliAnalysisTaskJetServices::Notify()
215 {
216   //
217   // Implemented Notify() to read the cross sections
218   // and number of trials from pyxsec.root
219   // 
220
221   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
222   Float_t xsection = 0;
223   Float_t ftrials  = 1;
224
225   fAvgTrials = 1;
226   if(tree){
227     TFile *curfile = tree->GetCurrentFile();
228     if (!curfile) {
229       Error("Notify","No current file");
230       return kFALSE;
231     }
232     if(!fh1Xsec||!fh1Trials){
233       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
234       return kFALSE;
235     }
236     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
237     fh1Xsec->Fill("<#sigma>",xsection);
238     // construct a poor man average trials 
239     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
240     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
241   }  
242   return kTRUE;
243 }
244
245 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
246 {
247
248   //
249   // Create the output container
250   //
251
252   fRandomizer = new TRandom3(0);
253   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
254
255   OpenFile(1);
256   if(!fHistList)fHistList = new TList();
257   fHistList->SetOwner();
258   PostData(1, fHistList); // post data in any case once
259
260   Bool_t oldStatus = TH1::AddDirectoryStatus();
261   TH1::AddDirectory(kFALSE);
262   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
263   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
264   fHistList->Add(fh1Xsec);
265
266   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
267   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
268   fHistList->Add(fh1Trials);
269
270   const Int_t nBinPt = 100;
271   Double_t binLimitsPt[nBinPt+1];
272   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
273     if(iPt == 0){
274       binLimitsPt[iPt] = 0.0;
275     }
276     else {// 1.0
277       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
278     }
279   }
280   
281
282   fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
283   fHistList->Add(fh1CentralityESD);
284   
285   fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
286   fHistList->Add(fh1Centrality);
287
288   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
289   fHistList->Add(fh1RP);
290
291   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5); 
292   fHistList->Add(fh2TriggerCount);
293
294   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5); 
295   fHistList->Add(fh2ESDTriggerCount);
296   const Int_t nBins = 6*kConstraints;
297   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
298   fHistList->Add(fh2TriggerVtx);
299
300   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
301   fHistList->Add(fh2ESDTriggerVtx);
302   
303
304   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
305   fHistList->Add(fh1PtHard);
306   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
307   fHistList->Add(fh1PtHardTrials);
308   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
309                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
310   fHistList->Add(fh1SelectionInfoESD);
311
312   fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
313                                 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
314   fHistList->Add(fh1EventCutInfoESD);
315
316   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
317   // 3 triggers BB BE/EB EE
318
319   fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Eventclass vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
320   fHistList->Add(fh2ESDTriggerRun);
321
322   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
323   fHistList->Add(fh2VtxXY);
324   // =========== Switch on Sumw2 for all histos ===========
325   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
326     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
327     if (h1){
328       h1->Sumw2();
329       continue;
330     }
331     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
332     if(hn)hn->Sumw2();
333   }
334
335   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
336   fHistList->Add(fh1NCosmicsPerEvent),
337
338
339   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle from tracks" , 20, 0.,100., 180, 0, TMath::Pi());
340   fHistList->Add(fh2RPCentrality);
341
342
343   fh2RPACentrality = new TH2F("fh2RPACentrality" ,"Reaction Plane Angle from vzero A" , 20, 0.,100., 180, 0, TMath::Pi());
344   fHistList->Add(fh2RPACentrality);
345
346   fh2RPCCentrality = new TH2F("fh2RPCCentrality" ,"Reaction Plane Angle from vzero C" , 20, 0.,100., 180, 0, TMath::Pi());
347   fHistList->Add(fh2RPCCentrality);
348
349   fh2RPAC = new TH2F("fh2RPAC" ,"Reaction Plane Angle vzero a vs c" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
350   fHistList->Add( fh2RPAC);
351
352   fh2RPAT = new TH2F("fh2RPAT" ,"Reaction Plane Angle vzero a vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
353   fHistList->Add( fh2RPAT);
354
355   fh2RPCT = new TH2F("fh2RPCT" ,"Reaction Plane Angle vzero c vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
356   fHistList->Add( fh2RPCT);
357
358   fh2XYA = new TH2F("fh2XYA" ,"XY vzeroa ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
359   fHistList->Add(fh2XYA);
360   fh2XYC = new TH2F("fh2XYC" ,"XY vzeroc ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
361   fHistList->Add(fh2XYC);
362
363   // profiles for mean 
364   fp1RPXA = new TProfile("fp1RPXA","mean vzeroa x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
365   fHistList->Add(fp1RPXA);
366
367   fp1RPYA = new TProfile("fp1RPYA","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
368   fHistList->Add(fp1RPYA);
369
370
371   fp1RPXC = new TProfile("fp1RPXC","mean vzeroc x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
372   fHistList->Add(fp1RPXC);
373
374   fp1RPYC = new TProfile("fp1RPYC","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
375   fHistList->Add(fp1RPYC);
376
377
378   TH1::AddDirectory(oldStatus);
379
380   // Add an AOD branch for replication
381   if(fNonStdFile.Length()){
382      if (fDebug > 1) AliInfo("Replicating header");
383      fgAODHeader = new AliAODHeader;
384      AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
385      if (fDebug > 1) AliInfo("Replicating vzeros");
386      fgAODVZERO = new AliAODVZERO;
387      AddAODBranch("AliAODVZERO",&fgAODVZERO,fNonStdFile.Data());
388      if (fDebug > 1) AliInfo("Replicating primary vertices");
389      fgAODVertices = new TClonesArray("AliAODVertex",3);
390      fgAODVertices->SetName("vertices");
391      AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
392   }
393 }
394
395 void AliAnalysisTaskJetServices::Init()
396 {
397   //
398   // Initialization
399   //
400   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
401
402 }
403
404 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
405 {
406
407   //
408   // Execute analysis for current event
409   //
410  
411   AliAODEvent *aod = 0;
412   AliESDEvent *esd = 0;
413
414
415   
416   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
417   AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
418   AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
419   fSelectionInfoESD = 0; // reset
420   fEventCutInfoESD = 0; // reset
421   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
422
423
424   static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
425    
426
427
428   if(fUseAODInput){    
429     aod = dynamic_cast<AliAODEvent*>(InputEvent());
430     if(!aod){
431       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
432       return;
433     }
434     // fethc the header
435   }
436   else{
437     //  assume that the AOD is in the general output...
438     aod  = AODEvent();
439     if(!aod){
440       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
441       return;
442     }
443     esd = dynamic_cast<AliESDEvent*>(InputEvent());
444   }
445   if(aod&&fDebug>2){
446     aod->Print();
447     Printf("Vertices %d",aod->GetNumberOfVertices());
448     Printf("tracks %d",aod->GetNumberOfTracks());
449     Printf("jets %d",aod->GetNJets());
450   }
451   fSelectionInfoESD |= kNoEventCut;
452   fEventCutInfoESD |= kNoEventCut;
453
454   Bool_t esdVtxValid = false;
455   Bool_t esdVtxIn = false;
456   Bool_t aodVtxValid = false;
457   Bool_t aodVtxIn = false;
458
459   if(esd){
460     // trigger analyisis
461     if(!fTriggerAnalysis){
462       fTriggerAnalysis = new AliTriggerAnalysis;
463       fTriggerAnalysis->SetAnalyzeMC(fMC);
464       fTriggerAnalysis->SetSPDGFOThreshhold(1);
465     }
466     //    fTriggerAnalysis->FillTriggerClasses(esd);
467     Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
468     Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
469     Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG);
470     Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG);
471     Bool_t spdFO      = fTriggerAnalysis->SPDFiredChips(esd, 0);;
472     if(v0A)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0A;
473     if(v0C)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0C;
474     if(!(v0ABG||v0CBG))fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kNoV0BG;
475     if(spdFO)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kSPDFO;
476   }
477  
478   // Apply additional constraints
479   Bool_t esdEventSelected = IsEventSelected(esd);
480   Bool_t esdEventPileUp = IsEventPileUp(esd);
481   Bool_t esdEventCosmic = IsEventCosmic(esd);
482
483   Bool_t aodEventSelected = IsEventSelected(aod);
484
485   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
486
487
488   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
489   fh1EventCutInfoESD->Fill(fEventCutInfoESD);
490
491   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
492   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
493   if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
494   if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
495
496
497   // here we have all selection information, fill histogram
498   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
499     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
500   }
501   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); 
502
503   if(esd&&aod&&fDebug){
504     if(esdEventSelected&&!aodEventSelected){
505       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
506       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
507       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
508       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
509       vtxESD->Print();
510       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
511       vtxAOD->Print();
512     }
513   }
514
515   // loop over all possible triggers for esd 
516
517   Float_t cent = 0;
518   if(fCollisionType==kPbPb){
519     if(aod)cent = aod->GetHeader()->GetCentrality();
520     if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
521     if(cent<0)cent = 101;
522   }
523   fCentrality = cent;
524   fRPAngle = 0;
525   fPsiVZEROA = 0;
526   fPsiVZEROC = 0;
527
528
529   if(esd){
530     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
531     esdVtxValid = IsVertexValid(vtxESD);
532     esdVtxIn = IsVertexIn(vtxESD);
533     if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
534       if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
535       if(cent<=80&&esdVtxIn){
536         aodH->SetFillAOD(kTRUE);
537         aodH->SetFillExtension(kTRUE);
538       }
539     }
540
541
542     Float_t zvtx = vtxESD->GetZ();
543     Int_t  iCl = GetEventClass(esd);
544     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
545     Bool_t cand = physicsSelection;
546
547     if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
548     fh2ESDTriggerCount->Fill(0.,kAllTriggered); 
549     fh2ESDTriggerCount->Fill(iCl,kAllTriggered); 
550     if(cand){
551       fh2ESDTriggerCount->Fill(0.,kSelectedALICE); 
552       fh2ESDTriggerCount->Fill(iCl,kSelectedALICE); 
553       fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
554     }
555     //    if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
556     if(esdVtxValid){
557       fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
558       fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
559       fh2ESDTriggerVtx->Fill(iCl,zvtx);
560       if(esdVtxIn){
561         fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
562         fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
563         fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
564       }
565       if(cand){
566         fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
567         fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
568         fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
569       }
570     }
571
572     if(cand&&esdVtxIn&&iCl<5){
573       fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
574       fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
575       fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
576       fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
577       fh2ESDTriggerCount->Fill(iCl,kSelected);
578       fh2ESDTriggerCount->Fill(0.,kSelected);
579       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
580       if(esd->GetCentrality()){
581         Float_t tmpCent = 100;
582         tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
583         if(tmpCent<0)tmpCent = 101;
584         fh1CentralityESD->Fill(tmpCent);
585       }
586     }
587   }
588
589
590
591   if(aod){
592     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
593     aodVtxValid = IsVertexValid(vtxAOD);
594     aodVtxIn = IsVertexIn(vtxAOD);
595     Float_t zvtx = vtxAOD->GetZ();
596     Int_t  iCl = GetEventClass(aod);
597     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
598     Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
599     if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
600     fh2TriggerCount->Fill(0.,kAllTriggered); 
601     fh2TriggerCount->Fill(iCl,kAllTriggered); 
602     if(cand){
603       fh2TriggerCount->Fill(0.,kSelectedALICE); 
604       fh2TriggerCount->Fill(iCl,kSelectedALICE); 
605       fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
606     }
607     if(aodVtxValid){
608       fh2TriggerCount->Fill(0.,kTriggeredVertex);
609       fh2TriggerCount->Fill(iCl,kTriggeredVertex);
610       fh2TriggerVtx->Fill(iCl,zvtx);
611       if(aodVtxIn){
612         fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
613         fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
614         fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
615       }
616       if(cand){
617         fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
618         fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
619         fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
620       }
621     }
622
623     if(cand&&aodVtxIn&&iCl<5){
624       fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
625       fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
626       fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
627       fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
628       fh2TriggerCount->Fill(iCl,kSelected);
629       fh2TriggerCount->Fill(0.,kSelected);
630       fh1Centrality->Fill(cent);
631       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
632       if(aodH&&cand&&fFilterAODCollisions&&!esd){
633         if(fCentrality<=80&&aodVtxIn){
634           aodH->SetFillAOD(kTRUE);
635           aodH->SetFillExtension(kTRUE);
636         }
637       }
638
639       TList recTracks;
640       GetListOfTracks(&recTracks);
641       CalculateReactionPlaneAngleVZERO(aod);
642       fRPAngle = aod->GetHeader()->GetEventplane();
643       fh1RP->Fill(fRPAngle);
644       fh2RPCentrality->Fill(fCentrality,fRPAngle);
645       fh2RPACentrality->Fill(fCentrality,fPsiVZEROA);
646       fh2RPCCentrality->Fill(fCentrality,fPsiVZEROC);
647       fh2RPAC->Fill(fPsiVZEROA,fPsiVZEROC);
648       fh2RPAT->Fill(fPsiVZEROA,fRPAngle);
649       fh2RPCT->Fill(fPsiVZEROC,fRPAngle);
650       if(fRPMethod==kRPTracks)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
651       else if(fRPMethod==kRPVZEROA)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
652       else if(fRPMethod==kRPVZEROC)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
653
654       if(fUseAODInput&&fCentrality<=80){
655         if(fFilterAODCollisions&&aod){
656           aodH->SetFillAOD(kTRUE);
657           aodH->SetFillExtension(kTRUE);
658         }
659       }
660     }    
661   }
662
663   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
664
665   
666   Double_t ptHard = 0; 
667   Double_t nTrials = 1; // Trials for MC trigger 
668
669   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
670   AliMCEvent* mcEvent = MCEvent();
671   //    AliStack *pStack = 0; 
672   if(mcEvent){
673     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
674     if(pythiaGenHeader){
675       nTrials = pythiaGenHeader->Trials();
676       ptHard  = pythiaGenHeader->GetPtHard();
677       int iProcessType = pythiaGenHeader->ProcessType();
678       // 11 f+f -> f+f
679       // 12 f+barf -> f+barf
680       // 13 f+barf -> g+g
681       // 28 f+g -> f+g
682       // 53 g+g -> f+barf
683       // 68 g+g -> g+g
684       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
685       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
686       fh1PtHard->Fill(ptHard);
687       fh1PtHardTrials->Fill(ptHard,nTrials);
688
689     }// if pythia gen header
690   }
691
692   // trigger selection
693   
694   // replication of 
695   if(fNonStdFile.Length()&&aod){
696     if (fgAODHeader){
697       *fgAODHeader =  *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
698       Double_t q[kRPMethods] = {fRPAngle,fPsiVZEROA,fPsiVZEROC};
699       fgAODHeader->SetQTheta(q,kRPMethods);
700     }
701     if (fgAODVZERO){
702       *fgAODVZERO =  *(dynamic_cast<AliAODVZERO*>(aod->GetVZEROData()));
703     }
704     if(fgAODVertices){
705       fgAODVertices->Delete();
706       TClonesArray &vertices = *fgAODVertices;
707       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
708       // we only use some basic information, 
709       
710
711       Double_t pos[3];
712       Double_t covVtx[6];
713
714       vtxAOD->GetXYZ(pos); // position                                                                
715       vtxAOD->GetCovMatrix(covVtx); //covariance matrix                                                      
716       Int_t jVertices = 0;
717       AliAODVertex * vtx = new(vertices[jVertices++])
718         AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
719       vtx->SetName(vtxAOD->GetName());
720       vtx->SetTitle(vtxAOD->GetTitle());
721       TString vtitle = vtxAOD->GetTitle();
722       vtx->SetNContributors(vtxAOD->GetNContributors());
723
724       // Add SPD "main" vertex                                                                    
725       const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
726       vtxS->GetXYZ(pos); // position
727       vtxS->GetCovMatrix(covVtx); //covariance matrix                                                      
728       AliAODVertex * mVSPD = new(vertices[jVertices++])
729         AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
730       mVSPD->SetName(vtxS->GetName());
731       mVSPD->SetTitle(vtxS->GetTitle());
732       mVSPD->SetNContributors(vtxS->GetNContributors());
733       
734       // Add tpc only vertex
735       if(esd){
736         const AliESDVertex *vtxT =  esd->GetPrimaryVertexTPC();
737         vtxT->GetXYZ(pos); // position                                                    
738         vtxT->GetCovMatrix(covVtx); //covariance matrix                                               
739         AliAODVertex * mVTPC = new(vertices[jVertices++])
740           AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
741         mVTPC->SetName(vtxT->GetName());
742         mVTPC->SetTitle(vtxT->GetTitle());
743         mVTPC->SetNContributors(vtxT->GetNContributors());
744       }
745     }
746   }
747   
748   PostData(1, fHistList);
749 }
750
751 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
752   if(!esd)return kFALSE;
753   const AliESDVertex *vtx = esd->GetPrimaryVertex();
754   return IsVertexIn(vtx); // vertex in calls vertex valid
755 }
756
757 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
758   if(fgAODVertices){
759     fgAODVertices->Delete();
760     delete fgAODVertices;
761   }
762   delete fRandomizer;
763   if(fgAODHeader)delete fgAODHeader;
764   if(fgAODVZERO)delete fgAODVZERO;
765   delete fp1CalibRPXA;
766   delete fp1CalibRPYA;
767   delete fp1CalibRPXC;
768   delete fp1CalibRPYC;
769
770 }
771
772
773 void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfile *xc, TProfile *yc){
774
775   if(xa){
776     if(fp1CalibRPXA)delete fp1CalibRPXA;
777     fp1CalibRPXA =  (TProfile*)xa->Clone(Form("%sCalib",xa->GetName()));
778   }
779   else{
780     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
781   }
782   if(ya){
783     if(fp1CalibRPYA)delete fp1CalibRPYA;
784     fp1CalibRPYA =  (TProfile*)ya->Clone(Form("%sCalib",ya->GetName()));
785   }
786   else{
787     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
788   }
789   if(xc){
790     if(fp1CalibRPXC)delete fp1CalibRPXC;
791     fp1CalibRPXC =  (TProfile*)xc->Clone(Form("%sCalib",xc->GetName()));
792   }
793   else{
794     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
795   }
796   if(ya){
797     if(fp1CalibRPYC)delete fp1CalibRPYC;
798     fp1CalibRPYC =  (TProfile*)yc->Clone(Form("%sCalib",yc->GetName()));
799   }
800   else{
801     Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
802   }
803 }
804
805 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
806   if(!aod)return kFALSE;
807   const AliAODVertex *vtx = aod->GetPrimaryVertex();
808   return IsVertexIn(vtx); // VertexIn calls VertexValid
809 }
810
811 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
812
813   // We can check the type of the vertex by its name and title
814   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
815   // vtx       name                  title
816   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
817   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
818   // TPC       TPCVertex             VertexerTracksNoConstraint
819   // TPC       TPCVertex             VertexerTracksWithConstraint
820   // SPD       SPDVertex             vertexer: 3D
821   // SPD       SPDVertex             vertexer: Z
822   
823   if(!vtx){
824     Printf("%s:%d No ESD vertex found",(char*)__FILE__,__LINE__);
825     return kFALSE;
826   }
827   Int_t nCont = vtx->GetNContributors();
828   if(nCont>=1){
829     fEventCutInfoESD |= kContributorsCut1;    
830     if(nCont>=2){
831       fEventCutInfoESD |= kContributorsCut2;    
832       if(nCont>=3){
833         fEventCutInfoESD |= kContributorsCut3;    
834       }
835     }
836   }
837   
838   if(nCont<3)return kFALSE;
839
840   // do not want tpc only primary vertex
841   TString vtxName(vtx->GetName());
842   if(vtxName.Contains("TPCVertex")){
843     fEventCutInfoESD |= kVertexTPC;
844     return kFALSE;
845   }
846   if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
847   if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
848
849
850   TString vtxTitle(vtx->GetTitle());
851   if(vtxTitle.Contains("vertexer: Z")){
852     if(vtx->GetDispersion()>0.02)return kFALSE;   
853   }
854   fEventCutInfoESD |= kSPDDispersionCut;
855   return kTRUE;
856 }
857
858
859 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
860
861   // We can check the type of the vertex by its name and title
862   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
863   // vtx       name                  title
864   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
865   // TPC       TPCVertex             VertexerTracksNoConstraint
866   // SPD       SPDVertex             vertexer: 3D
867   // SPD       SPDVertex             vertexer: Z
868
869   if(!vtx){
870     Printf("%s:%d No AOD vertex found",(char*)__FILE__,__LINE__);
871     return kFALSE;
872   }
873
874
875   if(fDebug){
876     Printf(" n contrib %d",vtx->GetNContributors());
877     vtx->Print();
878   }
879   
880   //  if(vtx->GetNContributors()<3)return kFALSE;
881   // do not want tpc only primary vertex
882   TString vtxName(vtx->GetName());
883   if(vtxName.Contains("TPCVertex"))return kFALSE;
884
885   // no dispersion yet...
886   /* 
887   TString vtxTitle(vtx->GetTitle());
888   if(vtxTitle.Contains("vertexer: Z")){
889     if(vtx->GetDispersion()>0.02)return kFALSE;
890   }
891   */
892   return kTRUE;
893 }
894
895
896 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
897
898   if(!IsVertexValid(vtx))return kFALSE;
899
900   Float_t zvtx = vtx->GetZ();
901   Float_t yvtx = vtx->GetY();
902   Float_t xvtx = vtx->GetX();
903
904   xvtx -= fVtxXMean;
905   yvtx -= fVtxYMean;
906   zvtx -= fVtxZMean;
907
908
909
910   if(TMath::Abs(zvtx)>fVtxZCut){
911     return kFALSE;
912   }
913   fEventCutInfoESD |= kVertexZCut;  
914   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
915   if(r2>(fVtxRCut*fVtxRCut)){
916     return kFALSE;
917   }
918   fEventCutInfoESD |= kVertexRCut;  
919   return kTRUE;
920 }
921
922
923 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
924
925   if(!IsVertexValid(vtx))return kFALSE;
926
927   Float_t zvtx = vtx->GetZ();
928   Float_t yvtx = vtx->GetY();
929   Float_t xvtx = vtx->GetX();
930
931   xvtx -= fVtxXMean;
932   yvtx -= fVtxYMean;
933   zvtx -= fVtxZMean;
934
935   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
936
937   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
938   return vertexIn;
939 }
940
941 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
942   if(!esd)return kFALSE;
943   return esd->IsPileupFromSPD();
944 }
945
946 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
947   if(!esd)return kFALSE;
948   // add track cuts for which we look for cosmics...
949
950   Bool_t isCosmic = kFALSE;
951   Int_t nTracks = esd->GetNumberOfTracks();
952   Int_t nCosmicCandidates = 0;
953
954   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
955     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
956     if (!track1)  continue;
957     UInt_t status1 = track1->GetStatus();
958     //If track is ITS stand alone track, skip the track
959     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
960     if(track1->Pt()<fPtMinCosmic) continue;
961     //Start 2nd track loop to look for correlations
962     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
963       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
964       if(!track2) continue;
965       UInt_t status2 = track2->GetStatus();
966       //If track is ITS stand alone track, skip the track
967       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
968       if(track2->Pt()<fPtMinCosmic) continue;
969       //Check if back-to-back
970       Double_t mom1[3],mom2[3];
971       track1->GetPxPyPz(mom1);
972       track2->GetPxPyPz(mom2);
973       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
974       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
975       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
976       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
977
978       Float_t deltaPhi = track1->Phi()-track2->Phi();
979       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
980
981       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
982       if(rIsol<fRIsolMinCosmic) continue;
983
984       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
985         nCosmicCandidates+=1;
986         isCosmic = kTRUE;
987       }
988       
989     }
990   }
991
992   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
993
994   return isCosmic;
995 }
996
997
998 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
999
1000   Float_t cent = 0;
1001   if(fCollisionType==kPbPb){
1002     if(esd->GetCentrality()){
1003       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
1004     }
1005     if(cent<0)cent = 101;
1006     if(cent>80||cent<0)return 5;
1007     if(cent>50)return 4;
1008     if(cent>30)return 3;
1009     if(cent>10)return 2;
1010     return 1;
1011   }
1012   return 1;
1013 }
1014
1015
1016 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
1017
1018   if(fCollisionType==kPbPb){
1019     Float_t cent = aod->GetHeader()->GetCentrality();
1020     if(cent>80||cent<0)return 5;
1021     if(cent>50)return 4;
1022     if(cent>30)return 3;
1023     if(cent>10)return 2;
1024     return 1;
1025   }
1026   return 1;
1027
1028 }
1029
1030
1031 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
1032 {
1033   // Terminate analysis
1034   //
1035 }
1036
1037 Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngleVZERO(AliAODEvent *aod){
1038
1039   //  const Double_t arr_eta[11]={-3.7, -3.2, -2.7, -2.2, -1.7,0, 2.8, 3.4, 3.9, 4.5,5.1};
1040   if(!aod)return kFALSE;
1041   static bool bFirst = true;
1042   static Double_t v0phi[64] = {0,};
1043
1044   if(bFirst){
1045     int is=0;
1046     for(int iArm = 0; iArm<2; iArm++){
1047       for(int iRing = 0; iRing<4 ; iRing++){
1048         for(int iSec = 0; iSec<8 ; iSec++){
1049           v0phi[is] = 22.5 + 45. * iSec;
1050           v0phi[is] *= TMath::Pi()/180; 
1051           // cout<< is <<" "<< v0phi[is]<<endl;
1052           is++;
1053         }
1054       }
1055     }
1056     bFirst = false;
1057   }
1058
1059   // 
1060   const AliAODVZERO *aodVZERO = aod->GetVZEROData();
1061   Double_t numYZNA = 0,numXZNA = 0,sumZNA = 0;
1062   Double_t meanXA = 0,meanYA = 0;
1063
1064   Double_t numYZNC = 0,numXZNC = 0,sumZNC = 0;
1065   Double_t meanXC = 0,meanYC = 0;
1066
1067
1068
1069   static Int_t iOldRun = -1;
1070   static Int_t iFoundBin = -1;
1071
1072   if(aod->GetRunNumber()!=iOldRun&&(fp1CalibRPYA)){
1073     // search only or the bin in case of new runs
1074     iFoundBin = -1;
1075     Int_t ib = fp1CalibRPYA->FindBin(aod->GetRunNumber());
1076     Float_t err = fp1CalibRPYA->GetBinError(ib);
1077     if(err>0){// value can be zero...
1078       iFoundBin = ib;
1079     }
1080     else{
1081       Int_t ibLo = ib-1;
1082       Int_t ibUp = ib+1;
1083       while(iFoundBin<0&&(ibLo>0||ibUp<=fp1CalibRPYA->GetNbinsX())){
1084         err = fp1CalibRPYA->GetBinError(ibLo);
1085         if(err>0){
1086           iFoundBin = ibLo;
1087         }
1088         else{
1089           err = fp1CalibRPYA->GetBinError(ibUp);
1090           if(err>0)iFoundBin = ibUp;
1091         }
1092         ibUp++;
1093         ibLo--;
1094       }
1095     }
1096     iOldRun = aod->GetRunNumber();
1097   }
1098
1099   if(fDebug)Printf("%s:%d iFoundBin %d",(char*)__FILE__,__LINE__,iFoundBin);
1100
1101   if(iFoundBin>0&&(fp1CalibRPYA)){
1102     meanXA = fp1CalibRPXA->GetBinContent(iFoundBin);
1103     meanYA = fp1CalibRPYA->GetBinContent(iFoundBin);
1104     meanXC = fp1CalibRPXC->GetBinContent(iFoundBin);
1105     meanYC = fp1CalibRPYC->GetBinContent(iFoundBin);
1106   }
1107
1108   if(fDebug)Printf("%s:%d iFoundBin %1.3E %1.3E %1.3E %1.3E",(char*)__FILE__,__LINE__,meanXA,meanYA,meanXC,meanYC);
1109
1110   for (int i=0; i<64; i++) {  
1111     Double_t mult = aodVZERO->GetMultiplicity(i);
1112     Double_t phi = v0phi[i];
1113     if (mult>0) {
1114       if (i<32) { //C-side
1115         Double_t wZNC= mult;
1116         numYZNC += sin(2.*phi)*wZNC; 
1117         numXZNC += cos(2.*phi)*wZNC;
1118         sumZNC+=wZNC;
1119       }
1120       else if(i>31){ //A-side
1121         Double_t wZNA=mult;
1122         numYZNA += sin(2.*phi)*wZNA; 
1123         numXZNA += cos(2.*phi)*wZNA;
1124         sumZNA+=wZNA; 
1125       } 
1126     }// mult>0
1127   }// 64 sectors
1128
1129   
1130   Double_t XC = 0;
1131   Double_t YC = 0;
1132   Double_t XA = 0;
1133   Double_t YA = 0;
1134
1135   if(sumZNC!=0){
1136
1137     XC = numXZNC/sumZNC; 
1138     YC = numYZNC/sumZNC; 
1139     fPsiVZEROC = 0.5*TMath::ATan2(YC-meanYC, XC-meanXC);
1140     if(fPsiVZEROC>TMath::Pi()){fPsiVZEROC-=TMath::Pi();}
1141     if(fPsiVZEROC<0){fPsiVZEROC+=TMath::Pi();}
1142
1143   
1144     fh2XYC->Fill(XC-meanXC,YC-meanYC); // control
1145     fp1RPXC->Fill(aod->GetRunNumber(),XC);
1146     fp1RPYC->Fill(aod->GetRunNumber(),YC);
1147     
1148
1149   }
1150   
1151   if(sumZNA!=0){
1152     XA = numXZNA/sumZNA;
1153     YA = numYZNA/sumZNA;
1154     fPsiVZEROA = 0.5*TMath::ATan2(YA-meanYA, XA-meanXA);
1155     if(fPsiVZEROA>TMath::Pi()){fPsiVZEROA-=TMath::Pi();}
1156     if(fPsiVZEROA<0){fPsiVZEROA+=TMath::Pi();}
1157
1158     fh2XYA->Fill(XA-meanXA,YA-meanYA); // control
1159     fp1RPXA->Fill(aod->GetRunNumber(),XA);
1160     fp1RPYA->Fill(aod->GetRunNumber(),YA);
1161
1162   }
1163   return kTRUE;
1164
1165 }
1166
1167 Int_t  AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1168   Int_t iCount = 0;
1169   AliAODEvent *aod = 0;
1170   if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1171   else aod = AODEvent();
1172   if(!aod){
1173     return iCount;
1174   }
1175   for(int it = 0;it < aod->GetNumberOfTracks();++it){
1176     AliAODTrack *tr = aod->GetTrack(it);
1177     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1178     if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1179     if(tr->Pt()<fMinTrackPt)continue;
1180     list->Add(tr);
1181     iCount++;
1182   }
1183   list->Sort();
1184   return iCount;
1185
1186 }
1187