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