38331a59104b9afa2f7426f6c5d1cf7088e15fca
[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(aodTrig&&aodEventSelected){
427           fh2TriggerVtx->Fill(it,zvtx);
428           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
429         }
430       }
431     }
432   }
433
434   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
435
436   
437   Double_t ptHard = 0; 
438   Double_t nTrials = 1; // Trials for MC trigger 
439
440   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
441   AliMCEvent* mcEvent = MCEvent();
442   //    AliStack *pStack = 0; 
443   if(mcEvent){
444     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
445     if(pythiaGenHeader){
446       nTrials = pythiaGenHeader->Trials();
447       ptHard  = pythiaGenHeader->GetPtHard();
448       int iProcessType = pythiaGenHeader->ProcessType();
449       // 11 f+f -> f+f
450       // 12 f+barf -> f+barf
451       // 13 f+barf -> g+g
452       // 28 f+g -> f+g
453       // 53 g+g -> f+barf
454       // 68 g+g -> g+g
455       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
456       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
457       fh1PtHard->Fill(ptHard);
458       fh1PtHardTrials->Fill(ptHard,nTrials);
459
460     }// if pythia gen header
461   }
462
463   // trigger selection
464   
465
466   PostData(1, fHistList);
467 }
468
469 Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
470   if(!esd)return kFALSE;
471   const AliESDVertex *vtx = esd->GetPrimaryVertex();
472   // We can check the type of the vertex by its name and title
473   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
474   // vtx       name                  title
475   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
476   // TPC       TPCVertex             VertexerTracksNoConstraint
477   // SPD       SPDVertex             vertexer: 3D
478   // SPD       SPDVertex             vertexer: Z
479   
480
481
482   if(vtx->GetNContributors()<3)return kFALSE;
483
484   // do not want tpc only primary vertex
485   TString vtxName(vtx->GetName());
486   if(vtxName.Contains("TPCVertex"))return kFALSE;
487   Float_t zvtx = vtx->GetZ();
488   Float_t yvtx = vtx->GetY();
489   Float_t xvtx = vtx->GetX();
490
491   // here we may fill the histos for run dependence....
492
493   xvtx -= fVtxXMean;
494   yvtx -= fVtxYMean;
495   zvtx -= fVtxZMean;
496
497   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
498
499   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
500   return eventSel;
501 }
502
503 Bool_t AliAnalysisTaskJetServices::IsEventPileUpESD(AliESDEvent* esd){
504   if(!esd)return kFALSE;
505   return esd->IsPileupFromSPD();
506 }
507
508 Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
509   if(!esd)return kFALSE;
510   // add track cuts for which we look for cosmics...
511
512   Bool_t isCosmic = kFALSE;
513   Int_t nTracks = esd->GetNumberOfTracks();
514   Int_t nCosmicCandidates = 0;
515
516   for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
517     AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
518     if (!track1)  continue;
519     UInt_t status1 = track1->GetStatus();
520     //If track is ITS stand alone track, skip the track
521     if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
522     if(track1->Pt()<fPtMinCosmic) continue;
523     //Start 2nd track loop to look for correlations
524     for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
525       AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
526       if(!track2) continue;
527       UInt_t status2 = track2->GetStatus();
528       //If track is ITS stand alone track, skip the track
529       if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
530       if(track2->Pt()<fPtMinCosmic) continue;
531       //Check if back-to-back
532       Double_t mom1[3],mom2[3];
533       track1->GetPxPyPz(mom1);
534       track2->GetPxPyPz(mom2);
535       TVector3 momv1(mom1[0],mom1[1],mom1[2]);
536       TVector3 momv2(mom2[0],mom2[1],mom2[2]);
537       Float_t theta = (float)(momv1.Phi()-momv2.Phi());
538       if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
539
540       Float_t deltaPhi = track1->Phi()-track2->Phi();
541       if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
542
543       Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
544       if(rIsol<fRIsolMinCosmic) continue;
545
546       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
547         nCosmicCandidates+=1;
548         isCosmic = kTRUE;
549       }
550       
551     }
552   }
553
554   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
555
556   return isCosmic;
557 }
558
559
560 Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
561   if(!aod)return kFALSE;
562   const AliAODVertex *vtx = aod->GetPrimaryVertex();
563   // We can check the type of the vertex by its name and title
564   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
565   // vtx       name                  title
566   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
567   // TPC       TPCVertex             VertexerTracksNoConstraint
568   // SPD       SPDVertex             vertexer: 3D
569   // SPD       SPDVertex             vertexer: Z
570   
571
572
573   if(vtx->GetNContributors()<3)return kFALSE;
574
575   // do not want tpc only primary vertex
576   TString vtxName(vtx->GetName());
577   if(vtxName.Contains("TPCVertex"))return kFALSE;
578
579   Float_t zvtx = vtx->GetZ();
580   Float_t yvtx = vtx->GetY();
581   Float_t xvtx = vtx->GetX();
582
583   // here we may fill the histos for run dependence....
584
585   xvtx -= fVtxXMean;
586   yvtx -= fVtxYMean;
587   zvtx -= fVtxZMean;
588
589   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
590   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
591   return eventSel;
592 }
593
594
595
596
597 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
598 {
599   // Terminate analysis
600   //
601 }