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