updates for differnt jet configs (including background), reduced number of UA1 Jets
[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      if (fDebug > 1) AliInfo("Replicating primary vertices");
272
273      fgAODVertices = new TClonesArray("AliAODVertex",3);
274      fgAODVertices->SetName("vertices");
275      AddAODBranch("TClonesArray",&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       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
533       // we only use some basic information, 
534       
535
536       Double_t pos[3];
537       Double_t covVtx[6];
538
539       vtxAOD->GetXYZ(pos); // position                                                                
540       vtxAOD->GetCovMatrix(covVtx); //covariance matrix                                                      
541       Int_t jVertices = 0;
542       AliAODVertex * vtx = new(vertices[jVertices++])
543         AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
544       vtx->SetName(vtxAOD->GetName());
545       vtx->SetTitle(vtxAOD->GetTitle());
546
547       TString vtitle = vtxAOD->GetTitle();
548       if (!vtitle.Contains("VertexerTracks"))
549         vtx->SetNContributors(vtxAOD->GetNContributors());
550
551       // Add SPD "main" vertex                                                                    
552       const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
553       vtxS->GetXYZ(pos); // position
554       vtxS->GetCovMatrix(covVtx); //covariance matrix                                                      
555       AliAODVertex * mVSPD = new(vertices[jVertices++])
556         AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
557       mVSPD->SetName(vtxS->GetName());
558       mVSPD->SetTitle(vtxS->GetTitle());
559       mVSPD->SetNContributors(vtxS->GetNContributors());
560       
561       // Add tpc only vertex
562       if(esd){
563         const AliESDVertex *vtxT =  esd->GetPrimaryVertexTPC();
564         vtxT->GetXYZ(pos); // position                                                    
565         vtxT->GetCovMatrix(covVtx); //covariance matrix                                               
566         AliAODVertex * mVTPC = new(vertices[jVertices++])
567           AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
568         mVTPC->SetName(vtxT->GetName());
569         mVTPC->SetTitle(vtxT->GetTitle());
570         mVTPC->SetNContributors(vtxT->GetNContributors());
571       }
572     }
573   }
574   
575   PostData(1, fHistList);
576 }
577
578 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
579   if(!esd)return kFALSE;
580   const AliESDVertex *vtx = esd->GetPrimaryVertex();
581   return IsVertexIn(vtx); // vertex in calls vertex valid
582 }
583
584 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
585   if(fgAODVertices){
586     fgAODVertices->Delete();
587     delete fgAODVertices;
588   }
589   if(fgAODHeader)delete fgAODHeader;
590 }
591
592
593 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
594   if(!aod)return kFALSE;
595   const AliAODVertex *vtx = aod->GetPrimaryVertex();
596   return IsVertexIn(vtx); // VertexIn calls VertexValid
597 }
598
599 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
600
601   // We can check the type of the vertex by its name and title
602   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
603   // vtx       name                  title
604   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
605   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
606   // TPC       TPCVertex             VertexerTracksNoConstraint
607   // TPC       TPCVertex             VertexerTracksWithConstraint
608   // SPD       SPDVertex             vertexer: 3D
609   // SPD       SPDVertex             vertexer: Z
610   
611   Int_t nCont = vtx->GetNContributors();
612
613   if(nCont>=1){
614     fEventCutInfoESD |= kContributorsCut1;    
615     if(nCont>=2){
616       fEventCutInfoESD |= kContributorsCut2;    
617       if(nCont>=3){
618         fEventCutInfoESD |= kContributorsCut3;    
619       }
620     }
621   }
622   
623   if(nCont<3)return kFALSE;
624
625   // do not want tpc only primary vertex
626   TString vtxName(vtx->GetName());
627   if(vtxName.Contains("TPCVertex")){
628     fEventCutInfoESD |= kVertexTPC;
629     return kFALSE;
630   }
631   if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
632   if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
633
634
635   TString vtxTitle(vtx->GetTitle());
636   if(vtxTitle.Contains("vertexer: Z")){
637     if(vtx->GetDispersion()>0.02)return kFALSE;   
638   }
639   fEventCutInfoESD |= kSPDDispersionCut;
640   return kTRUE;
641 }
642
643
644 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
645
646   // We can check the type of the vertex by its name and title
647   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
648   // vtx       name                  title
649   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
650   // TPC       TPCVertex             VertexerTracksNoConstraint
651   // SPD       SPDVertex             vertexer: 3D
652   // SPD       SPDVertex             vertexer: Z
653   
654   if(vtx->GetNContributors()<3)return kFALSE;
655   // do not want tpc only primary vertex
656   TString vtxName(vtx->GetName());
657   if(vtxName.Contains("TPCVertex"))return kFALSE;
658
659   // no dispersion yet...
660   /* 
661   TString vtxTitle(vtx->GetTitle());
662   if(vtxTitle.Contains("vertexer: Z")){
663     if(vtx->GetDispersion()>0.02)return kFALSE;
664   }
665   */
666   return kTRUE;
667 }
668
669
670 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
671
672   if(!IsVertexValid(vtx))return kFALSE;
673
674   Float_t zvtx = vtx->GetZ();
675   Float_t yvtx = vtx->GetY();
676   Float_t xvtx = vtx->GetX();
677
678   xvtx -= fVtxXMean;
679   yvtx -= fVtxYMean;
680   zvtx -= fVtxZMean;
681
682
683
684   if(TMath::Abs(zvtx)>fVtxZCut){
685     return kFALSE;
686   }
687   fEventCutInfoESD |= kVertexZCut;  
688   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
689   if(r2>(fVtxRCut*fVtxRCut)){
690     return kFALSE;
691   }
692   fEventCutInfoESD |= kVertexRCut;  
693   return kTRUE;
694 }
695
696
697 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
698
699   if(!IsVertexValid(vtx))return kFALSE;
700
701   Float_t zvtx = vtx->GetZ();
702   Float_t yvtx = vtx->GetY();
703   Float_t xvtx = vtx->GetX();
704
705   xvtx -= fVtxXMean;
706   yvtx -= fVtxYMean;
707   zvtx -= fVtxZMean;
708
709   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
710
711   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
712   return vertexIn;
713 }
714
715 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
716   if(!esd)return kFALSE;
717   return esd->IsPileupFromSPD();
718 }
719
720 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
721   if(!esd)return kFALSE;
722   // add track cuts for which we look for cosmics...
723
724   Bool_t isCosmic = kFALSE;
725   Int_t nTracks = esd->GetNumberOfTracks();
726   Int_t nCosmicCandidates = 0;
727
728   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
729     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
730     if (!track1)  continue;
731     UInt_t status1 = track1->GetStatus();
732     //If track is ITS stand alone track, skip the track
733     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
734     if(track1->Pt()<fPtMinCosmic) continue;
735     //Start 2nd track loop to look for correlations
736     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
737       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
738       if(!track2) continue;
739       UInt_t status2 = track2->GetStatus();
740       //If track is ITS stand alone track, skip the track
741       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
742       if(track2->Pt()<fPtMinCosmic) continue;
743       //Check if back-to-back
744       Double_t mom1[3],mom2[3];
745       track1->GetPxPyPz(mom1);
746       track2->GetPxPyPz(mom2);
747       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
748       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
749       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
750       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
751
752       Float_t deltaPhi = track1->Phi()-track2->Phi();
753       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
754
755       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
756       if(rIsol<fRIsolMinCosmic) continue;
757
758       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
759         nCosmicCandidates+=1;
760         isCosmic = kTRUE;
761       }
762       
763     }
764   }
765
766   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
767
768   return isCosmic;
769 }
770
771
772 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
773
774   Float_t cent = 999;
775   if(esd->GetCentrality()){
776     cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
777   }
778   if(cent>50)return 4;
779   if(cent>30)return 3;
780   if(cent>10)return 2;
781   return 1;
782
783 }
784
785
786 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
787
788   Float_t cent = aod->GetHeader()->GetCentrality();
789
790   if(cent>50)return 4;
791   if(cent>30)return 3;
792   if(cent>10)return 2;
793   return 1;
794
795 }
796
797
798 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
799 {
800   // Terminate analysis
801   //
802 }
803