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