38cfd1c4633ddf7b34d3a37c4be4f664b20ab867
[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 "AliESDCentrality.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   fh2TriggerCount(0x0),
98   fh2ESDTriggerCount(0x0),
99   fh2TriggerVtx(0x0),
100   fh2ESDTriggerVtx(0x0),
101   fh2ESDTriggerRun(0x0),
102   fh2VtxXY(0x0),
103   fh1NCosmicsPerEvent(0x0),
104   fTriggerAnalysis(0x0),
105   fHistList(0x0)  
106 {
107   fRunRange[0] = fRunRange[1] = 0; 
108 }
109
110 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
111   AliAnalysisTaskSE(name),
112   fUseAODInput(kFALSE),
113   fUsePhysicsSelection(kFALSE),
114   fMC(kFALSE),
115   fFilterAODCollisions(kFALSE),
116   fPhysicsSelectionFlag(AliVEvent::kMB),
117   fSelectionInfoESD(0),
118   fEventCutInfoESD(0),
119   fAvgTrials(1),
120   fVtxXMean(0),
121   fVtxYMean(0),
122   fVtxZMean(0),
123   fVtxRCut(1.),
124   fVtxZCut(8.),
125   fPtMinCosmic(5.),
126   fRIsolMinCosmic(3.),
127   fMaxCosmicAngle(0.01),
128   fNonStdFile(""),
129   fh1Xsec(0x0),
130   fh1Trials(0x0),
131   fh1PtHard(0x0),
132   fh1PtHardTrials(0x0),
133   fh1SelectionInfoESD(0x0),
134   fh1EventCutInfoESD(0),
135   fh2TriggerCount(0x0),
136   fh2ESDTriggerCount(0x0),
137   fh2TriggerVtx(0x0),
138   fh2ESDTriggerVtx(0x0),
139   fh2ESDTriggerRun(0x0),
140   fh2VtxXY(0x0),
141   fh1NCosmicsPerEvent(0x0),
142   fTriggerAnalysis(0x0),
143   fHistList(0x0)  
144 {
145   fRunRange[0] = fRunRange[1] = 0; 
146   DefineOutput(1,TList::Class());
147 }
148
149
150
151 Bool_t AliAnalysisTaskJetServices::Notify()
152 {
153   //
154   // Implemented Notify() to read the cross sections
155   // and number of trials from pyxsec.root
156   // 
157
158   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
159   Float_t xsection = 0;
160   Float_t ftrials  = 1;
161
162   fAvgTrials = 1;
163   if(tree){
164     TFile *curfile = tree->GetCurrentFile();
165     if (!curfile) {
166       Error("Notify","No current file");
167       return kFALSE;
168     }
169     if(!fh1Xsec||!fh1Trials){
170       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
171       return kFALSE;
172     }
173     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
174     fh1Xsec->Fill("<#sigma>",xsection);
175     // construct a poor man average trials 
176     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
177     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
178   }  
179   return kTRUE;
180 }
181
182 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
183 {
184
185   //
186   // Create the output container
187   //
188
189
190   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
191
192   OpenFile(1);
193   if(!fHistList)fHistList = new TList();
194
195   Bool_t oldStatus = TH1::AddDirectoryStatus();
196   TH1::AddDirectory(kFALSE);
197   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
198   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
199   fHistList->Add(fh1Xsec);
200
201   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
202   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
203   fHistList->Add(fh1Trials);
204
205   const Int_t nBinPt = 100;
206   Double_t binLimitsPt[nBinPt+1];
207   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
208     if(iPt == 0){
209       binLimitsPt[iPt] = 0.0;
210     }
211     else {// 1.0
212       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
213     }
214   }
215   
216   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
217   fHistList->Add(fh2TriggerCount);
218
219   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
220   fHistList->Add(fh2ESDTriggerCount);
221   const Int_t nBins = AliAnalysisHelperJetTasks::kTrigger*kConstraints;
222   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
223   fHistList->Add(fh2TriggerVtx);
224
225   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
226   fHistList->Add(fh2ESDTriggerVtx);
227   
228
229   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
230   fHistList->Add(fh1PtHard);
231   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
232   fHistList->Add(fh1PtHardTrials);
233   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
234                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
235   fHistList->Add(fh1SelectionInfoESD);
236
237   fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
238                                 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
239   fHistList->Add(fh1EventCutInfoESD);
240
241   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
242   // 3 triggers BB BE/EB EE
243
244   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);
245   fHistList->Add(fh2ESDTriggerRun);
246
247   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
248   fHistList->Add(fh2VtxXY);
249   // =========== Switch on Sumw2 for all histos ===========
250   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
251     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
252     if (h1){
253       h1->Sumw2();
254       continue;
255     }
256     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
257     if(hn)hn->Sumw2();
258   }
259
260   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
261   fHistList->Add(fh1NCosmicsPerEvent),
262
263
264   TH1::AddDirectory(oldStatus);
265
266   // Add an AOD branch for replication
267   if(fNonStdFile.Length()){
268      if (fDebug > 1) AliInfo("Replicating header");
269      fgAODHeader = new AliAODHeader;
270      AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
271      /*
272      if (fDebug > 1) AliInfo("Replicating vertices");
273      fgAODVertices = new TClonesArray("AliAODVertex",2);
274      fgAODVertices->SetName("vertices");
275      AddAODBranch("AliAODHeader",&fgAODVertices,fNonStdFile.Data());
276      */
277   }
278 }
279
280 void AliAnalysisTaskJetServices::Init()
281 {
282   //
283   // Initialization
284   //
285   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
286
287 }
288
289 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
290 {
291
292   //
293   // Execute analysis for current event
294   //
295  
296   AliAODEvent *aod = 0;
297   AliESDEvent *esd = 0;
298
299
300   
301   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
302   AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
303   fSelectionInfoESD = 0; // reset
304   fEventCutInfoESD = 0; // reset
305   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
306
307
308   static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
309    
310
311
312   if(fUseAODInput){    
313     aod = dynamic_cast<AliAODEvent*>(InputEvent());
314     if(!aod){
315       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
316       return;
317     }
318     // fethc the header
319   }
320   else{
321     //  assume that the AOD is in the general output...
322     aod  = AODEvent();
323     if(!aod){
324       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
325       return;
326     }
327     esd = dynamic_cast<AliESDEvent*>(InputEvent());
328   }
329   if(aod&&fDebug>2){
330     aod->Print();
331     Printf("Vertices %d",aod->GetNumberOfVertices());
332     Printf("tracks %d",aod->GetNumberOfTracks());
333     Printf("jets %d",aod->GetNJets());
334   }
335   fSelectionInfoESD |= kNoEventCut;
336   fEventCutInfoESD |= kNoEventCut;
337
338   Bool_t esdVtxValid = false;
339   Bool_t esdVtxIn = false;
340   Bool_t aodVtxValid = false;
341   Bool_t aodVtxIn = false;
342
343   if(esd){
344     // trigger analyisis
345     if(!fTriggerAnalysis){
346       fTriggerAnalysis = new AliTriggerAnalysis;
347       fTriggerAnalysis->SetAnalyzeMC(fMC);
348       fTriggerAnalysis->SetSPDGFOThreshhold(1);
349     }
350     //    fTriggerAnalysis->FillTriggerClasses(esd);
351     Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
352     Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
353     Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG);
354     Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG);
355     Bool_t spdFO      = fTriggerAnalysis->SPDFiredChips(esd, 0);;
356     if(v0A)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0A;
357     if(v0C)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0C;
358     if(!(v0ABG||v0CBG))fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kNoV0BG;
359     if(spdFO)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kSPDFO;
360   }
361  
362   // Apply additional constraints
363   Bool_t esdEventSelected = IsEventSelected(esd);
364   Bool_t esdEventPileUp = IsEventPileUp(esd);
365   Bool_t esdEventCosmic = IsEventCosmic(esd);
366
367   Bool_t aodEventSelected = IsEventSelected(aod);
368
369   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
370   if(aodH&&physicsSelection&&fFilterAODCollisions&&fgAODHeader){
371     if(fgAODHeader->GetCentrality()<=80){
372       aodH->SetFillAOD(kTRUE);
373     }
374   }
375
376
377   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
378   fh1EventCutInfoESD->Fill(fEventCutInfoESD);
379
380   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
381   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
382   if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
383   if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
384
385
386   // here we have all selection information, fill histogram
387   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
388     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
389   }
390   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); 
391
392   if(esd&&aod&&fDebug){
393     if(esdEventSelected&&!aodEventSelected){
394       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
395       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
396       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
397       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
398       vtxESD->Print();
399       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
400       vtxAOD->Print();
401     }
402   }
403
404   // loop over all possible triggers for esd 
405
406   if(esd){
407     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
408     esdVtxValid = IsVertexValid(vtxESD);
409     esdVtxIn = IsVertexIn(vtxESD);
410     Float_t zvtx = vtxESD->GetZ();
411     Int_t  iCl = GetEventClass(esd);
412     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
413     Bool_t cand = physicsSelection;
414
415     if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
416
417     if(cand){
418       fh2ESDTriggerCount->Fill(0.,kSelectedALICE); 
419       fh2ESDTriggerCount->Fill(iCl,kSelectedALICE); 
420       fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
421     }
422     //    if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
423     if(esdVtxValid){
424       fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
425       fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
426       fh2ESDTriggerVtx->Fill(iCl,zvtx);
427       if(esdVtxIn){
428         fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
429         fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
430         fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
431       }
432       if(cand){
433         fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
434         fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
435         fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
436       }
437     }
438     if(cand&&esdVtxIn){
439       fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
440       fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
441       fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
442       fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
443       fh2ESDTriggerCount->Fill(iCl,kSelected);
444       fh2ESDTriggerCount->Fill(0.,kSelected);
445       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
446     }
447   }
448
449
450
451   if(aod){
452     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
453     aodVtxValid = IsVertexValid(vtxAOD);
454     aodVtxIn = IsVertexIn(vtxAOD);
455     Float_t zvtx = vtxAOD->GetZ();
456     Int_t  iCl = GetEventClass(aod);
457     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
458     Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
459     if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
460     if(cand){
461       fh2TriggerCount->Fill(0.,kSelectedALICE); 
462       fh2TriggerCount->Fill(iCl,kSelectedALICE); 
463       fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
464     }
465     if(aodVtxValid){
466       fh2TriggerCount->Fill(0.,kTriggeredVertex);
467       fh2TriggerCount->Fill(iCl,kTriggeredVertex);
468       fh2TriggerVtx->Fill(iCl,zvtx);
469       if(aodVtxIn){
470         fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
471         fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
472         fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
473       }
474       if(cand){
475         fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
476         fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
477         fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
478       }
479     }
480     if(cand&&aodVtxIn){
481       fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
482       fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
483       fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
484       fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
485       fh2TriggerCount->Fill(iCl,kSelected);
486       fh2TriggerCount->Fill(0.,kSelected);
487       if(fUseAODInput)AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
488     }
489   }
490
491   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
492
493   
494   Double_t ptHard = 0; 
495   Double_t nTrials = 1; // Trials for MC trigger 
496
497   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
498   AliMCEvent* mcEvent = MCEvent();
499   //    AliStack *pStack = 0; 
500   if(mcEvent){
501     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
502     if(pythiaGenHeader){
503       nTrials = pythiaGenHeader->Trials();
504       ptHard  = pythiaGenHeader->GetPtHard();
505       int iProcessType = pythiaGenHeader->ProcessType();
506       // 11 f+f -> f+f
507       // 12 f+barf -> f+barf
508       // 13 f+barf -> g+g
509       // 28 f+g -> f+g
510       // 53 g+g -> f+barf
511       // 68 g+g -> g+g
512       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
513       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
514       fh1PtHard->Fill(ptHard);
515       fh1PtHardTrials->Fill(ptHard,nTrials);
516
517     }// if pythia gen header
518   }
519
520   // trigger selection
521   
522   // replication of 
523   if(fNonStdFile.Length()&&aod){
524     if (fgAODHeader){
525       *fgAODHeader =  *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
526     }
527     if(fgAODVertices){
528       Printf("%p",fgAODVertices);
529       Printf("%p",fgAODVertices->At(0));
530       fgAODVertices->Delete();
531       TClonesArray &vertices = *fgAODVertices;
532
533       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
534       AliAODVertex *primary = new (&vertices[0]) AliAODVertex();
535       primary->SetName(vtxAOD->GetName());
536       primary->SetTitle(vtxAOD->GetTitle());
537     }
538   }
539   
540   PostData(1, fHistList);
541 }
542
543 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
544   if(!esd)return kFALSE;
545   const AliESDVertex *vtx = esd->GetPrimaryVertex();
546   return IsVertexIn(vtx); // vertex in calls vertex valid
547 }
548
549 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
550   if(fgAODVertices){
551     fgAODVertices->Delete();
552     delete fgAODVertices;
553   }
554 }
555
556
557 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
558   if(!aod)return kFALSE;
559   const AliAODVertex *vtx = aod->GetPrimaryVertex();
560   return IsVertexIn(vtx); // VertexIn calls VertexValid
561 }
562
563 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
564
565   // We can check the type of the vertex by its name and title
566   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
567   // vtx       name                  title
568   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
569   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
570   // TPC       TPCVertex             VertexerTracksNoConstraint
571   // TPC       TPCVertex             VertexerTracksWithConstraint
572   // SPD       SPDVertex             vertexer: 3D
573   // SPD       SPDVertex             vertexer: Z
574   
575   Int_t nCont = vtx->GetNContributors();
576
577   if(nCont>=1){
578     fEventCutInfoESD |= kContributorsCut1;    
579     if(nCont>=2){
580       fEventCutInfoESD |= kContributorsCut2;    
581       if(nCont>=3){
582         fEventCutInfoESD |= kContributorsCut3;    
583       }
584     }
585   }
586   
587   if(nCont<3)return kFALSE;
588
589   // do not want tpc only primary vertex
590   TString vtxName(vtx->GetName());
591   if(vtxName.Contains("TPCVertex")){
592     fEventCutInfoESD |= kVertexTPC;
593     return kFALSE;
594   }
595   if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
596   if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
597
598
599   TString vtxTitle(vtx->GetTitle());
600   if(vtxTitle.Contains("vertexer: Z")){
601     if(vtx->GetDispersion()>0.02)return kFALSE;   
602   }
603   fEventCutInfoESD |= kSPDDispersionCut;
604   return kTRUE;
605 }
606
607
608 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
609
610   // We can check the type of the vertex by its name and title
611   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
612   // vtx       name                  title
613   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
614   // TPC       TPCVertex             VertexerTracksNoConstraint
615   // SPD       SPDVertex             vertexer: 3D
616   // SPD       SPDVertex             vertexer: Z
617   
618   if(vtx->GetNContributors()<3)return kFALSE;
619   // do not want tpc only primary vertex
620   TString vtxName(vtx->GetName());
621   if(vtxName.Contains("TPCVertex"))return kFALSE;
622
623   // no dispersion yet...
624   /* 
625   TString vtxTitle(vtx->GetTitle());
626   if(vtxTitle.Contains("vertexer: Z")){
627     if(vtx->GetDispersion()>0.02)return kFALSE;
628   }
629   */
630   return kTRUE;
631 }
632
633
634 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
635
636   if(!IsVertexValid(vtx))return kFALSE;
637
638   Float_t zvtx = vtx->GetZ();
639   Float_t yvtx = vtx->GetY();
640   Float_t xvtx = vtx->GetX();
641
642   xvtx -= fVtxXMean;
643   yvtx -= fVtxYMean;
644   zvtx -= fVtxZMean;
645
646
647
648   if(TMath::Abs(zvtx)>fVtxZCut){
649     return kFALSE;
650   }
651   fEventCutInfoESD |= kVertexZCut;  
652   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
653   if(r2>(fVtxRCut*fVtxRCut)){
654     return kFALSE;
655   }
656   fEventCutInfoESD |= kVertexRCut;  
657   return kTRUE;
658 }
659
660
661 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
662
663   if(!IsVertexValid(vtx))return kFALSE;
664
665   Float_t zvtx = vtx->GetZ();
666   Float_t yvtx = vtx->GetY();
667   Float_t xvtx = vtx->GetX();
668
669   xvtx -= fVtxXMean;
670   yvtx -= fVtxYMean;
671   zvtx -= fVtxZMean;
672
673   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
674
675   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
676   return vertexIn;
677 }
678
679 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
680   if(!esd)return kFALSE;
681   return esd->IsPileupFromSPD();
682 }
683
684 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
685   if(!esd)return kFALSE;
686   // add track cuts for which we look for cosmics...
687
688   Bool_t isCosmic = kFALSE;
689   Int_t nTracks = esd->GetNumberOfTracks();
690   Int_t nCosmicCandidates = 0;
691
692   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
693     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
694     if (!track1)  continue;
695     UInt_t status1 = track1->GetStatus();
696     //If track is ITS stand alone track, skip the track
697     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
698     if(track1->Pt()<fPtMinCosmic) continue;
699     //Start 2nd track loop to look for correlations
700     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
701       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
702       if(!track2) continue;
703       UInt_t status2 = track2->GetStatus();
704       //If track is ITS stand alone track, skip the track
705       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
706       if(track2->Pt()<fPtMinCosmic) continue;
707       //Check if back-to-back
708       Double_t mom1[3],mom2[3];
709       track1->GetPxPyPz(mom1);
710       track2->GetPxPyPz(mom2);
711       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
712       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
713       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
714       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
715
716       Float_t deltaPhi = track1->Phi()-track2->Phi();
717       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
718
719       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
720       if(rIsol<fRIsolMinCosmic) continue;
721
722       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
723         nCosmicCandidates+=1;
724         isCosmic = kTRUE;
725       }
726       
727     }
728   }
729
730   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
731
732   return isCosmic;
733 }
734
735
736 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
737
738   Float_t cent = 999;
739   if(esd->GetCentrality()){
740     cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
741   }
742   if(cent>50)return 4;
743   if(cent>30)return 3;
744   if(cent>10)return 2;
745   return 1;
746
747 }
748
749
750 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
751
752   Float_t cent = aod->GetHeader()->GetCentrality();
753
754   if(cent>50)return 4;
755   if(cent>30)return 3;
756   if(cent>10)return 2;
757   return 1;
758
759 }
760
761
762 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
763 {
764   // Terminate analysis
765   //
766 }
767