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