Class updated to read the threshold file generated online and propagated to offline...
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetServices.cxx
CommitLineData
d3a3f33d 1
bf7b8731 2// **************************************
3// Task used for the correction of determiantion of reconstructed jet spectra
4// Compares input (gen) and output (rec) jets
5// *******************************************
6
7
8/**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
13 * *
14 * Permission to use, copy, modify and distribute this software and its *
15 * documentation strictly for non-commercial purposes is hereby granted *
16 * without fee, provided that the above copyright notice appears in all *
17 * copies and that both the copyright notice and this permission notice *
18 * appear in the supporting documentation. The authors make no claims *
19 * about the suitability of this software for any purpose. It is *
20 * provided "as is" without express or implied warranty. *
21 **************************************************************************/
22
23
24#include <TROOT.h>
25#include <TRandom.h>
b1cd90b5 26#include <TString.h>
bf7b8731 27#include <TSystem.h>
28#include <TInterpreter.h>
29#include <TChain.h>
30#include <TFile.h>
31#include <TKey.h>
32#include <TH1F.h>
33#include <TH2F.h>
34#include <TH3F.h>
35#include <TProfile.h>
36#include <TList.h>
37#include <TLorentzVector.h>
40445651 38#include <TRandom3.h>
bf7b8731 39#include <TClonesArray.h>
c0997643 40#include "TDatabasePDG.h"
bf7b8731 41
42#include "AliAnalysisTaskJetServices.h"
b6dd6ad2 43#include "AliCentrality.h"
fe669ac6 44#include "AliAnalysisDataContainer.h"
45#include "AliAnalysisDataSlot.h"
bf7b8731 46#include "AliAnalysisManager.h"
47#include "AliJetFinder.h"
48#include "AliJetHeader.h"
49#include "AliJetReader.h"
50#include "AliJetReaderHeader.h"
51#include "AliUA1JetHeaderV1.h"
52#include "AliESDEvent.h"
53#include "AliAODEvent.h"
54#include "AliAODHandler.h"
55#include "AliAODTrack.h"
2eded560 56#include "AliAODVZERO.h"
bf7b8731 57#include "AliAODJet.h"
58#include "AliAODMCParticle.h"
59#include "AliMCEventHandler.h"
60#include "AliMCEvent.h"
61#include "AliStack.h"
62#include "AliGenPythiaEventHeader.h"
63#include "AliJetKineReaderHeader.h"
64#include "AliGenCocktailEventHeader.h"
65#include "AliInputEventHandler.h"
fe669ac6 66#include "AliPhysicsSelection.h"
3493c3a9 67#include "AliTriggerAnalysis.h"
bf7b8731 68
69#include "AliAnalysisHelperJetTasks.h"
70
71ClassImp(AliAnalysisTaskJetServices)
72
2eded560 73AliAODVZERO* AliAnalysisTaskJetServices::fgAODVZERO = NULL;
f4132e7d 74AliAODHeader* AliAnalysisTaskJetServices::fgAODHeader = NULL;
80a790fd 75TClonesArray* AliAnalysisTaskJetServices::fgAODVertices = NULL;
f4132e7d 76
40445651 77AliAnalysisTaskJetServices::AliAnalysisTaskJetServices():
78 AliAnalysisTaskSE(),
bf7b8731 79 fUseAODInput(kFALSE),
cc0649e4 80 fUsePhysicsSelection(kFALSE),
3493c3a9 81 fMC(kFALSE),
f4132e7d 82 fFilterAODCollisions(kFALSE),
39e7e8ab 83 fPhysicsSelectionFlag(AliVEvent::kMB),
f2dd0695 84 fSelectionInfoESD(0),
6d597ca2 85 fEventCutInfoESD(0),
40445651 86 fFilterMask(0),
2eded560 87 fRPMethod(0),
21df3cc7 88 fCollisionType(kPbPb),
3d1dba29 89 fNTrigger(0),
bf7b8731 90 fAvgTrials(1),
c0997643 91 fVtxXMean(0),
92 fVtxYMean(0),
93 fVtxZMean(0),
94 fVtxRCut(1.),
95 fVtxZCut(8.),
df7848fc 96 fPtMinCosmic(5.),
97 fRIsolMinCosmic(3.),
98 fMaxCosmicAngle(0.01),
40445651 99 fCentrality(101),
100 fTrackRecEtaWindow(0.9),
101 fMinTrackPt(0.15),
102 fRPAngle(0),
2eded560 103 fPsiVZEROA(0),
104 fPsiVZEROC(0),
3d1dba29 105 fTriggerBit(0x0),
40445651 106 fRandomizer(0),
f4132e7d 107 fNonStdFile(""),
3d1dba29 108 fTriggerName(0x0),
bf7b8731 109 fh1Xsec(0x0),
110 fh1Trials(0x0),
432abb2b 111 fh1AvgTrials(0x0),
bf7b8731 112 fh1PtHard(0x0),
113 fh1PtHardTrials(0x0),
f2dd0695 114 fh1SelectionInfoESD(0x0),
6d597ca2 115 fh1EventCutInfoESD(0),
742ee86c 116 fh1CentralityESD(0),
117 fh1Centrality(0),
40445651 118 fh1RP(0),
e5697f4d 119 fh2ReducedTrigger(0),
3d1dba29 120 fh2CentralityTriggerESD(0),
121 fh2CentralityTrigger(0),
bf7b8731 122 fh2TriggerCount(0x0),
123 fh2ESDTriggerCount(0x0),
124 fh2TriggerVtx(0x0),
125 fh2ESDTriggerVtx(0x0),
b5a3f310 126 fh2ESDTriggerRun(0x0),
127 fh2VtxXY(0x0),
df7848fc 128 fh1NCosmicsPerEvent(0x0),
2eded560 129 fp1RPXA(0x0),
130 fp1RPYA(0x0),
131 fp1RPXC(0x0),
132 fp1RPYC(0x0),
b360e4be 133 fp1CalibRPXA(0x0),
134 fp1CalibRPYA(0x0),
135 fp1CalibRPXC(0x0),
136 fp1CalibRPYC(0x0),
2eded560 137 fh2RPAC(0x0),
138 fh2RPAT(0x0),
139 fh2RPCT(0x0),
140 fh2XYA(0x0),
141 fh2XYC(0x0),
40445651 142 fh2RPCentrality(0x0),
2eded560 143 fh2RPACentrality(0x0),
144 fh2RPCCentrality(0x0),
3493c3a9 145 fTriggerAnalysis(0x0),
bf7b8731 146 fHistList(0x0)
147{
b5a3f310 148 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 149}
150
151AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
152 AliAnalysisTaskSE(name),
153 fUseAODInput(kFALSE),
cc0649e4 154 fUsePhysicsSelection(kFALSE),
3493c3a9 155 fMC(kFALSE),
f4132e7d 156 fFilterAODCollisions(kFALSE),
39e7e8ab 157 fPhysicsSelectionFlag(AliVEvent::kMB),
f2dd0695 158 fSelectionInfoESD(0),
6d597ca2 159 fEventCutInfoESD(0),
40445651 160 fFilterMask(0),
2eded560 161 fRPMethod(0),
21df3cc7 162 fCollisionType(kPbPb),
3d1dba29 163 fNTrigger(0),
bf7b8731 164 fAvgTrials(1),
c0997643 165 fVtxXMean(0),
166 fVtxYMean(0),
167 fVtxZMean(0),
168 fVtxRCut(1.),
169 fVtxZCut(8.),
df7848fc 170 fPtMinCosmic(5.),
171 fRIsolMinCosmic(3.),
172 fMaxCosmicAngle(0.01),
40445651 173 fCentrality(101),
174 fTrackRecEtaWindow(0.9),
175 fMinTrackPt(0.15),
176 fRPAngle(0),
2eded560 177 fPsiVZEROA(0),
178 fPsiVZEROC(0),
3d1dba29 179 fTriggerBit(0x0),
40445651 180 fRandomizer(0),
f4132e7d 181 fNonStdFile(""),
3d1dba29 182 fTriggerName(0x0),
bf7b8731 183 fh1Xsec(0x0),
184 fh1Trials(0x0),
432abb2b 185 fh1AvgTrials(0x0),
bf7b8731 186 fh1PtHard(0x0),
187 fh1PtHardTrials(0x0),
f2dd0695 188 fh1SelectionInfoESD(0x0),
6d597ca2 189 fh1EventCutInfoESD(0),
742ee86c 190 fh1CentralityESD(0),
191 fh1Centrality(0),
40445651 192 fh1RP(0),
e5697f4d 193 fh2ReducedTrigger(0),
3d1dba29 194 fh2CentralityTriggerESD(0),
195 fh2CentralityTrigger(0),
bf7b8731 196 fh2TriggerCount(0x0),
197 fh2ESDTriggerCount(0x0),
198 fh2TriggerVtx(0x0),
199 fh2ESDTriggerVtx(0x0),
b5a3f310 200 fh2ESDTriggerRun(0x0),
201 fh2VtxXY(0x0),
df7848fc 202 fh1NCosmicsPerEvent(0x0),
2eded560 203 fp1RPXA(0x0),
204 fp1RPYA(0x0),
205 fp1RPXC(0x0),
206 fp1RPYC(0x0),
b360e4be 207 fp1CalibRPXA(0x0),
208 fp1CalibRPYA(0x0),
209 fp1CalibRPXC(0x0),
210 fp1CalibRPYC(0x0),
2eded560 211 fh2RPAC(0x0),
212 fh2RPAT(0x0),
213 fh2RPCT(0x0),
214 fh2XYA(0x0),
215 fh2XYC(0x0),
40445651 216 fh2RPCentrality(0x0),
2eded560 217 fh2RPACentrality(0x0),
218 fh2RPCCentrality(0x0),
3493c3a9 219 fTriggerAnalysis(0x0),
bf7b8731 220 fHistList(0x0)
221{
b5a3f310 222 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 223 DefineOutput(1,TList::Class());
224}
225
226
227
228Bool_t AliAnalysisTaskJetServices::Notify()
229{
230 //
231 // Implemented Notify() to read the cross sections
232 // and number of trials from pyxsec.root
233 //
234
235 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
236 Float_t xsection = 0;
237 Float_t ftrials = 1;
238
239 fAvgTrials = 1;
240 if(tree){
241 TFile *curfile = tree->GetCurrentFile();
242 if (!curfile) {
243 Error("Notify","No current file");
244 return kFALSE;
245 }
246 if(!fh1Xsec||!fh1Trials){
247 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
248 return kFALSE;
249 }
250 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
251 fh1Xsec->Fill("<#sigma>",xsection);
252 // construct a poor man average trials
253 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
254 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
432abb2b 255 fh1Trials->Fill("#sum{avg ntrials}",ftrials);
bf7b8731 256 }
257 return kTRUE;
258}
259
260void AliAnalysisTaskJetServices::UserCreateOutputObjects()
261{
262
263 //
264 // Create the output container
265 //
266
40445651 267 fRandomizer = new TRandom3(0);
bf7b8731 268 if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
269
270 OpenFile(1);
271 if(!fHistList)fHistList = new TList();
40440b28 272 fHistList->SetOwner();
273 PostData(1, fHistList); // post data in any case once
bf7b8731 274
275 Bool_t oldStatus = TH1::AddDirectoryStatus();
276 TH1::AddDirectory(kFALSE);
bf7b8731 277 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
278 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
279 fHistList->Add(fh1Xsec);
280
281 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
282 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
283 fHistList->Add(fh1Trials);
284
432abb2b 285 fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
286 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
287 fHistList->Add(fh1AvgTrials);
288
66fd4b2b 289 const Int_t nBinPt = 125;
bf7b8731 290 Double_t binLimitsPt[nBinPt+1];
291 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
292 if(iPt == 0){
293 binLimitsPt[iPt] = 0.0;
294 }
295 else {// 1.0
66fd4b2b 296 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
bf7b8731 297 }
298 }
299
742ee86c 300
eb1739b7 301 fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
742ee86c 302 fHistList->Add(fh1CentralityESD);
303
eb1739b7 304 fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
742ee86c 305 fHistList->Add(fh1Centrality);
306
40445651 307 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
308 fHistList->Add(fh1RP);
309
e5697f4d 310 fh2ReducedTrigger = new TH2F("fh2ReducedTrigger","red trigger;red trigger;centrality",1<<fNTrigger,-0.5,(1<<fNTrigger)-0.5,100,-0.5,99.5);
311 fHistList->Add(fh2ReducedTrigger);
3d1dba29 312
313 fh2CentralityTriggerESD = new TH2F("fh2CentralityTriggerESD",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
314 fHistList->Add(fh2CentralityTriggerESD);
315
316 fh2CentralityTrigger = new TH2F("fh2CentralityTrigger",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
317 fHistList->Add(fh2CentralityTrigger);
318
319 for(int i = 0;i<fNTrigger;++i){
5f5ae725 320 fh2CentralityTriggerESD->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
321 fh2CentralityTrigger->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
3d1dba29 322 }
323
324
d95fc15a 325 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
bf7b8731 326 fHistList->Add(fh2TriggerCount);
327
d95fc15a 328 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
bf7b8731 329 fHistList->Add(fh2ESDTriggerCount);
d95fc15a 330 const Int_t nBins = 6*kConstraints;
57ca1193 331 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
bf7b8731 332 fHistList->Add(fh2TriggerVtx);
333
57ca1193 334 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
bf7b8731 335 fHistList->Add(fh2ESDTriggerVtx);
336
337
338 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
339 fHistList->Add(fh1PtHard);
340 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
341 fHistList->Add(fh1PtHardTrials);
f2dd0695 342 fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
343 AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
f2dd0695 344 fHistList->Add(fh1SelectionInfoESD);
6d597ca2 345
346 fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
347 kTotalEventCuts,0.5,kTotalEventCuts+0.5);
348 fHistList->Add(fh1EventCutInfoESD);
349
b5a3f310 350 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
351 // 3 triggers BB BE/EB EE
352
f4132e7d 353 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 354 fHistList->Add(fh2ESDTriggerRun);
355
356 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
357 fHistList->Add(fh2VtxXY);
bf7b8731 358 // =========== Switch on Sumw2 for all histos ===========
359 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
360 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
361 if (h1){
362 h1->Sumw2();
363 continue;
364 }
365 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
366 if(hn)hn->Sumw2();
367 }
368
df7848fc 369 fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
2eded560 370 fHistList->Add(fh1NCosmicsPerEvent),
40445651 371
40445651 372
2eded560 373 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle from tracks" , 20, 0.,100., 180, 0, TMath::Pi());
40445651 374 fHistList->Add(fh2RPCentrality);
375
40445651 376
2eded560 377 fh2RPACentrality = new TH2F("fh2RPACentrality" ,"Reaction Plane Angle from vzero A" , 20, 0.,100., 180, 0, TMath::Pi());
378 fHistList->Add(fh2RPACentrality);
40445651 379
2eded560 380 fh2RPCCentrality = new TH2F("fh2RPCCentrality" ,"Reaction Plane Angle from vzero C" , 20, 0.,100., 180, 0, TMath::Pi());
381 fHistList->Add(fh2RPCCentrality);
40445651 382
2eded560 383 fh2RPAC = new TH2F("fh2RPAC" ,"Reaction Plane Angle vzero a vs c" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
384 fHistList->Add( fh2RPAC);
40445651 385
2eded560 386 fh2RPAT = new TH2F("fh2RPAT" ,"Reaction Plane Angle vzero a vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
387 fHistList->Add( fh2RPAT);
388
389 fh2RPCT = new TH2F("fh2RPCT" ,"Reaction Plane Angle vzero c vs tracks" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
390 fHistList->Add( fh2RPCT);
391
392 fh2XYA = new TH2F("fh2XYA" ,"XY vzeroa ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
393 fHistList->Add(fh2XYA);
394 fh2XYC = new TH2F("fh2XYC" ,"XY vzeroc ;X;Y;" ,100,-0.3,0.3,100,-0.3,0.3);
395 fHistList->Add(fh2XYC);
396
397 // profiles for mean
398 fp1RPXA = new TProfile("fp1RPXA","mean vzeroa x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
399 fHistList->Add(fp1RPXA);
400
401 fp1RPYA = new TProfile("fp1RPYA","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
402 fHistList->Add(fp1RPYA);
403
404
405 fp1RPXC = new TProfile("fp1RPXC","mean vzeroc x vs run number;run;x",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
406 fHistList->Add(fp1RPXC);
407
408 fp1RPYC = new TProfile("fp1RPYC","mean vzeroa y vs run number;run;y",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5);
409 fHistList->Add(fp1RPYC);
df7848fc 410
bf7b8731 411
412 TH1::AddDirectory(oldStatus);
f4132e7d 413
414 // Add an AOD branch for replication
415 if(fNonStdFile.Length()){
416 if (fDebug > 1) AliInfo("Replicating header");
417 fgAODHeader = new AliAODHeader;
418 AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
2eded560 419 if (fDebug > 1) AliInfo("Replicating vzeros");
420 fgAODVZERO = new AliAODVZERO;
421 AddAODBranch("AliAODVZERO",&fgAODVZERO,fNonStdFile.Data());
070f06d2 422 if (fDebug > 1) AliInfo("Replicating primary vertices");
070f06d2 423 fgAODVertices = new TClonesArray("AliAODVertex",3);
80a790fd 424 fgAODVertices->SetName("vertices");
070f06d2 425 AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
f4132e7d 426 }
bf7b8731 427}
428
429void AliAnalysisTaskJetServices::Init()
430{
431 //
432 // Initialization
433 //
bf7b8731 434 if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
435
436}
437
438void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
439{
440
441 //
442 // Execute analysis for current event
443 //
444
445 AliAODEvent *aod = 0;
446 AliESDEvent *esd = 0;
f4132e7d 447
448
3493c3a9 449
bf7b8731 450 AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
f4132e7d 451 AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
40445651 452 AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
f2dd0695 453 fSelectionInfoESD = 0; // reset
6d597ca2 454 fEventCutInfoESD = 0; // reset
f2dd0695 455 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
456
bf7b8731 457
f4132e7d 458 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
459
460
461
bf7b8731 462 if(fUseAODInput){
463 aod = dynamic_cast<AliAODEvent*>(InputEvent());
464 if(!aod){
465 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
466 return;
467 }
468 // fethc the header
469 }
470 else{
471 // assume that the AOD is in the general output...
472 aod = AODEvent();
473 if(!aod){
474 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
475 return;
476 }
477 esd = dynamic_cast<AliESDEvent*>(InputEvent());
478 }
ffab794c 479 if(aod&&fDebug>2){
480 aod->Print();
481 Printf("Vertices %d",aod->GetNumberOfVertices());
482 Printf("tracks %d",aod->GetNumberOfTracks());
483 Printf("jets %d",aod->GetNJets());
484 }
3493c3a9 485 fSelectionInfoESD |= kNoEventCut;
486 fEventCutInfoESD |= kNoEventCut;
6d597ca2 487
488 Bool_t esdVtxValid = false;
489 Bool_t esdVtxIn = false;
490 Bool_t aodVtxValid = false;
491 Bool_t aodVtxIn = false;
f2dd0695 492
c0997643 493 // Apply additional constraints
6d597ca2 494 Bool_t esdEventSelected = IsEventSelected(esd);
495 Bool_t esdEventPileUp = IsEventPileUp(esd);
496 Bool_t esdEventCosmic = IsEventCosmic(esd);
f2dd0695 497
6d597ca2 498 Bool_t aodEventSelected = IsEventSelected(aod);
f2dd0695 499
39e7e8ab 500 Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
20b170eb 501
f4132e7d 502
6d597ca2 503 fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
504 fh1EventCutInfoESD->Fill(fEventCutInfoESD);
505
f2dd0695 506 if(esdEventSelected) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kVertexIn;
507 if(esdEventPileUp) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsPileUp;
508 if(esdEventCosmic) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsCosmic;
509 if(physicsSelection) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kPhysicsSelection;
510
511
3d1dba29 512
513
f2dd0695 514 // here we have all selection information, fill histogram
515 for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
516 if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
517 }
6d597ca2 518 AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
c0997643 519
520 if(esd&&aod&&fDebug){
521 if(esdEventSelected&&!aodEventSelected){
522 Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
523 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
bf7b8731 524 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
c0997643 525 Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
526 vtxESD->Print();
527 Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
528 vtxAOD->Print();
bf7b8731 529 }
c0997643 530 }
531
532 // loop over all possible triggers for esd
533
aa3ba8d2 534 Float_t cent = 0;
21df3cc7 535 if(fCollisionType==kPbPb){
536 if(aod)cent = aod->GetHeader()->GetCentrality();
537 if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
538 if(cent<0)cent = 101;
539 }
40445651 540 fCentrality = cent;
541 fRPAngle = 0;
2eded560 542 fPsiVZEROA = 0;
543 fPsiVZEROC = 0;
544
40445651 545
c0997643 546 if(esd){
547 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
f4132e7d 548 esdVtxValid = IsVertexValid(vtxESD);
549 esdVtxIn = IsVertexIn(vtxESD);
d8f21f85 550 if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
40440b28 551 if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
d8f21f85 552 if(cent<=80&&esdVtxIn){
553 aodH->SetFillAOD(kTRUE);
c9931fda 554 aodH->SetFillExtension(kTRUE);
d8f21f85 555 }
556 }
557
558
57ca1193 559 Float_t zvtx = vtxESD->GetZ();
f4132e7d 560 Int_t iCl = GetEventClass(esd);
561 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
562 Bool_t cand = physicsSelection;
563
ffab794c 564 if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
d8f21f85 565 fh2ESDTriggerCount->Fill(0.,kAllTriggered);
566 fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
f4132e7d 567 if(cand){
568 fh2ESDTriggerCount->Fill(0.,kSelectedALICE);
569 fh2ESDTriggerCount->Fill(iCl,kSelectedALICE);
570 fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
571 }
572 // if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
573 if(esdVtxValid){
574 fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
575 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
576 fh2ESDTriggerVtx->Fill(iCl,zvtx);
577 if(esdVtxIn){
578 fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
579 fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
580 fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
bf7b8731 581 }
f4132e7d 582 if(cand){
583 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
584 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
585 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
c0997643 586 }
bf7b8731 587 }
742ee86c 588
c08c3ad2 589 if(cand&&esdVtxIn&&iCl<5){
f4132e7d 590 fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
591 fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
592 fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
593 fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
594 fh2ESDTriggerCount->Fill(iCl,kSelected);
595 fh2ESDTriggerCount->Fill(0.,kSelected);
596 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
742ee86c 597 if(esd->GetCentrality()){
598 Float_t tmpCent = 100;
599 tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
f5c02503 600 if(tmpCent<0)tmpCent = 101;
742ee86c 601 fh1CentralityESD->Fill(tmpCent);
3d1dba29 602 UInt_t ir = 0;
603 for(int it = 0;it<fNTrigger;++it){
604 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
605 fh2CentralityTriggerESD->Fill(tmpCent,it);
606 ir |= (1<<it);
607 }
608 }
e5697f4d 609 fh2ReducedTrigger->Fill(ir,cent);
742ee86c 610 }
f4132e7d 611 }
bf7b8731 612 }
613
ffab794c 614
615
c0997643 616 if(aod){
617 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
6d597ca2 618 aodVtxValid = IsVertexValid(vtxAOD);
619 aodVtxIn = IsVertexIn(vtxAOD);
ffab794c 620 Float_t zvtx = vtxAOD->GetZ();
621 Int_t iCl = GetEventClass(aod);
622 AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
623 Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
624 if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
d8f21f85 625 fh2TriggerCount->Fill(0.,kAllTriggered);
626 fh2TriggerCount->Fill(iCl,kAllTriggered);
ffab794c 627 if(cand){
628 fh2TriggerCount->Fill(0.,kSelectedALICE);
629 fh2TriggerCount->Fill(iCl,kSelectedALICE);
630 fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
631 }
632 if(aodVtxValid){
633 fh2TriggerCount->Fill(0.,kTriggeredVertex);
634 fh2TriggerCount->Fill(iCl,kTriggeredVertex);
635 fh2TriggerVtx->Fill(iCl,zvtx);
636 if(aodVtxIn){
637 fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
638 fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
639 fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
c0997643 640 }
ffab794c 641 if(cand){
642 fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
643 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
644 fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
645 }
646 }
90a006c2 647
c08c3ad2 648 if(cand&&aodVtxIn&&iCl<5){
ffab794c 649 fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
650 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
651 fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
652 fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
653 fh2TriggerCount->Fill(iCl,kSelected);
654 fh2TriggerCount->Fill(0.,kSelected);
742ee86c 655 fh1Centrality->Fill(cent);
3d1dba29 656 for(int it = 0;it<fNTrigger;++it){
657 if(fInputHandler->IsEventSelected()&fTriggerBit[it])fh2CentralityTrigger->Fill(cent,it);
658 }
742ee86c 659 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
90a006c2 660 if(aodH&&cand&&fFilterAODCollisions&&!esd){
21df3cc7 661 if(fCentrality<=80&&aodVtxIn){
90a006c2 662 aodH->SetFillAOD(kTRUE);
663 aodH->SetFillExtension(kTRUE);
664 }
665 }
666
40445651 667 TList recTracks;
668 GetListOfTracks(&recTracks);
2eded560 669 CalculateReactionPlaneAngleVZERO(aod);
d3a3f33d 670 fRPAngle = aod->GetHeader()->GetEventplane();
671 fh1RP->Fill(fRPAngle);
672 fh2RPCentrality->Fill(fCentrality,fRPAngle);
2eded560 673 fh2RPACentrality->Fill(fCentrality,fPsiVZEROA);
674 fh2RPCCentrality->Fill(fCentrality,fPsiVZEROC);
675 fh2RPAC->Fill(fPsiVZEROA,fPsiVZEROC);
676 fh2RPAT->Fill(fPsiVZEROA,fRPAngle);
677 fh2RPCT->Fill(fPsiVZEROC,fRPAngle);
678 if(fRPMethod==kRPTracks)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
679 else if(fRPMethod==kRPVZEROA)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
680 else if(fRPMethod==kRPVZEROC)AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fPsiVZEROA); // set slection to false
d3a3f33d 681
21df3cc7 682 if(fUseAODInput&&fCentrality<=80){
221d8aa9 683 if(fFilterAODCollisions&&aod){
742ee86c 684 aodH->SetFillAOD(kTRUE);
eedaa08e 685 aodH->SetFillExtension(kTRUE);
221d8aa9 686 }
687 }
221d8aa9 688 }
c0997643 689 }
bf7b8731 690
691 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
692
693
694 Double_t ptHard = 0;
695 Double_t nTrials = 1; // Trials for MC trigger
696
432abb2b 697 fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
bf7b8731 698 AliMCEvent* mcEvent = MCEvent();
699 // AliStack *pStack = 0;
700 if(mcEvent){
701 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
702 if(pythiaGenHeader){
703 nTrials = pythiaGenHeader->Trials();
704 ptHard = pythiaGenHeader->GetPtHard();
705 int iProcessType = pythiaGenHeader->ProcessType();
706 // 11 f+f -> f+f
707 // 12 f+barf -> f+barf
708 // 13 f+barf -> g+g
709 // 28 f+g -> f+g
710 // 53 g+g -> f+barf
711 // 68 g+g -> g+g
712 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
713 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
714 fh1PtHard->Fill(ptHard);
715 fh1PtHardTrials->Fill(ptHard,nTrials);
716
717 }// if pythia gen header
718 }
719
720 // trigger selection
721
f4132e7d 722 // replication of
723 if(fNonStdFile.Length()&&aod){
724 if (fgAODHeader){
725 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
2eded560 726 Double_t q[kRPMethods] = {fRPAngle,fPsiVZEROA,fPsiVZEROC};
727 fgAODHeader->SetQTheta(q,kRPMethods);
728 }
729 if (fgAODVZERO){
730 *fgAODVZERO = *(dynamic_cast<AliAODVZERO*>(aod->GetVZEROData()));
f4132e7d 731 }
80a790fd 732 if(fgAODVertices){
80a790fd 733 fgAODVertices->Delete();
734 TClonesArray &vertices = *fgAODVertices;
80a790fd 735 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
070f06d2 736 // we only use some basic information,
737
738
739 Double_t pos[3];
740 Double_t covVtx[6];
741
742 vtxAOD->GetXYZ(pos); // position
743 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
744 Int_t jVertices = 0;
745 AliAODVertex * vtx = new(vertices[jVertices++])
746 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
747 vtx->SetName(vtxAOD->GetName());
748 vtx->SetTitle(vtxAOD->GetTitle());
070f06d2 749 TString vtitle = vtxAOD->GetTitle();
d8f21f85 750 vtx->SetNContributors(vtxAOD->GetNContributors());
070f06d2 751
752 // Add SPD "main" vertex
753 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
754 vtxS->GetXYZ(pos); // position
755 vtxS->GetCovMatrix(covVtx); //covariance matrix
756 AliAODVertex * mVSPD = new(vertices[jVertices++])
757 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
758 mVSPD->SetName(vtxS->GetName());
759 mVSPD->SetTitle(vtxS->GetTitle());
760 mVSPD->SetNContributors(vtxS->GetNContributors());
761
762 // Add tpc only vertex
763 if(esd){
764 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
765 vtxT->GetXYZ(pos); // position
766 vtxT->GetCovMatrix(covVtx); //covariance matrix
767 AliAODVertex * mVTPC = new(vertices[jVertices++])
768 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
769 mVTPC->SetName(vtxT->GetName());
770 mVTPC->SetTitle(vtxT->GetTitle());
771 mVTPC->SetNContributors(vtxT->GetNContributors());
772 }
80a790fd 773 }
f4132e7d 774 }
775
bf7b8731 776 PostData(1, fHistList);
80a790fd 777}
bf7b8731 778
6d597ca2 779Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
c0997643 780 if(!esd)return kFALSE;
781 const AliESDVertex *vtx = esd->GetPrimaryVertex();
6d597ca2 782 return IsVertexIn(vtx); // vertex in calls vertex valid
783}
784
80a790fd 785AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
786 if(fgAODVertices){
787 fgAODVertices->Delete();
788 delete fgAODVertices;
789 }
40445651 790 delete fRandomizer;
070f06d2 791 if(fgAODHeader)delete fgAODHeader;
b360e4be 792 if(fgAODVZERO)delete fgAODVZERO;
793 delete fp1CalibRPXA;
794 delete fp1CalibRPYA;
795 delete fp1CalibRPXC;
796 delete fp1CalibRPYC;
3d1dba29 797 delete [] fTriggerBit;
798 delete [] fTriggerName;
b360e4be 799
80a790fd 800}
801
6d597ca2 802
b360e4be 803void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfile *xc, TProfile *yc){
804
805 if(xa){
806 if(fp1CalibRPXA)delete fp1CalibRPXA;
807 fp1CalibRPXA = (TProfile*)xa->Clone(Form("%sCalib",xa->GetName()));
808 }
809 else{
810 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
811 }
812 if(ya){
813 if(fp1CalibRPYA)delete fp1CalibRPYA;
814 fp1CalibRPYA = (TProfile*)ya->Clone(Form("%sCalib",ya->GetName()));
815 }
816 else{
817 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
818 }
819 if(xc){
820 if(fp1CalibRPXC)delete fp1CalibRPXC;
821 fp1CalibRPXC = (TProfile*)xc->Clone(Form("%sCalib",xc->GetName()));
822 }
823 else{
824 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
825 }
826 if(ya){
827 if(fp1CalibRPYC)delete fp1CalibRPYC;
828 fp1CalibRPYC = (TProfile*)yc->Clone(Form("%sCalib",yc->GetName()));
829 }
830 else{
831 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
832 }
833}
834
6d597ca2 835Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
836 if(!aod)return kFALSE;
837 const AliAODVertex *vtx = aod->GetPrimaryVertex();
838 return IsVertexIn(vtx); // VertexIn calls VertexValid
839}
840
841Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
842
c0997643 843 // We can check the type of the vertex by its name and title
844 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
845 // vtx name title
846 // Tracks PrimaryVertex VertexerTracksNoConstraint
6d597ca2 847 // Tracks PrimaryVertex VertexerTracksWithConstraint
c0997643 848 // TPC TPCVertex VertexerTracksNoConstraint
6d597ca2 849 // TPC TPCVertex VertexerTracksWithConstraint
c0997643 850 // SPD SPDVertex vertexer: 3D
851 // SPD SPDVertex vertexer: Z
852
d3a3f33d 853 if(!vtx){
854 Printf("%s:%d No ESD vertex found",(char*)__FILE__,__LINE__);
855 return kFALSE;
856 }
6d597ca2 857 Int_t nCont = vtx->GetNContributors();
6d597ca2 858 if(nCont>=1){
859 fEventCutInfoESD |= kContributorsCut1;
860 if(nCont>=2){
861 fEventCutInfoESD |= kContributorsCut2;
862 if(nCont>=3){
863 fEventCutInfoESD |= kContributorsCut3;
864 }
865 }
866 }
867
868 if(nCont<3)return kFALSE;
c0997643 869
6d597ca2 870 // do not want tpc only primary vertex
871 TString vtxName(vtx->GetName());
872 if(vtxName.Contains("TPCVertex")){
873 fEventCutInfoESD |= kVertexTPC;
874 return kFALSE;
875 }
876 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
877 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
c0997643 878
c0997643 879
6d597ca2 880 TString vtxTitle(vtx->GetTitle());
881 if(vtxTitle.Contains("vertexer: Z")){
882 if(vtx->GetDispersion()>0.02)return kFALSE;
883 }
884 fEventCutInfoESD |= kSPDDispersionCut;
885 return kTRUE;
886}
887
888
889Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
890
891 // We can check the type of the vertex by its name and title
892 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
893 // vtx name title
894 // Tracks PrimaryVertex VertexerTracksNoConstraint
895 // TPC TPCVertex VertexerTracksNoConstraint
896 // SPD SPDVertex vertexer: 3D
897 // SPD SPDVertex vertexer: Z
d8f21f85 898
d3a3f33d 899 if(!vtx){
900 Printf("%s:%d No AOD vertex found",(char*)__FILE__,__LINE__);
901 return kFALSE;
902 }
903
9410e6a2 904 TString vtxName(vtx->GetName());
d8f21f85 905 if(fDebug){
906 Printf(" n contrib %d",vtx->GetNContributors());
9410e6a2 907 Printf("vtxName: %s",vtxName.Data());
d8f21f85 908 vtx->Print();
909 }
6d597ca2 910
d8f21f85 911 // if(vtx->GetNContributors()<3)return kFALSE;
c0997643 912 // do not want tpc only primary vertex
9410e6a2 913
c0997643 914 if(vtxName.Contains("TPCVertex"))return kFALSE;
6d597ca2 915
916 // no dispersion yet...
917 /*
918 TString vtxTitle(vtx->GetTitle());
919 if(vtxTitle.Contains("vertexer: Z")){
920 if(vtx->GetDispersion()>0.02)return kFALSE;
921 }
922 */
923 return kTRUE;
924}
925
926
927Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
928
929 if(!IsVertexValid(vtx))return kFALSE;
930
c0997643 931 Float_t zvtx = vtx->GetZ();
932 Float_t yvtx = vtx->GetY();
933 Float_t xvtx = vtx->GetX();
934
6d597ca2 935 xvtx -= fVtxXMean;
936 yvtx -= fVtxYMean;
937 zvtx -= fVtxZMean;
938
939
940
941 if(TMath::Abs(zvtx)>fVtxZCut){
942 return kFALSE;
943 }
944 fEventCutInfoESD |= kVertexZCut;
945 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
946 if(r2>(fVtxRCut*fVtxRCut)){
947 return kFALSE;
948 }
949 fEventCutInfoESD |= kVertexRCut;
950 return kTRUE;
951}
952
953
954Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
955
956 if(!IsVertexValid(vtx))return kFALSE;
957
958 Float_t zvtx = vtx->GetZ();
959 Float_t yvtx = vtx->GetY();
960 Float_t xvtx = vtx->GetX();
c0997643 961
962 xvtx -= fVtxXMean;
963 yvtx -= fVtxYMean;
964 zvtx -= fVtxZMean;
965
966 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
967
6d597ca2 968 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
969 return vertexIn;
c0997643 970}
971
6d597ca2 972Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
f2dd0695 973 if(!esd)return kFALSE;
974 return esd->IsPileupFromSPD();
975}
976
6d597ca2 977Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
f2dd0695 978 if(!esd)return kFALSE;
979 // add track cuts for which we look for cosmics...
df7848fc 980
981 Bool_t isCosmic = kFALSE;
982 Int_t nTracks = esd->GetNumberOfTracks();
983 Int_t nCosmicCandidates = 0;
984
985 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
986 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
987 if (!track1) continue;
75946d5d 988 UInt_t status1 = track1->GetStatus();
989 //If track is ITS stand alone track, skip the track
990 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
df7848fc 991 if(track1->Pt()<fPtMinCosmic) continue;
992 //Start 2nd track loop to look for correlations
993 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
994 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
995 if(!track2) continue;
75946d5d 996 UInt_t status2 = track2->GetStatus();
997 //If track is ITS stand alone track, skip the track
998 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
df7848fc 999 if(track2->Pt()<fPtMinCosmic) continue;
1000 //Check if back-to-back
1001 Double_t mom1[3],mom2[3];
1002 track1->GetPxPyPz(mom1);
1003 track2->GetPxPyPz(mom2);
1004 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
1005 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
1006 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
1007 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
1008
1009 Float_t deltaPhi = track1->Phi()-track2->Phi();
1010 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
1011
1012 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
1013 if(rIsol<fRIsolMinCosmic) continue;
1014
1015 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
75946d5d 1016 nCosmicCandidates+=1;
df7848fc 1017 isCosmic = kTRUE;
1018 }
1019
1020 }
1021 }
1022
1023 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
1024
1025 return isCosmic;
f2dd0695 1026}
1027
1028
f4132e7d 1029Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
1030
aa3ba8d2 1031 Float_t cent = 0;
21df3cc7 1032 if(fCollisionType==kPbPb){
1033 if(esd->GetCentrality()){
1034 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
1035 }
1036 if(cent<0)cent = 101;
1037 if(cent>80||cent<0)return 5;
1038 if(cent>50)return 4;
1039 if(cent>30)return 3;
1040 if(cent>10)return 2;
1041 return 1;
f4132e7d 1042 }
f4132e7d 1043 return 1;
f4132e7d 1044}
f2dd0695 1045
ffab794c 1046
1047Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
1048
21df3cc7 1049 if(fCollisionType==kPbPb){
1050 Float_t cent = aod->GetHeader()->GetCentrality();
1051 if(cent>80||cent<0)return 5;
1052 if(cent>50)return 4;
1053 if(cent>30)return 3;
1054 if(cent>10)return 2;
1055 return 1;
1056 }
ffab794c 1057 return 1;
1058
1059}
1060
1061
bf7b8731 1062void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
1063{
1064 // Terminate analysis
1065 //
bf7b8731 1066}
f4132e7d 1067
2eded560 1068Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngleVZERO(AliAODEvent *aod){
40445651 1069
2eded560 1070 // const Double_t arr_eta[11]={-3.7, -3.2, -2.7, -2.2, -1.7,0, 2.8, 3.4, 3.9, 4.5,5.1};
1071 if(!aod)return kFALSE;
1072 static bool bFirst = true;
1073 static Double_t v0phi[64] = {0,};
1074
1075 if(bFirst){
1076 int is=0;
1077 for(int iArm = 0; iArm<2; iArm++){
1078 for(int iRing = 0; iRing<4 ; iRing++){
1079 for(int iSec = 0; iSec<8 ; iSec++){
1080 v0phi[is] = 22.5 + 45. * iSec;
1081 v0phi[is] *= TMath::Pi()/180;
1082 // cout<< is <<" "<< v0phi[is]<<endl;
1083 is++;
1084 }
1085 }
40445651 1086 }
2eded560 1087 bFirst = false;
40445651 1088 }
1089
2eded560 1090 //
1091 const AliAODVZERO *aodVZERO = aod->GetVZEROData();
1092 Double_t numYZNA = 0,numXZNA = 0,sumZNA = 0;
1093 Double_t meanXA = 0,meanYA = 0;
1094
1095 Double_t numYZNC = 0,numXZNC = 0,sumZNC = 0;
1096 Double_t meanXC = 0,meanYC = 0;
1097
b360e4be 1098
1099
1100 static Int_t iOldRun = -1;
1101 static Int_t iFoundBin = -1;
2e90089b 1102
1103 if(aod->GetRunNumber()!=iOldRun&&(fp1CalibRPYA)){
b360e4be 1104 // search only or the bin in case of new runs
1105 iFoundBin = -1;
1106 Int_t ib = fp1CalibRPYA->FindBin(aod->GetRunNumber());
1107 Float_t err = fp1CalibRPYA->GetBinError(ib);
1108 if(err>0){// value can be zero...
1109 iFoundBin = ib;
1110 }
1111 else{
1112 Int_t ibLo = ib-1;
1113 Int_t ibUp = ib+1;
1114 while(iFoundBin<0&&(ibLo>0||ibUp<=fp1CalibRPYA->GetNbinsX())){
1115 err = fp1CalibRPYA->GetBinError(ibLo);
1116 if(err>0){
1117 iFoundBin = ibLo;
1118 }
1119 else{
1120 err = fp1CalibRPYA->GetBinError(ibUp);
1121 if(err>0)iFoundBin = ibUp;
1122 }
1123 ibUp++;
1124 ibLo--;
1125 }
1126 }
1127 iOldRun = aod->GetRunNumber();
1128 }
1129
2e90089b 1130 if(fDebug)Printf("%s:%d iFoundBin %d",(char*)__FILE__,__LINE__,iFoundBin);
b360e4be 1131
d889ce29 1132 if(iFoundBin>0&&(fp1CalibRPYA)){
b360e4be 1133 meanXA = fp1CalibRPXA->GetBinContent(iFoundBin);
1134 meanYA = fp1CalibRPYA->GetBinContent(iFoundBin);
1135 meanXC = fp1CalibRPXC->GetBinContent(iFoundBin);
1136 meanYC = fp1CalibRPYC->GetBinContent(iFoundBin);
1137 }
1138
2e90089b 1139 if(fDebug)Printf("%s:%d iFoundBin %1.3E %1.3E %1.3E %1.3E",(char*)__FILE__,__LINE__,meanXA,meanYA,meanXC,meanYC);
b360e4be 1140
2eded560 1141 for (int i=0; i<64; i++) {
1142 Double_t mult = aodVZERO->GetMultiplicity(i);
1143 Double_t phi = v0phi[i];
1144 if (mult>0) {
1145 if (i<32) { //C-side
1146 Double_t wZNC= mult;
1147 numYZNC += sin(2.*phi)*wZNC;
1148 numXZNC += cos(2.*phi)*wZNC;
1149 sumZNC+=wZNC;
1150 }
1151 else if(i>31){ //A-side
1152 Double_t wZNA=mult;
1153 numYZNA += sin(2.*phi)*wZNA;
1154 numXZNA += cos(2.*phi)*wZNA;
1155 sumZNA+=wZNA;
1156 }
1157 }// mult>0
1158 }// 64 sectors
1159
40445651 1160
a374d56f 1161 Double_t XC = 0;
1162 Double_t YC = 0;
1163 Double_t XA = 0;
1164 Double_t YA = 0;
1165
1166 if(sumZNC!=0){
2eded560 1167
a374d56f 1168 XC = numXZNC/sumZNC;
1169 YC = numYZNC/sumZNC;
6622701d 1170 fPsiVZEROC = 0.5*TMath::ATan2(YC-meanYC, XC-meanXC);
a374d56f 1171 if(fPsiVZEROC>TMath::Pi()){fPsiVZEROC-=TMath::Pi();}
1172 if(fPsiVZEROC<0){fPsiVZEROC+=TMath::Pi();}
1173
1174
1175 fh2XYC->Fill(XC-meanXC,YC-meanYC); // control
1176 fp1RPXC->Fill(aod->GetRunNumber(),XC);
1177 fp1RPYC->Fill(aod->GetRunNumber(),YC);
1178
2eded560 1179
a374d56f 1180 }
40445651 1181
a374d56f 1182 if(sumZNA!=0){
1183 XA = numXZNA/sumZNA;
1184 YA = numYZNA/sumZNA;
1185 fPsiVZEROA = 0.5*TMath::ATan2(YA-meanYA, XA-meanXA);
1186 if(fPsiVZEROA>TMath::Pi()){fPsiVZEROA-=TMath::Pi();}
1187 if(fPsiVZEROA<0){fPsiVZEROA+=TMath::Pi();}
1188
1189 fh2XYA->Fill(XA-meanXA,YA-meanYA); // control
1190 fp1RPXA->Fill(aod->GetRunNumber(),XA);
1191 fp1RPYA->Fill(aod->GetRunNumber(),YA);
1192
1193 }
40445651 1194 return kTRUE;
40445651 1195
40445651 1196}
1197
40445651 1198Int_t AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1199 Int_t iCount = 0;
1200 AliAODEvent *aod = 0;
1201 if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1202 else aod = AODEvent();
1203 if(!aod){
1204 return iCount;
1205 }
1206 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1207 AliAODTrack *tr = aod->GetTrack(it);
1208 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1209 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1210 if(tr->Pt()<fMinTrackPt)continue;
1211 list->Add(tr);
1212 iCount++;
1213 }
1214 list->Sort();
1215 return iCount;
1216
1217}
1218
3d1dba29 1219void AliAnalysisTaskJetServices::SetNTrigger(Int_t n){
1220 if(n>0){
1221 fNTrigger = n;
1222 delete [] fTriggerName;
1223 fTriggerName = new TString [fNTrigger];
1224 delete [] fTriggerBit;fTriggerBit = 0;
1225 fTriggerBit = new UInt_t [fNTrigger];
1226 }
1227 else{
1228 fNTrigger = 0;
1229 }
1230}
1231
1232void AliAnalysisTaskJetServices::SetTrigger(Int_t i,UInt_t it,const char* c){
1233 if(i<fNTrigger){
1234 Printf("%d",it);
1235 Printf("%p",c);
1236 Printf("%s",c);
1237 Printf("%p",&fTriggerName[i]);
1238
1239 fTriggerBit[i] = it;
1240 fTriggerName[i] = c; // placement new
1241 // new(&fTriggerName[i]) TString(c); // placement new
1242 Printf("%s",fTriggerName[i].Data());
1243
1244 }
1245}