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