Fixes for centrality < 0, adding jet finder setting with higher pT cut
[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),
742ee86c 97 fh1CentralityESD(0),
98 fh1Centrality(0),
bf7b8731 99 fh2TriggerCount(0x0),
100 fh2ESDTriggerCount(0x0),
101 fh2TriggerVtx(0x0),
102 fh2ESDTriggerVtx(0x0),
b5a3f310 103 fh2ESDTriggerRun(0x0),
104 fh2VtxXY(0x0),
df7848fc 105 fh1NCosmicsPerEvent(0x0),
3493c3a9 106 fTriggerAnalysis(0x0),
bf7b8731 107 fHistList(0x0)
108{
b5a3f310 109 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 110}
111
112AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
113 AliAnalysisTaskSE(name),
114 fUseAODInput(kFALSE),
cc0649e4 115 fUsePhysicsSelection(kFALSE),
3493c3a9 116 fMC(kFALSE),
f4132e7d 117 fFilterAODCollisions(kFALSE),
39e7e8ab 118 fPhysicsSelectionFlag(AliVEvent::kMB),
f2dd0695 119 fSelectionInfoESD(0),
6d597ca2 120 fEventCutInfoESD(0),
bf7b8731 121 fAvgTrials(1),
c0997643 122 fVtxXMean(0),
123 fVtxYMean(0),
124 fVtxZMean(0),
125 fVtxRCut(1.),
126 fVtxZCut(8.),
df7848fc 127 fPtMinCosmic(5.),
128 fRIsolMinCosmic(3.),
129 fMaxCosmicAngle(0.01),
f4132e7d 130 fNonStdFile(""),
bf7b8731 131 fh1Xsec(0x0),
132 fh1Trials(0x0),
133 fh1PtHard(0x0),
134 fh1PtHardTrials(0x0),
f2dd0695 135 fh1SelectionInfoESD(0x0),
6d597ca2 136 fh1EventCutInfoESD(0),
742ee86c 137 fh1CentralityESD(0),
138 fh1Centrality(0),
bf7b8731 139 fh2TriggerCount(0x0),
140 fh2ESDTriggerCount(0x0),
141 fh2TriggerVtx(0x0),
142 fh2ESDTriggerVtx(0x0),
b5a3f310 143 fh2ESDTriggerRun(0x0),
144 fh2VtxXY(0x0),
df7848fc 145 fh1NCosmicsPerEvent(0x0),
3493c3a9 146 fTriggerAnalysis(0x0),
bf7b8731 147 fHistList(0x0)
148{
b5a3f310 149 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 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
bf7b8731 194 if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
195
196 OpenFile(1);
197 if(!fHistList)fHistList = new TList();
40440b28 198 fHistList->SetOwner();
199 PostData(1, fHistList); // post data in any case once
bf7b8731 200
201 Bool_t oldStatus = TH1::AddDirectoryStatus();
202 TH1::AddDirectory(kFALSE);
bf7b8731 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
742ee86c 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
bf7b8731 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);
57ca1193 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.);
bf7b8731 236 fHistList->Add(fh2TriggerVtx);
237
57ca1193 238 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
bf7b8731 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);
f2dd0695 246 fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
247 AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
f2dd0695 248 fHistList->Add(fh1SelectionInfoESD);
6d597ca2 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
b5a3f310 254 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
255 // 3 triggers BB BE/EB EE
256
f4132e7d 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);
b5a3f310 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);
bf7b8731 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
df7848fc 273 fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
274 fHistList->Add(fh1NCosmicsPerEvent),
275
bf7b8731 276
277 TH1::AddDirectory(oldStatus);
f4132e7d 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());
070f06d2 284 if (fDebug > 1) AliInfo("Replicating primary vertices");
070f06d2 285 fgAODVertices = new TClonesArray("AliAODVertex",3);
80a790fd 286 fgAODVertices->SetName("vertices");
070f06d2 287 AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
f4132e7d 288 }
bf7b8731 289}
290
291void AliAnalysisTaskJetServices::Init()
292{
293 //
294 // Initialization
295 //
bf7b8731 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;
f4132e7d 309
310
3493c3a9 311
bf7b8731 312 AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
f4132e7d 313 AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
f2dd0695 314 fSelectionInfoESD = 0; // reset
6d597ca2 315 fEventCutInfoESD = 0; // reset
f2dd0695 316 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
317
bf7b8731 318
f4132e7d 319 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
320
321
322
bf7b8731 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 }
ffab794c 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 }
3493c3a9 346 fSelectionInfoESD |= kNoEventCut;
347 fEventCutInfoESD |= kNoEventCut;
6d597ca2 348
349 Bool_t esdVtxValid = false;
350 Bool_t esdVtxIn = false;
351 Bool_t aodVtxValid = false;
352 Bool_t aodVtxIn = false;
f2dd0695 353
b5a3f310 354 if(esd){
3493c3a9 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;
b5a3f310 371 }
f4132e7d 372
c0997643 373 // Apply additional constraints
6d597ca2 374 Bool_t esdEventSelected = IsEventSelected(esd);
375 Bool_t esdEventPileUp = IsEventPileUp(esd);
376 Bool_t esdEventCosmic = IsEventCosmic(esd);
f2dd0695 377
6d597ca2 378 Bool_t aodEventSelected = IsEventSelected(aod);
f2dd0695 379
39e7e8ab 380 Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
20b170eb 381
f4132e7d 382
6d597ca2 383 fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
384 fh1EventCutInfoESD->Fill(fEventCutInfoESD);
385
f2dd0695 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 }
6d597ca2 396 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
c0997643 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();
bf7b8731 402 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
c0997643 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();
bf7b8731 407 }
c0997643 408 }
409
410 // loop over all possible triggers for esd
411
742ee86c 412 Float_t cent = 100;
413 if(aod)cent = aod->GetHeader()->GetCentrality();
414 if(cent<0)cent = 101;
c0997643 415 if(esd){
416 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
f4132e7d 417 esdVtxValid = IsVertexValid(vtxESD);
418 esdVtxIn = IsVertexIn(vtxESD);
d8f21f85 419 if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
40440b28 420 if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
d8f21f85 421 if(cent<=80&&esdVtxIn){
422 aodH->SetFillAOD(kTRUE);
c9931fda 423 aodH->SetFillExtension(kTRUE);
d8f21f85 424 }
425 }
426
427
57ca1193 428 Float_t zvtx = vtxESD->GetZ();
f4132e7d 429 Int_t iCl = GetEventClass(esd);
430 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
431 Bool_t cand = physicsSelection;
432
ffab794c 433 if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
d8f21f85 434 fh2ESDTriggerCount->Fill(0.,kAllTriggered);
435 fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
f4132e7d 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);
bf7b8731 450 }
f4132e7d 451 if(cand){
452 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
453 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
454 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
c0997643 455 }
bf7b8731 456 }
742ee86c 457
c08c3ad2 458 if(cand&&esdVtxIn&&iCl<5){
f4132e7d 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
742ee86c 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 }
f4132e7d 472 }
bf7b8731 473 }
474
ffab794c 475
476
c0997643 477 if(aod){
478 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
6d597ca2 479 aodVtxValid = IsVertexValid(vtxAOD);
480 aodVtxIn = IsVertexIn(vtxAOD);
ffab794c 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());
d8f21f85 486 fh2TriggerCount->Fill(0.,kAllTriggered);
487 fh2TriggerCount->Fill(iCl,kAllTriggered);
ffab794c 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);
c0997643 501 }
ffab794c 502 if(cand){
503 fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
504 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
505 fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
506 }
507 }
c08c3ad2 508 if(cand&&aodVtxIn&&iCl<5){
ffab794c 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);
742ee86c 515 fh1Centrality->Fill(cent);
516 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
517 if(fUseAODInput&&cent<=80){
221d8aa9 518 if(fFilterAODCollisions&&aod){
742ee86c 519 aodH->SetFillAOD(kTRUE);
221d8aa9 520 }
521 }
221d8aa9 522 }
c0997643 523 }
bf7b8731 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
f4132e7d 556 // replication of
557 if(fNonStdFile.Length()&&aod){
558 if (fgAODHeader){
559 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
560 }
80a790fd 561 if(fgAODVertices){
80a790fd 562 fgAODVertices->Delete();
563 TClonesArray &vertices = *fgAODVertices;
80a790fd 564 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
070f06d2 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());
070f06d2 578 TString vtitle = vtxAOD->GetTitle();
d8f21f85 579 vtx->SetNContributors(vtxAOD->GetNContributors());
070f06d2 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 }
80a790fd 602 }
f4132e7d 603 }
604
bf7b8731 605 PostData(1, fHistList);
80a790fd 606}
bf7b8731 607
6d597ca2 608Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
c0997643 609 if(!esd)return kFALSE;
610 const AliESDVertex *vtx = esd->GetPrimaryVertex();
6d597ca2 611 return IsVertexIn(vtx); // vertex in calls vertex valid
612}
613
80a790fd 614AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
615 if(fgAODVertices){
616 fgAODVertices->Delete();
617 delete fgAODVertices;
618 }
070f06d2 619 if(fgAODHeader)delete fgAODHeader;
80a790fd 620}
621
6d597ca2 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
c0997643 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
6d597ca2 635 // Tracks PrimaryVertex VertexerTracksWithConstraint
c0997643 636 // TPC TPCVertex VertexerTracksNoConstraint
6d597ca2 637 // TPC TPCVertex VertexerTracksWithConstraint
c0997643 638 // SPD SPDVertex vertexer: 3D
639 // SPD SPDVertex vertexer: Z
640
6d597ca2 641 Int_t nCont = vtx->GetNContributors();
6d597ca2 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;
c0997643 653
6d597ca2 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;
c0997643 662
c0997643 663
6d597ca2 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
d8f21f85 682
683 if(fDebug){
684 Printf(" n contrib %d",vtx->GetNContributors());
685 vtx->Print();
686 }
6d597ca2 687
d8f21f85 688 // if(vtx->GetNContributors()<3)return kFALSE;
c0997643 689 // do not want tpc only primary vertex
690 TString vtxName(vtx->GetName());
691 if(vtxName.Contains("TPCVertex"))return kFALSE;
6d597ca2 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
c0997643 708 Float_t zvtx = vtx->GetZ();
709 Float_t yvtx = vtx->GetY();
710 Float_t xvtx = vtx->GetX();
711
6d597ca2 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();
c0997643 738
739 xvtx -= fVtxXMean;
740 yvtx -= fVtxYMean;
741 zvtx -= fVtxZMean;
742
743 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
744
6d597ca2 745 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
746 return vertexIn;
c0997643 747}
748
6d597ca2 749Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
f2dd0695 750 if(!esd)return kFALSE;
751 return esd->IsPileupFromSPD();
752}
753
6d597ca2 754Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
f2dd0695 755 if(!esd)return kFALSE;
756 // add track cuts for which we look for cosmics...
df7848fc 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;
75946d5d 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;
df7848fc 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;
75946d5d 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;
df7848fc 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) {
75946d5d 793 nCosmicCandidates+=1;
df7848fc 794 isCosmic = kTRUE;
795 }
796
797 }
798 }
799
800 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
801
802 return isCosmic;
f2dd0695 803}
804
805
f4132e7d 806Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
807
808 Float_t cent = 999;
809 if(esd->GetCentrality()){
810 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
811 }
742ee86c 812 if(cent>80||cent<0)return 5;
f4132e7d 813 if(cent>50)return 4;
814 if(cent>30)return 3;
815 if(cent>10)return 2;
816 return 1;
bf7b8731 817
f4132e7d 818}
f2dd0695 819
ffab794c 820
821Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
822
823 Float_t cent = aod->GetHeader()->GetCentrality();
742ee86c 824 if(cent>80||cent<0)return 5;
ffab794c 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
bf7b8731 833void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
834{
835 // Terminate analysis
836 //
bf7b8731 837}
f4132e7d 838