fix to avoid crash when runnign with option 1
[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();
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());
070f06d2 271 if (fDebug > 1) AliInfo("Replicating primary vertices");
272
273 fgAODVertices = new TClonesArray("AliAODVertex",3);
80a790fd 274 fgAODVertices->SetName("vertices");
070f06d2 275 AddAODBranch("TClonesArray",&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);
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();
407 if(cent<=80&&esdVtxIn){
408 aodH->SetFillAOD(kTRUE);
409 }
410 }
411
412
57ca1193 413 Float_t zvtx = vtxESD->GetZ();
f4132e7d 414 Int_t iCl = GetEventClass(esd);
415 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
416 Bool_t cand = physicsSelection;
417
ffab794c 418 if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
d8f21f85 419 fh2ESDTriggerCount->Fill(0.,kAllTriggered);
420 fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
f4132e7d 421 if(cand){
422 fh2ESDTriggerCount->Fill(0.,kSelectedALICE);
423 fh2ESDTriggerCount->Fill(iCl,kSelectedALICE);
424 fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
425 }
426 // if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
427 if(esdVtxValid){
428 fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
429 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
430 fh2ESDTriggerVtx->Fill(iCl,zvtx);
431 if(esdVtxIn){
432 fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
433 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
434 fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
bf7b8731 435 }
f4132e7d 436 if(cand){
437 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
438 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
439 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
c0997643 440 }
bf7b8731 441 }
f4132e7d 442 if(cand&&esdVtxIn){
443 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
444 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
445 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
446 fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
447 fh2ESDTriggerCount->Fill(iCl,kSelected);
448 fh2ESDTriggerCount->Fill(0.,kSelected);
449 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
450 }
bf7b8731 451 }
452
ffab794c 453
454
c0997643 455 if(aod){
456 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
6d597ca2 457 aodVtxValid = IsVertexValid(vtxAOD);
458 aodVtxIn = IsVertexIn(vtxAOD);
ffab794c 459 Float_t zvtx = vtxAOD->GetZ();
460 Int_t iCl = GetEventClass(aod);
461 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
462 Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
463 if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
d8f21f85 464 fh2TriggerCount->Fill(0.,kAllTriggered);
465 fh2TriggerCount->Fill(iCl,kAllTriggered);
ffab794c 466 if(cand){
467 fh2TriggerCount->Fill(0.,kSelectedALICE);
468 fh2TriggerCount->Fill(iCl,kSelectedALICE);
469 fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
470 }
471 if(aodVtxValid){
472 fh2TriggerCount->Fill(0.,kTriggeredVertex);
473 fh2TriggerCount->Fill(iCl,kTriggeredVertex);
474 fh2TriggerVtx->Fill(iCl,zvtx);
475 if(aodVtxIn){
476 fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
477 fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
478 fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
c0997643 479 }
ffab794c 480 if(cand){
481 fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
482 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
483 fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
484 }
485 }
486 if(cand&&aodVtxIn){
487 fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
488 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
489 fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
490 fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
491 fh2TriggerCount->Fill(iCl,kSelected);
492 fh2TriggerCount->Fill(0.,kSelected);
221d8aa9 493 if(fUseAODInput){
494 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
495 if(fFilterAODCollisions&&aod){
496 Float_t cent = aod->GetHeader()->GetCentrality();
497 if(cent<=80){
498 aodH->SetFillAOD(kTRUE);
499 }
500 }
501 }
502
503
504 }
c0997643 505 }
bf7b8731 506
507 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
508
509
510 Double_t ptHard = 0;
511 Double_t nTrials = 1; // Trials for MC trigger
512
513 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
514 AliMCEvent* mcEvent = MCEvent();
515 // AliStack *pStack = 0;
516 if(mcEvent){
517 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
518 if(pythiaGenHeader){
519 nTrials = pythiaGenHeader->Trials();
520 ptHard = pythiaGenHeader->GetPtHard();
521 int iProcessType = pythiaGenHeader->ProcessType();
522 // 11 f+f -> f+f
523 // 12 f+barf -> f+barf
524 // 13 f+barf -> g+g
525 // 28 f+g -> f+g
526 // 53 g+g -> f+barf
527 // 68 g+g -> g+g
528 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
529 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
530 fh1PtHard->Fill(ptHard);
531 fh1PtHardTrials->Fill(ptHard,nTrials);
532
533 }// if pythia gen header
534 }
535
536 // trigger selection
537
f4132e7d 538 // replication of
539 if(fNonStdFile.Length()&&aod){
540 if (fgAODHeader){
541 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
542 }
80a790fd 543 if(fgAODVertices){
80a790fd 544 fgAODVertices->Delete();
545 TClonesArray &vertices = *fgAODVertices;
80a790fd 546 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
070f06d2 547 // we only use some basic information,
548
549
550 Double_t pos[3];
551 Double_t covVtx[6];
552
553 vtxAOD->GetXYZ(pos); // position
554 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
555 Int_t jVertices = 0;
556 AliAODVertex * vtx = new(vertices[jVertices++])
557 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
558 vtx->SetName(vtxAOD->GetName());
559 vtx->SetTitle(vtxAOD->GetTitle());
560
561 TString vtitle = vtxAOD->GetTitle();
d8f21f85 562 vtx->SetNContributors(vtxAOD->GetNContributors());
070f06d2 563
564 // Add SPD "main" vertex
565 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
566 vtxS->GetXYZ(pos); // position
567 vtxS->GetCovMatrix(covVtx); //covariance matrix
568 AliAODVertex * mVSPD = new(vertices[jVertices++])
569 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
570 mVSPD->SetName(vtxS->GetName());
571 mVSPD->SetTitle(vtxS->GetTitle());
572 mVSPD->SetNContributors(vtxS->GetNContributors());
573
574 // Add tpc only vertex
575 if(esd){
576 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
577 vtxT->GetXYZ(pos); // position
578 vtxT->GetCovMatrix(covVtx); //covariance matrix
579 AliAODVertex * mVTPC = new(vertices[jVertices++])
580 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
581 mVTPC->SetName(vtxT->GetName());
582 mVTPC->SetTitle(vtxT->GetTitle());
583 mVTPC->SetNContributors(vtxT->GetNContributors());
584 }
80a790fd 585 }
f4132e7d 586 }
587
bf7b8731 588 PostData(1, fHistList);
80a790fd 589}
bf7b8731 590
6d597ca2 591Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
c0997643 592 if(!esd)return kFALSE;
593 const AliESDVertex *vtx = esd->GetPrimaryVertex();
6d597ca2 594 return IsVertexIn(vtx); // vertex in calls vertex valid
595}
596
80a790fd 597AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
598 if(fgAODVertices){
599 fgAODVertices->Delete();
600 delete fgAODVertices;
601 }
070f06d2 602 if(fgAODHeader)delete fgAODHeader;
80a790fd 603}
604
6d597ca2 605
606Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
607 if(!aod)return kFALSE;
608 const AliAODVertex *vtx = aod->GetPrimaryVertex();
609 return IsVertexIn(vtx); // VertexIn calls VertexValid
610}
611
612Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
613
c0997643 614 // We can check the type of the vertex by its name and title
615 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
616 // vtx name title
617 // Tracks PrimaryVertex VertexerTracksNoConstraint
6d597ca2 618 // Tracks PrimaryVertex VertexerTracksWithConstraint
c0997643 619 // TPC TPCVertex VertexerTracksNoConstraint
6d597ca2 620 // TPC TPCVertex VertexerTracksWithConstraint
c0997643 621 // SPD SPDVertex vertexer: 3D
622 // SPD SPDVertex vertexer: Z
623
6d597ca2 624 Int_t nCont = vtx->GetNContributors();
6d597ca2 625 if(nCont>=1){
626 fEventCutInfoESD |= kContributorsCut1;
627 if(nCont>=2){
628 fEventCutInfoESD |= kContributorsCut2;
629 if(nCont>=3){
630 fEventCutInfoESD |= kContributorsCut3;
631 }
632 }
633 }
634
635 if(nCont<3)return kFALSE;
c0997643 636
6d597ca2 637 // do not want tpc only primary vertex
638 TString vtxName(vtx->GetName());
639 if(vtxName.Contains("TPCVertex")){
640 fEventCutInfoESD |= kVertexTPC;
641 return kFALSE;
642 }
643 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
644 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
c0997643 645
c0997643 646
6d597ca2 647 TString vtxTitle(vtx->GetTitle());
648 if(vtxTitle.Contains("vertexer: Z")){
649 if(vtx->GetDispersion()>0.02)return kFALSE;
650 }
651 fEventCutInfoESD |= kSPDDispersionCut;
652 return kTRUE;
653}
654
655
656Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
657
658 // We can check the type of the vertex by its name and title
659 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
660 // vtx name title
661 // Tracks PrimaryVertex VertexerTracksNoConstraint
662 // TPC TPCVertex VertexerTracksNoConstraint
663 // SPD SPDVertex vertexer: 3D
664 // SPD SPDVertex vertexer: Z
d8f21f85 665
666 if(fDebug){
667 Printf(" n contrib %d",vtx->GetNContributors());
668 vtx->Print();
669 }
6d597ca2 670
d8f21f85 671 // if(vtx->GetNContributors()<3)return kFALSE;
c0997643 672 // do not want tpc only primary vertex
673 TString vtxName(vtx->GetName());
674 if(vtxName.Contains("TPCVertex"))return kFALSE;
6d597ca2 675
676 // no dispersion yet...
677 /*
678 TString vtxTitle(vtx->GetTitle());
679 if(vtxTitle.Contains("vertexer: Z")){
680 if(vtx->GetDispersion()>0.02)return kFALSE;
681 }
682 */
683 return kTRUE;
684}
685
686
687Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
688
689 if(!IsVertexValid(vtx))return kFALSE;
690
c0997643 691 Float_t zvtx = vtx->GetZ();
692 Float_t yvtx = vtx->GetY();
693 Float_t xvtx = vtx->GetX();
694
6d597ca2 695 xvtx -= fVtxXMean;
696 yvtx -= fVtxYMean;
697 zvtx -= fVtxZMean;
698
699
700
701 if(TMath::Abs(zvtx)>fVtxZCut){
702 return kFALSE;
703 }
704 fEventCutInfoESD |= kVertexZCut;
705 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
706 if(r2>(fVtxRCut*fVtxRCut)){
707 return kFALSE;
708 }
709 fEventCutInfoESD |= kVertexRCut;
710 return kTRUE;
711}
712
713
714Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
715
716 if(!IsVertexValid(vtx))return kFALSE;
717
718 Float_t zvtx = vtx->GetZ();
719 Float_t yvtx = vtx->GetY();
720 Float_t xvtx = vtx->GetX();
c0997643 721
722 xvtx -= fVtxXMean;
723 yvtx -= fVtxYMean;
724 zvtx -= fVtxZMean;
725
726 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
727
6d597ca2 728 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
729 return vertexIn;
c0997643 730}
731
6d597ca2 732Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
f2dd0695 733 if(!esd)return kFALSE;
734 return esd->IsPileupFromSPD();
735}
736
6d597ca2 737Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
f2dd0695 738 if(!esd)return kFALSE;
739 // add track cuts for which we look for cosmics...
df7848fc 740
741 Bool_t isCosmic = kFALSE;
742 Int_t nTracks = esd->GetNumberOfTracks();
743 Int_t nCosmicCandidates = 0;
744
745 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
746 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
747 if (!track1) continue;
75946d5d 748 UInt_t status1 = track1->GetStatus();
749 //If track is ITS stand alone track, skip the track
750 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
df7848fc 751 if(track1->Pt()<fPtMinCosmic) continue;
752 //Start 2nd track loop to look for correlations
753 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
754 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
755 if(!track2) continue;
75946d5d 756 UInt_t status2 = track2->GetStatus();
757 //If track is ITS stand alone track, skip the track
758 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
df7848fc 759 if(track2->Pt()<fPtMinCosmic) continue;
760 //Check if back-to-back
761 Double_t mom1[3],mom2[3];
762 track1->GetPxPyPz(mom1);
763 track2->GetPxPyPz(mom2);
764 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
765 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
766 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
767 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
768
769 Float_t deltaPhi = track1->Phi()-track2->Phi();
770 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
771
772 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
773 if(rIsol<fRIsolMinCosmic) continue;
774
775 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
75946d5d 776 nCosmicCandidates+=1;
df7848fc 777 isCosmic = kTRUE;
778 }
779
780 }
781 }
782
783 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
784
785 return isCosmic;
f2dd0695 786}
787
788
f4132e7d 789Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
790
791 Float_t cent = 999;
792 if(esd->GetCentrality()){
793 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
794 }
d8f21f85 795 if(cent>80)return 5;
f4132e7d 796 if(cent>50)return 4;
797 if(cent>30)return 3;
798 if(cent>10)return 2;
799 return 1;
bf7b8731 800
f4132e7d 801}
f2dd0695 802
ffab794c 803
804Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
805
806 Float_t cent = aod->GetHeader()->GetCentrality();
d8f21f85 807 if(cent>80)return 5;
ffab794c 808 if(cent>50)return 4;
809 if(cent>30)return 3;
810 if(cent>10)return 2;
811 return 1;
812
813}
814
815
bf7b8731 816void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
817{
818 // Terminate analysis
819 //
bf7b8731 820}
f4132e7d 821