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