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