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