Fixes for centrality < 0, adding jet finder setting with higher pT cut
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx
1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets   
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22  
23 #include <TROOT.h>
24 #include <TRandom.h>
25 #include <TString.h>
26 #include <TSystem.h>
27 #include <TInterpreter.h>
28 #include <TChain.h>
29 #include <TFile.h>
30 #include <TKey.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3F.h>
34 #include <TProfile.h>
35 #include <TList.h>
36 #include <TLorentzVector.h>
37 #include <TClonesArray.h>
38 #include "TDatabasePDG.h"
39
40 #include "AliAnalysisTaskJetServices.h"
41 #include "AliCentrality.h"
42 #include "AliAnalysisDataContainer.h"
43 #include "AliAnalysisDataSlot.h"
44 #include "AliAnalysisManager.h"
45 #include "AliJetFinder.h"
46 #include "AliJetHeader.h"
47 #include "AliJetReader.h"
48 #include "AliJetReaderHeader.h"
49 #include "AliUA1JetHeaderV1.h"
50 #include "AliESDEvent.h"
51 #include "AliAODEvent.h"
52 #include "AliAODHandler.h"
53 #include "AliAODTrack.h"
54 #include "AliAODJet.h"
55 #include "AliAODMCParticle.h"
56 #include "AliMCEventHandler.h"
57 #include "AliMCEvent.h"
58 #include "AliStack.h"
59 #include "AliGenPythiaEventHeader.h"
60 #include "AliJetKineReaderHeader.h"
61 #include "AliGenCocktailEventHeader.h"
62 #include "AliInputEventHandler.h"
63 #include "AliPhysicsSelection.h"
64 #include "AliTriggerAnalysis.h"
65
66 #include "AliAnalysisHelperJetTasks.h"
67
68 ClassImp(AliAnalysisTaskJetServices)
69
70 AliAODHeader*  AliAnalysisTaskJetServices::fgAODHeader = NULL;
71 TClonesArray*   AliAnalysisTaskJetServices::fgAODVertices = NULL;
72
73 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
74   fUseAODInput(kFALSE),
75   fUsePhysicsSelection(kFALSE),
76   fMC(kFALSE),
77   fFilterAODCollisions(kFALSE),
78   fPhysicsSelectionFlag(AliVEvent::kMB),
79   fSelectionInfoESD(0),
80   fEventCutInfoESD(0),
81   fAvgTrials(1),
82   fVtxXMean(0),
83   fVtxYMean(0),
84   fVtxZMean(0),
85   fVtxRCut(1.),
86   fVtxZCut(8.),
87   fPtMinCosmic(5.),
88   fRIsolMinCosmic(3.),
89   fMaxCosmicAngle(0.01),
90   fNonStdFile(""),
91   fh1Xsec(0x0),
92   fh1Trials(0x0),
93   fh1PtHard(0x0),
94   fh1PtHardTrials(0x0),
95   fh1SelectionInfoESD(0x0),
96   fh1EventCutInfoESD(0),
97   fh1CentralityESD(0),
98   fh1Centrality(0),
99   fh2TriggerCount(0x0),
100   fh2ESDTriggerCount(0x0),
101   fh2TriggerVtx(0x0),
102   fh2ESDTriggerVtx(0x0),
103   fh2ESDTriggerRun(0x0),
104   fh2VtxXY(0x0),
105   fh1NCosmicsPerEvent(0x0),
106   fTriggerAnalysis(0x0),
107   fHistList(0x0)  
108 {
109   fRunRange[0] = fRunRange[1] = 0; 
110 }
111
112 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
113   AliAnalysisTaskSE(name),
114   fUseAODInput(kFALSE),
115   fUsePhysicsSelection(kFALSE),
116   fMC(kFALSE),
117   fFilterAODCollisions(kFALSE),
118   fPhysicsSelectionFlag(AliVEvent::kMB),
119   fSelectionInfoESD(0),
120   fEventCutInfoESD(0),
121   fAvgTrials(1),
122   fVtxXMean(0),
123   fVtxYMean(0),
124   fVtxZMean(0),
125   fVtxRCut(1.),
126   fVtxZCut(8.),
127   fPtMinCosmic(5.),
128   fRIsolMinCosmic(3.),
129   fMaxCosmicAngle(0.01),
130   fNonStdFile(""),
131   fh1Xsec(0x0),
132   fh1Trials(0x0),
133   fh1PtHard(0x0),
134   fh1PtHardTrials(0x0),
135   fh1SelectionInfoESD(0x0),
136   fh1EventCutInfoESD(0),
137   fh1CentralityESD(0),
138   fh1Centrality(0),
139   fh2TriggerCount(0x0),
140   fh2ESDTriggerCount(0x0),
141   fh2TriggerVtx(0x0),
142   fh2ESDTriggerVtx(0x0),
143   fh2ESDTriggerRun(0x0),
144   fh2VtxXY(0x0),
145   fh1NCosmicsPerEvent(0x0),
146   fTriggerAnalysis(0x0),
147   fHistList(0x0)  
148 {
149   fRunRange[0] = fRunRange[1] = 0; 
150   DefineOutput(1,TList::Class());
151 }
152
153
154
155 Bool_t AliAnalysisTaskJetServices::Notify()
156 {
157   //
158   // Implemented Notify() to read the cross sections
159   // and number of trials from pyxsec.root
160   // 
161
162   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
163   Float_t xsection = 0;
164   Float_t ftrials  = 1;
165
166   fAvgTrials = 1;
167   if(tree){
168     TFile *curfile = tree->GetCurrentFile();
169     if (!curfile) {
170       Error("Notify","No current file");
171       return kFALSE;
172     }
173     if(!fh1Xsec||!fh1Trials){
174       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
175       return kFALSE;
176     }
177     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
178     fh1Xsec->Fill("<#sigma>",xsection);
179     // construct a poor man average trials 
180     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
181     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
182   }  
183   return kTRUE;
184 }
185
186 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
187 {
188
189   //
190   // Create the output container
191   //
192
193
194   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
195
196   OpenFile(1);
197   if(!fHistList)fHistList = new TList();
198   fHistList->SetOwner();
199   PostData(1, fHistList); // post data in any case once
200
201   Bool_t oldStatus = TH1::AddDirectoryStatus();
202   TH1::AddDirectory(kFALSE);
203   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
204   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
205   fHistList->Add(fh1Xsec);
206
207   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
208   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
209   fHistList->Add(fh1Trials);
210
211   const Int_t nBinPt = 100;
212   Double_t binLimitsPt[nBinPt+1];
213   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
214     if(iPt == 0){
215       binLimitsPt[iPt] = 0.0;
216     }
217     else {// 1.0
218       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
219     }
220   }
221   
222
223   fh1CentralityESD = new TH1F("fh1CentralityESD","cent",102,-0.5,101.5);
224   fHistList->Add(fh1CentralityESD);
225   
226   fh1Centrality = new TH1F("fh1Centrality","cent",102,-0.5,101.5);
227   fHistList->Add(fh1Centrality);
228
229   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
230   fHistList->Add(fh2TriggerCount);
231
232   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
233   fHistList->Add(fh2ESDTriggerCount);
234   const Int_t nBins = AliAnalysisHelperJetTasks::kTrigger*kConstraints;
235   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
236   fHistList->Add(fh2TriggerVtx);
237
238   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
239   fHistList->Add(fh2ESDTriggerVtx);
240   
241
242   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
243   fHistList->Add(fh1PtHard);
244   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
245   fHistList->Add(fh1PtHardTrials);
246   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
247                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
248   fHistList->Add(fh1SelectionInfoESD);
249
250   fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
251                                 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
252   fHistList->Add(fh1EventCutInfoESD);
253
254   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
255   // 3 triggers BB BE/EB EE
256
257   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);
258   fHistList->Add(fh2ESDTriggerRun);
259
260   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
261   fHistList->Add(fh2VtxXY);
262   // =========== Switch on Sumw2 for all histos ===========
263   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
264     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
265     if (h1){
266       h1->Sumw2();
267       continue;
268     }
269     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
270     if(hn)hn->Sumw2();
271   }
272
273   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
274   fHistList->Add(fh1NCosmicsPerEvent),
275
276
277   TH1::AddDirectory(oldStatus);
278
279   // Add an AOD branch for replication
280   if(fNonStdFile.Length()){
281      if (fDebug > 1) AliInfo("Replicating header");
282      fgAODHeader = new AliAODHeader;
283      AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
284      if (fDebug > 1) AliInfo("Replicating primary vertices");
285      fgAODVertices = new TClonesArray("AliAODVertex",3);
286      fgAODVertices->SetName("vertices");
287      AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
288   }
289 }
290
291 void AliAnalysisTaskJetServices::Init()
292 {
293   //
294   // Initialization
295   //
296   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
297
298 }
299
300 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
301 {
302
303   //
304   // Execute analysis for current event
305   //
306  
307   AliAODEvent *aod = 0;
308   AliESDEvent *esd = 0;
309
310
311   
312   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
313   AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
314   fSelectionInfoESD = 0; // reset
315   fEventCutInfoESD = 0; // reset
316   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
317
318
319   static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
320    
321
322
323   if(fUseAODInput){    
324     aod = dynamic_cast<AliAODEvent*>(InputEvent());
325     if(!aod){
326       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
327       return;
328     }
329     // fethc the header
330   }
331   else{
332     //  assume that the AOD is in the general output...
333     aod  = AODEvent();
334     if(!aod){
335       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
336       return;
337     }
338     esd = dynamic_cast<AliESDEvent*>(InputEvent());
339   }
340   if(aod&&fDebug>2){
341     aod->Print();
342     Printf("Vertices %d",aod->GetNumberOfVertices());
343     Printf("tracks %d",aod->GetNumberOfTracks());
344     Printf("jets %d",aod->GetNJets());
345   }
346   fSelectionInfoESD |= kNoEventCut;
347   fEventCutInfoESD |= kNoEventCut;
348
349   Bool_t esdVtxValid = false;
350   Bool_t esdVtxIn = false;
351   Bool_t aodVtxValid = false;
352   Bool_t aodVtxIn = false;
353
354   if(esd){
355     // trigger analyisis
356     if(!fTriggerAnalysis){
357       fTriggerAnalysis = new AliTriggerAnalysis;
358       fTriggerAnalysis->SetAnalyzeMC(fMC);
359       fTriggerAnalysis->SetSPDGFOThreshhold(1);
360     }
361     //    fTriggerAnalysis->FillTriggerClasses(esd);
362     Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
363     Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
364     Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG);
365     Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG);
366     Bool_t spdFO      = fTriggerAnalysis->SPDFiredChips(esd, 0);;
367     if(v0A)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0A;
368     if(v0C)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0C;
369     if(!(v0ABG||v0CBG))fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kNoV0BG;
370     if(spdFO)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kSPDFO;
371   }
372  
373   // Apply additional constraints
374   Bool_t esdEventSelected = IsEventSelected(esd);
375   Bool_t esdEventPileUp = IsEventPileUp(esd);
376   Bool_t esdEventCosmic = IsEventCosmic(esd);
377
378   Bool_t aodEventSelected = IsEventSelected(aod);
379
380   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
381
382
383   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
384   fh1EventCutInfoESD->Fill(fEventCutInfoESD);
385
386   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
387   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
388   if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
389   if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
390
391
392   // here we have all selection information, fill histogram
393   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
394     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
395   }
396   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); 
397
398   if(esd&&aod&&fDebug){
399     if(esdEventSelected&&!aodEventSelected){
400       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
401       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
402       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
403       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
404       vtxESD->Print();
405       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
406       vtxAOD->Print();
407     }
408   }
409
410   // loop over all possible triggers for esd 
411
412   Float_t cent = 100;
413   if(aod)cent = aod->GetHeader()->GetCentrality();
414   if(cent<0)cent = 101;
415   if(esd){
416     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
417     esdVtxValid = IsVertexValid(vtxESD);
418     esdVtxIn = IsVertexIn(vtxESD);
419     if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
420       if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
421       if(cent<=80&&esdVtxIn){
422         aodH->SetFillAOD(kTRUE);
423         aodH->SetFillExtension(kTRUE);
424       }
425     }
426
427
428     Float_t zvtx = vtxESD->GetZ();
429     Int_t  iCl = GetEventClass(esd);
430     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
431     Bool_t cand = physicsSelection;
432
433     if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
434     fh2ESDTriggerCount->Fill(0.,kAllTriggered); 
435     fh2ESDTriggerCount->Fill(iCl,kAllTriggered); 
436     if(cand){
437       fh2ESDTriggerCount->Fill(0.,kSelectedALICE); 
438       fh2ESDTriggerCount->Fill(iCl,kSelectedALICE); 
439       fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
440     }
441     //    if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
442     if(esdVtxValid){
443       fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
444       fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
445       fh2ESDTriggerVtx->Fill(iCl,zvtx);
446       if(esdVtxIn){
447         fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
448         fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
449         fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
450       }
451       if(cand){
452         fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
453         fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
454         fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
455       }
456     }
457
458     if(cand&&esdVtxIn&&iCl<5){
459       fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
460       fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
461       fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
462       fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
463       fh2ESDTriggerCount->Fill(iCl,kSelected);
464       fh2ESDTriggerCount->Fill(0.,kSelected);
465       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
466       if(esd->GetCentrality()){
467         Float_t tmpCent = 100;
468         tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
469         if(tmpCent<0)tmpCent = 101;
470         fh1CentralityESD->Fill(tmpCent);
471       }
472     }
473   }
474
475
476
477   if(aod){
478     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
479     aodVtxValid = IsVertexValid(vtxAOD);
480     aodVtxIn = IsVertexIn(vtxAOD);
481     Float_t zvtx = vtxAOD->GetZ();
482     Int_t  iCl = GetEventClass(aod);
483     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
484     Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
485     if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
486     fh2TriggerCount->Fill(0.,kAllTriggered); 
487     fh2TriggerCount->Fill(iCl,kAllTriggered); 
488     if(cand){
489       fh2TriggerCount->Fill(0.,kSelectedALICE); 
490       fh2TriggerCount->Fill(iCl,kSelectedALICE); 
491       fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
492     }
493     if(aodVtxValid){
494       fh2TriggerCount->Fill(0.,kTriggeredVertex);
495       fh2TriggerCount->Fill(iCl,kTriggeredVertex);
496       fh2TriggerVtx->Fill(iCl,zvtx);
497       if(aodVtxIn){
498         fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
499         fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
500         fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
501       }
502       if(cand){
503         fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
504         fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
505         fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
506       }
507     }
508     if(cand&&aodVtxIn&&iCl<5){
509       fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
510       fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
511       fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
512       fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
513       fh2TriggerCount->Fill(iCl,kSelected);
514       fh2TriggerCount->Fill(0.,kSelected);
515       fh1Centrality->Fill(cent);
516       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
517       if(fUseAODInput&&cent<=80){
518         if(fFilterAODCollisions&&aod){
519           aodH->SetFillAOD(kTRUE);
520         }
521       }
522     }    
523   }
524
525   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
526
527   
528   Double_t ptHard = 0; 
529   Double_t nTrials = 1; // Trials for MC trigger 
530
531   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
532   AliMCEvent* mcEvent = MCEvent();
533   //    AliStack *pStack = 0; 
534   if(mcEvent){
535     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
536     if(pythiaGenHeader){
537       nTrials = pythiaGenHeader->Trials();
538       ptHard  = pythiaGenHeader->GetPtHard();
539       int iProcessType = pythiaGenHeader->ProcessType();
540       // 11 f+f -> f+f
541       // 12 f+barf -> f+barf
542       // 13 f+barf -> g+g
543       // 28 f+g -> f+g
544       // 53 g+g -> f+barf
545       // 68 g+g -> g+g
546       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
547       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
548       fh1PtHard->Fill(ptHard);
549       fh1PtHardTrials->Fill(ptHard,nTrials);
550
551     }// if pythia gen header
552   }
553
554   // trigger selection
555   
556   // replication of 
557   if(fNonStdFile.Length()&&aod){
558     if (fgAODHeader){
559       *fgAODHeader =  *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
560     }
561     if(fgAODVertices){
562       fgAODVertices->Delete();
563       TClonesArray &vertices = *fgAODVertices;
564       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
565       // we only use some basic information, 
566       
567
568       Double_t pos[3];
569       Double_t covVtx[6];
570
571       vtxAOD->GetXYZ(pos); // position                                                                
572       vtxAOD->GetCovMatrix(covVtx); //covariance matrix                                                      
573       Int_t jVertices = 0;
574       AliAODVertex * vtx = new(vertices[jVertices++])
575         AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
576       vtx->SetName(vtxAOD->GetName());
577       vtx->SetTitle(vtxAOD->GetTitle());
578       TString vtitle = vtxAOD->GetTitle();
579       vtx->SetNContributors(vtxAOD->GetNContributors());
580
581       // Add SPD "main" vertex                                                                    
582       const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
583       vtxS->GetXYZ(pos); // position
584       vtxS->GetCovMatrix(covVtx); //covariance matrix                                                      
585       AliAODVertex * mVSPD = new(vertices[jVertices++])
586         AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
587       mVSPD->SetName(vtxS->GetName());
588       mVSPD->SetTitle(vtxS->GetTitle());
589       mVSPD->SetNContributors(vtxS->GetNContributors());
590       
591       // Add tpc only vertex
592       if(esd){
593         const AliESDVertex *vtxT =  esd->GetPrimaryVertexTPC();
594         vtxT->GetXYZ(pos); // position                                                    
595         vtxT->GetCovMatrix(covVtx); //covariance matrix                                               
596         AliAODVertex * mVTPC = new(vertices[jVertices++])
597           AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
598         mVTPC->SetName(vtxT->GetName());
599         mVTPC->SetTitle(vtxT->GetTitle());
600         mVTPC->SetNContributors(vtxT->GetNContributors());
601       }
602     }
603   }
604   
605   PostData(1, fHistList);
606 }
607
608 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
609   if(!esd)return kFALSE;
610   const AliESDVertex *vtx = esd->GetPrimaryVertex();
611   return IsVertexIn(vtx); // vertex in calls vertex valid
612 }
613
614 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
615   if(fgAODVertices){
616     fgAODVertices->Delete();
617     delete fgAODVertices;
618   }
619   if(fgAODHeader)delete fgAODHeader;
620 }
621
622
623 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
624   if(!aod)return kFALSE;
625   const AliAODVertex *vtx = aod->GetPrimaryVertex();
626   return IsVertexIn(vtx); // VertexIn calls VertexValid
627 }
628
629 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
630
631   // We can check the type of the vertex by its name and title
632   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
633   // vtx       name                  title
634   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
635   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
636   // TPC       TPCVertex             VertexerTracksNoConstraint
637   // TPC       TPCVertex             VertexerTracksWithConstraint
638   // SPD       SPDVertex             vertexer: 3D
639   // SPD       SPDVertex             vertexer: Z
640   
641   Int_t nCont = vtx->GetNContributors();
642   if(nCont>=1){
643     fEventCutInfoESD |= kContributorsCut1;    
644     if(nCont>=2){
645       fEventCutInfoESD |= kContributorsCut2;    
646       if(nCont>=3){
647         fEventCutInfoESD |= kContributorsCut3;    
648       }
649     }
650   }
651   
652   if(nCont<3)return kFALSE;
653
654   // do not want tpc only primary vertex
655   TString vtxName(vtx->GetName());
656   if(vtxName.Contains("TPCVertex")){
657     fEventCutInfoESD |= kVertexTPC;
658     return kFALSE;
659   }
660   if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
661   if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
662
663
664   TString vtxTitle(vtx->GetTitle());
665   if(vtxTitle.Contains("vertexer: Z")){
666     if(vtx->GetDispersion()>0.02)return kFALSE;   
667   }
668   fEventCutInfoESD |= kSPDDispersionCut;
669   return kTRUE;
670 }
671
672
673 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
674
675   // We can check the type of the vertex by its name and title
676   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
677   // vtx       name                  title
678   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
679   // TPC       TPCVertex             VertexerTracksNoConstraint
680   // SPD       SPDVertex             vertexer: 3D
681   // SPD       SPDVertex             vertexer: Z
682
683   if(fDebug){
684     Printf(" n contrib %d",vtx->GetNContributors());
685     vtx->Print();
686   }
687   
688   //  if(vtx->GetNContributors()<3)return kFALSE;
689   // do not want tpc only primary vertex
690   TString vtxName(vtx->GetName());
691   if(vtxName.Contains("TPCVertex"))return kFALSE;
692
693   // no dispersion yet...
694   /* 
695   TString vtxTitle(vtx->GetTitle());
696   if(vtxTitle.Contains("vertexer: Z")){
697     if(vtx->GetDispersion()>0.02)return kFALSE;
698   }
699   */
700   return kTRUE;
701 }
702
703
704 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
705
706   if(!IsVertexValid(vtx))return kFALSE;
707
708   Float_t zvtx = vtx->GetZ();
709   Float_t yvtx = vtx->GetY();
710   Float_t xvtx = vtx->GetX();
711
712   xvtx -= fVtxXMean;
713   yvtx -= fVtxYMean;
714   zvtx -= fVtxZMean;
715
716
717
718   if(TMath::Abs(zvtx)>fVtxZCut){
719     return kFALSE;
720   }
721   fEventCutInfoESD |= kVertexZCut;  
722   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
723   if(r2>(fVtxRCut*fVtxRCut)){
724     return kFALSE;
725   }
726   fEventCutInfoESD |= kVertexRCut;  
727   return kTRUE;
728 }
729
730
731 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
732
733   if(!IsVertexValid(vtx))return kFALSE;
734
735   Float_t zvtx = vtx->GetZ();
736   Float_t yvtx = vtx->GetY();
737   Float_t xvtx = vtx->GetX();
738
739   xvtx -= fVtxXMean;
740   yvtx -= fVtxYMean;
741   zvtx -= fVtxZMean;
742
743   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
744
745   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
746   return vertexIn;
747 }
748
749 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
750   if(!esd)return kFALSE;
751   return esd->IsPileupFromSPD();
752 }
753
754 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
755   if(!esd)return kFALSE;
756   // add track cuts for which we look for cosmics...
757
758   Bool_t isCosmic = kFALSE;
759   Int_t nTracks = esd->GetNumberOfTracks();
760   Int_t nCosmicCandidates = 0;
761
762   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
763     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
764     if (!track1)  continue;
765     UInt_t status1 = track1->GetStatus();
766     //If track is ITS stand alone track, skip the track
767     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
768     if(track1->Pt()<fPtMinCosmic) continue;
769     //Start 2nd track loop to look for correlations
770     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
771       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
772       if(!track2) continue;
773       UInt_t status2 = track2->GetStatus();
774       //If track is ITS stand alone track, skip the track
775       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
776       if(track2->Pt()<fPtMinCosmic) continue;
777       //Check if back-to-back
778       Double_t mom1[3],mom2[3];
779       track1->GetPxPyPz(mom1);
780       track2->GetPxPyPz(mom2);
781       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
782       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
783       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
784       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
785
786       Float_t deltaPhi = track1->Phi()-track2->Phi();
787       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
788
789       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
790       if(rIsol<fRIsolMinCosmic) continue;
791
792       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
793         nCosmicCandidates+=1;
794         isCosmic = kTRUE;
795       }
796       
797     }
798   }
799
800   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
801
802   return isCosmic;
803 }
804
805
806 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
807
808   Float_t cent = 999;
809   if(esd->GetCentrality()){
810     cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
811   }
812   if(cent>80||cent<0)return 5;
813   if(cent>50)return 4;
814   if(cent>30)return 3;
815   if(cent>10)return 2;
816   return 1;
817
818 }
819
820
821 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
822
823   Float_t cent = aod->GetHeader()->GetCentrality();
824   if(cent>80||cent<0)return 5;
825   if(cent>50)return 4;
826   if(cent>30)return 3;
827   if(cent>10)return 2;
828   return 1;
829
830 }
831
832
833 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
834 {
835   // Terminate analysis
836   //
837 }
838