coverity fixes
[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"
b6dd6ad2 41#include "AliCentrality.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();
40440b28 194 fHistList->SetOwner();
195 PostData(1, fHistList); // post data in any case once
bf7b8731 196
197 Bool_t oldStatus = TH1::AddDirectoryStatus();
198 TH1::AddDirectory(kFALSE);
bf7b8731 199 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
200 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
201 fHistList->Add(fh1Xsec);
202
203 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
204 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
205 fHistList->Add(fh1Trials);
206
207 const Int_t nBinPt = 100;
208 Double_t binLimitsPt[nBinPt+1];
209 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
210 if(iPt == 0){
211 binLimitsPt[iPt] = 0.0;
212 }
213 else {// 1.0
214 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
215 }
216 }
217
218 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
219 fHistList->Add(fh2TriggerCount);
220
221 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
222 fHistList->Add(fh2ESDTriggerCount);
57ca1193 223 const Int_t nBins = AliAnalysisHelperJetTasks::kTrigger*kConstraints;
224 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
bf7b8731 225 fHistList->Add(fh2TriggerVtx);
226
57ca1193 227 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
bf7b8731 228 fHistList->Add(fh2ESDTriggerVtx);
229
230
231 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
232 fHistList->Add(fh1PtHard);
233 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
234 fHistList->Add(fh1PtHardTrials);
f2dd0695 235 fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
236 AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
f2dd0695 237 fHistList->Add(fh1SelectionInfoESD);
6d597ca2 238
239 fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
240 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
241 fHistList->Add(fh1EventCutInfoESD);
242
b5a3f310 243 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
244 // 3 triggers BB BE/EB EE
245
f4132e7d 246 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 247 fHistList->Add(fh2ESDTriggerRun);
248
249 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
250 fHistList->Add(fh2VtxXY);
bf7b8731 251 // =========== Switch on Sumw2 for all histos ===========
252 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
253 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
254 if (h1){
255 h1->Sumw2();
256 continue;
257 }
258 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
259 if(hn)hn->Sumw2();
260 }
261
df7848fc 262 fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
263 fHistList->Add(fh1NCosmicsPerEvent),
264
bf7b8731 265
266 TH1::AddDirectory(oldStatus);
f4132e7d 267
268 // Add an AOD branch for replication
269 if(fNonStdFile.Length()){
270 if (fDebug > 1) AliInfo("Replicating header");
271 fgAODHeader = new AliAODHeader;
272 AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
070f06d2 273 if (fDebug > 1) AliInfo("Replicating primary vertices");
070f06d2 274 fgAODVertices = new TClonesArray("AliAODVertex",3);
80a790fd 275 fgAODVertices->SetName("vertices");
070f06d2 276 AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
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);
20b170eb 370
f4132e7d 371
6d597ca2 372 fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
373 fh1EventCutInfoESD->Fill(fEventCutInfoESD);
374
f2dd0695 375 if(esdEventSelected) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kVertexIn;
376 if(esdEventPileUp) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsPileUp;
377 if(esdEventCosmic) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsCosmic;
378 if(physicsSelection) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kPhysicsSelection;
379
380
381 // here we have all selection information, fill histogram
382 for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
383 if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
384 }
6d597ca2 385 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
c0997643 386
387 if(esd&&aod&&fDebug){
388 if(esdEventSelected&&!aodEventSelected){
389 Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
390 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
bf7b8731 391 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
c0997643 392 Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
393 vtxESD->Print();
394 Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
395 vtxAOD->Print();
bf7b8731 396 }
c0997643 397 }
398
399 // loop over all possible triggers for esd
400
401 if(esd){
402 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
f4132e7d 403 esdVtxValid = IsVertexValid(vtxESD);
404 esdVtxIn = IsVertexIn(vtxESD);
d8f21f85 405 if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
406 Float_t cent = aod->GetHeader()->GetCentrality();
40440b28 407 if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
d8f21f85 408 if(cent<=80&&esdVtxIn){
409 aodH->SetFillAOD(kTRUE);
c9931fda 410 aodH->SetFillExtension(kTRUE);
d8f21f85 411 }
412 }
413
414
57ca1193 415 Float_t zvtx = vtxESD->GetZ();
f4132e7d 416 Int_t iCl = GetEventClass(esd);
417 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
418 Bool_t cand = physicsSelection;
419
ffab794c 420 if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
d8f21f85 421 fh2ESDTriggerCount->Fill(0.,kAllTriggered);
422 fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
f4132e7d 423 if(cand){
424 fh2ESDTriggerCount->Fill(0.,kSelectedALICE);
425 fh2ESDTriggerCount->Fill(iCl,kSelectedALICE);
426 fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
427 }
428 // if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
429 if(esdVtxValid){
430 fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
431 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
432 fh2ESDTriggerVtx->Fill(iCl,zvtx);
433 if(esdVtxIn){
434 fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
435 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
436 fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
bf7b8731 437 }
f4132e7d 438 if(cand){
439 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
440 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
441 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
c0997643 442 }
bf7b8731 443 }
f4132e7d 444 if(cand&&esdVtxIn){
445 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
446 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
447 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
448 fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
449 fh2ESDTriggerCount->Fill(iCl,kSelected);
450 fh2ESDTriggerCount->Fill(0.,kSelected);
451 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
452 }
bf7b8731 453 }
454
ffab794c 455
456
c0997643 457 if(aod){
458 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
6d597ca2 459 aodVtxValid = IsVertexValid(vtxAOD);
460 aodVtxIn = IsVertexIn(vtxAOD);
ffab794c 461 Float_t zvtx = vtxAOD->GetZ();
462 Int_t iCl = GetEventClass(aod);
463 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
464 Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
465 if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
d8f21f85 466 fh2TriggerCount->Fill(0.,kAllTriggered);
467 fh2TriggerCount->Fill(iCl,kAllTriggered);
ffab794c 468 if(cand){
469 fh2TriggerCount->Fill(0.,kSelectedALICE);
470 fh2TriggerCount->Fill(iCl,kSelectedALICE);
471 fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
472 }
473 if(aodVtxValid){
474 fh2TriggerCount->Fill(0.,kTriggeredVertex);
475 fh2TriggerCount->Fill(iCl,kTriggeredVertex);
476 fh2TriggerVtx->Fill(iCl,zvtx);
477 if(aodVtxIn){
478 fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
479 fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
480 fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
c0997643 481 }
ffab794c 482 if(cand){
483 fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
484 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
485 fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
486 }
487 }
488 if(cand&&aodVtxIn){
489 fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
490 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
491 fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
492 fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
493 fh2TriggerCount->Fill(iCl,kSelected);
494 fh2TriggerCount->Fill(0.,kSelected);
221d8aa9 495 if(fUseAODInput){
496 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
497 if(fFilterAODCollisions&&aod){
498 Float_t cent = aod->GetHeader()->GetCentrality();
499 if(cent<=80){
500 aodH->SetFillAOD(kTRUE);
501 }
502 }
503 }
504
505
506 }
c0997643 507 }
bf7b8731 508
509 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
510
511
512 Double_t ptHard = 0;
513 Double_t nTrials = 1; // Trials for MC trigger
514
515 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
516 AliMCEvent* mcEvent = MCEvent();
517 // AliStack *pStack = 0;
518 if(mcEvent){
519 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
520 if(pythiaGenHeader){
521 nTrials = pythiaGenHeader->Trials();
522 ptHard = pythiaGenHeader->GetPtHard();
523 int iProcessType = pythiaGenHeader->ProcessType();
524 // 11 f+f -> f+f
525 // 12 f+barf -> f+barf
526 // 13 f+barf -> g+g
527 // 28 f+g -> f+g
528 // 53 g+g -> f+barf
529 // 68 g+g -> g+g
530 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
531 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
532 fh1PtHard->Fill(ptHard);
533 fh1PtHardTrials->Fill(ptHard,nTrials);
534
535 }// if pythia gen header
536 }
537
538 // trigger selection
539
f4132e7d 540 // replication of
541 if(fNonStdFile.Length()&&aod){
542 if (fgAODHeader){
543 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
544 }
80a790fd 545 if(fgAODVertices){
80a790fd 546 fgAODVertices->Delete();
547 TClonesArray &vertices = *fgAODVertices;
80a790fd 548 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
070f06d2 549 // we only use some basic information,
550
551
552 Double_t pos[3];
553 Double_t covVtx[6];
554
555 vtxAOD->GetXYZ(pos); // position
556 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
557 Int_t jVertices = 0;
558 AliAODVertex * vtx = new(vertices[jVertices++])
559 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
560 vtx->SetName(vtxAOD->GetName());
561 vtx->SetTitle(vtxAOD->GetTitle());
070f06d2 562 TString vtitle = vtxAOD->GetTitle();
d8f21f85 563 vtx->SetNContributors(vtxAOD->GetNContributors());
070f06d2 564
565 // Add SPD "main" vertex
566 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
567 vtxS->GetXYZ(pos); // position
568 vtxS->GetCovMatrix(covVtx); //covariance matrix
569 AliAODVertex * mVSPD = new(vertices[jVertices++])
570 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
571 mVSPD->SetName(vtxS->GetName());
572 mVSPD->SetTitle(vtxS->GetTitle());
573 mVSPD->SetNContributors(vtxS->GetNContributors());
574
575 // Add tpc only vertex
576 if(esd){
577 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
578 vtxT->GetXYZ(pos); // position
579 vtxT->GetCovMatrix(covVtx); //covariance matrix
580 AliAODVertex * mVTPC = new(vertices[jVertices++])
581 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
582 mVTPC->SetName(vtxT->GetName());
583 mVTPC->SetTitle(vtxT->GetTitle());
584 mVTPC->SetNContributors(vtxT->GetNContributors());
585 }
80a790fd 586 }
f4132e7d 587 }
588
bf7b8731 589 PostData(1, fHistList);
80a790fd 590}
bf7b8731 591
6d597ca2 592Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
c0997643 593 if(!esd)return kFALSE;
594 const AliESDVertex *vtx = esd->GetPrimaryVertex();
6d597ca2 595 return IsVertexIn(vtx); // vertex in calls vertex valid
596}
597
80a790fd 598AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
599 if(fgAODVertices){
600 fgAODVertices->Delete();
601 delete fgAODVertices;
602 }
070f06d2 603 if(fgAODHeader)delete fgAODHeader;
80a790fd 604}
605
6d597ca2 606
607Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
608 if(!aod)return kFALSE;
609 const AliAODVertex *vtx = aod->GetPrimaryVertex();
610 return IsVertexIn(vtx); // VertexIn calls VertexValid
611}
612
613Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
614
c0997643 615 // We can check the type of the vertex by its name and title
616 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
617 // vtx name title
618 // Tracks PrimaryVertex VertexerTracksNoConstraint
6d597ca2 619 // Tracks PrimaryVertex VertexerTracksWithConstraint
c0997643 620 // TPC TPCVertex VertexerTracksNoConstraint
6d597ca2 621 // TPC TPCVertex VertexerTracksWithConstraint
c0997643 622 // SPD SPDVertex vertexer: 3D
623 // SPD SPDVertex vertexer: Z
624
6d597ca2 625 Int_t nCont = vtx->GetNContributors();
6d597ca2 626 if(nCont>=1){
627 fEventCutInfoESD |= kContributorsCut1;
628 if(nCont>=2){
629 fEventCutInfoESD |= kContributorsCut2;
630 if(nCont>=3){
631 fEventCutInfoESD |= kContributorsCut3;
632 }
633 }
634 }
635
636 if(nCont<3)return kFALSE;
c0997643 637
6d597ca2 638 // do not want tpc only primary vertex
639 TString vtxName(vtx->GetName());
640 if(vtxName.Contains("TPCVertex")){
641 fEventCutInfoESD |= kVertexTPC;
642 return kFALSE;
643 }
644 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
645 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
c0997643 646
c0997643 647
6d597ca2 648 TString vtxTitle(vtx->GetTitle());
649 if(vtxTitle.Contains("vertexer: Z")){
650 if(vtx->GetDispersion()>0.02)return kFALSE;
651 }
652 fEventCutInfoESD |= kSPDDispersionCut;
653 return kTRUE;
654}
655
656
657Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
658
659 // We can check the type of the vertex by its name and title
660 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
661 // vtx name title
662 // Tracks PrimaryVertex VertexerTracksNoConstraint
663 // TPC TPCVertex VertexerTracksNoConstraint
664 // SPD SPDVertex vertexer: 3D
665 // SPD SPDVertex vertexer: Z
d8f21f85 666
667 if(fDebug){
668 Printf(" n contrib %d",vtx->GetNContributors());
669 vtx->Print();
670 }
6d597ca2 671
d8f21f85 672 // if(vtx->GetNContributors()<3)return kFALSE;
c0997643 673 // do not want tpc only primary vertex
674 TString vtxName(vtx->GetName());
675 if(vtxName.Contains("TPCVertex"))return kFALSE;
6d597ca2 676
677 // no dispersion yet...
678 /*
679 TString vtxTitle(vtx->GetTitle());
680 if(vtxTitle.Contains("vertexer: Z")){
681 if(vtx->GetDispersion()>0.02)return kFALSE;
682 }
683 */
684 return kTRUE;
685}
686
687
688Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
689
690 if(!IsVertexValid(vtx))return kFALSE;
691
c0997643 692 Float_t zvtx = vtx->GetZ();
693 Float_t yvtx = vtx->GetY();
694 Float_t xvtx = vtx->GetX();
695
6d597ca2 696 xvtx -= fVtxXMean;
697 yvtx -= fVtxYMean;
698 zvtx -= fVtxZMean;
699
700
701
702 if(TMath::Abs(zvtx)>fVtxZCut){
703 return kFALSE;
704 }
705 fEventCutInfoESD |= kVertexZCut;
706 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
707 if(r2>(fVtxRCut*fVtxRCut)){
708 return kFALSE;
709 }
710 fEventCutInfoESD |= kVertexRCut;
711 return kTRUE;
712}
713
714
715Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
716
717 if(!IsVertexValid(vtx))return kFALSE;
718
719 Float_t zvtx = vtx->GetZ();
720 Float_t yvtx = vtx->GetY();
721 Float_t xvtx = vtx->GetX();
c0997643 722
723 xvtx -= fVtxXMean;
724 yvtx -= fVtxYMean;
725 zvtx -= fVtxZMean;
726
727 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
728
6d597ca2 729 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
730 return vertexIn;
c0997643 731}
732
6d597ca2 733Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
f2dd0695 734 if(!esd)return kFALSE;
735 return esd->IsPileupFromSPD();
736}
737
6d597ca2 738Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
f2dd0695 739 if(!esd)return kFALSE;
740 // add track cuts for which we look for cosmics...
df7848fc 741
742 Bool_t isCosmic = kFALSE;
743 Int_t nTracks = esd->GetNumberOfTracks();
744 Int_t nCosmicCandidates = 0;
745
746 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
747 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
748 if (!track1) continue;
75946d5d 749 UInt_t status1 = track1->GetStatus();
750 //If track is ITS stand alone track, skip the track
751 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
df7848fc 752 if(track1->Pt()<fPtMinCosmic) continue;
753 //Start 2nd track loop to look for correlations
754 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
755 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
756 if(!track2) continue;
75946d5d 757 UInt_t status2 = track2->GetStatus();
758 //If track is ITS stand alone track, skip the track
759 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
df7848fc 760 if(track2->Pt()<fPtMinCosmic) continue;
761 //Check if back-to-back
762 Double_t mom1[3],mom2[3];
763 track1->GetPxPyPz(mom1);
764 track2->GetPxPyPz(mom2);
765 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
766 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
767 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
768 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
769
770 Float_t deltaPhi = track1->Phi()-track2->Phi();
771 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
772
773 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
774 if(rIsol<fRIsolMinCosmic) continue;
775
776 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
75946d5d 777 nCosmicCandidates+=1;
df7848fc 778 isCosmic = kTRUE;
779 }
780
781 }
782 }
783
784 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
785
786 return isCosmic;
f2dd0695 787}
788
789
f4132e7d 790Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
791
792 Float_t cent = 999;
793 if(esd->GetCentrality()){
794 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
795 }
d8f21f85 796 if(cent>80)return 5;
f4132e7d 797 if(cent>50)return 4;
798 if(cent>30)return 3;
799 if(cent>10)return 2;
800 return 1;
bf7b8731 801
f4132e7d 802}
f2dd0695 803
ffab794c 804
805Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
806
807 Float_t cent = aod->GetHeader()->GetCentrality();
d8f21f85 808 if(cent>80)return 5;
ffab794c 809 if(cent>50)return 4;
810 if(cent>30)return 3;
811 if(cent>10)return 2;
812 return 1;
813
814}
815
816
bf7b8731 817void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
818{
819 // Terminate analysis
820 //
bf7b8731 821}
f4132e7d 822