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