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