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