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