Protection added against files from DQM were found in the DAQ FXS, but no raw tag...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / 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);
221d8aa9 657 }
658 }
221d8aa9 659 }
c0997643 660 }
bf7b8731 661
662 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
663
664
665 Double_t ptHard = 0;
666 Double_t nTrials = 1; // Trials for MC trigger
667
668 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
669 AliMCEvent* mcEvent = MCEvent();
670 // AliStack *pStack = 0;
671 if(mcEvent){
672 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
673 if(pythiaGenHeader){
674 nTrials = pythiaGenHeader->Trials();
675 ptHard = pythiaGenHeader->GetPtHard();
676 int iProcessType = pythiaGenHeader->ProcessType();
677 // 11 f+f -> f+f
678 // 12 f+barf -> f+barf
679 // 13 f+barf -> g+g
680 // 28 f+g -> f+g
681 // 53 g+g -> f+barf
682 // 68 g+g -> g+g
683 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
684 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
685 fh1PtHard->Fill(ptHard);
686 fh1PtHardTrials->Fill(ptHard,nTrials);
687
688 }// if pythia gen header
689 }
690
691 // trigger selection
692
f4132e7d 693 // replication of
694 if(fNonStdFile.Length()&&aod){
695 if (fgAODHeader){
696 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
2eded560 697 Double_t q[kRPMethods] = {fRPAngle,fPsiVZEROA,fPsiVZEROC};
698 fgAODHeader->SetQTheta(q,kRPMethods);
699 }
700 if (fgAODVZERO){
701 *fgAODVZERO = *(dynamic_cast<AliAODVZERO*>(aod->GetVZEROData()));
f4132e7d 702 }
80a790fd 703 if(fgAODVertices){
80a790fd 704 fgAODVertices->Delete();
705 TClonesArray &vertices = *fgAODVertices;
80a790fd 706 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
070f06d2 707 // we only use some basic information,
708
709
710 Double_t pos[3];
711 Double_t covVtx[6];
712
713 vtxAOD->GetXYZ(pos); // position
714 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
715 Int_t jVertices = 0;
716 AliAODVertex * vtx = new(vertices[jVertices++])
717 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
718 vtx->SetName(vtxAOD->GetName());
719 vtx->SetTitle(vtxAOD->GetTitle());
070f06d2 720 TString vtitle = vtxAOD->GetTitle();
d8f21f85 721 vtx->SetNContributors(vtxAOD->GetNContributors());
070f06d2 722
723 // Add SPD "main" vertex
724 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
725 vtxS->GetXYZ(pos); // position
726 vtxS->GetCovMatrix(covVtx); //covariance matrix
727 AliAODVertex * mVSPD = new(vertices[jVertices++])
728 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
729 mVSPD->SetName(vtxS->GetName());
730 mVSPD->SetTitle(vtxS->GetTitle());
731 mVSPD->SetNContributors(vtxS->GetNContributors());
732
733 // Add tpc only vertex
734 if(esd){
735 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
736 vtxT->GetXYZ(pos); // position
737 vtxT->GetCovMatrix(covVtx); //covariance matrix
738 AliAODVertex * mVTPC = new(vertices[jVertices++])
739 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
740 mVTPC->SetName(vtxT->GetName());
741 mVTPC->SetTitle(vtxT->GetTitle());
742 mVTPC->SetNContributors(vtxT->GetNContributors());
743 }
80a790fd 744 }
f4132e7d 745 }
746
bf7b8731 747 PostData(1, fHistList);
80a790fd 748}
bf7b8731 749
6d597ca2 750Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
c0997643 751 if(!esd)return kFALSE;
752 const AliESDVertex *vtx = esd->GetPrimaryVertex();
6d597ca2 753 return IsVertexIn(vtx); // vertex in calls vertex valid
754}
755
80a790fd 756AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
757 if(fgAODVertices){
758 fgAODVertices->Delete();
759 delete fgAODVertices;
760 }
40445651 761 delete fRandomizer;
070f06d2 762 if(fgAODHeader)delete fgAODHeader;
b360e4be 763 if(fgAODVZERO)delete fgAODVZERO;
764 delete fp1CalibRPXA;
765 delete fp1CalibRPYA;
766 delete fp1CalibRPXC;
767 delete fp1CalibRPYC;
768
80a790fd 769}
770
6d597ca2 771
b360e4be 772void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfile *xc, TProfile *yc){
773
774 if(xa){
775 if(fp1CalibRPXA)delete fp1CalibRPXA;
776 fp1CalibRPXA = (TProfile*)xa->Clone(Form("%sCalib",xa->GetName()));
777 }
778 else{
779 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
780 }
781 if(ya){
782 if(fp1CalibRPYA)delete fp1CalibRPYA;
783 fp1CalibRPYA = (TProfile*)ya->Clone(Form("%sCalib",ya->GetName()));
784 }
785 else{
786 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
787 }
788 if(xc){
789 if(fp1CalibRPXC)delete fp1CalibRPXC;
790 fp1CalibRPXC = (TProfile*)xc->Clone(Form("%sCalib",xc->GetName()));
791 }
792 else{
793 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
794 }
795 if(ya){
796 if(fp1CalibRPYC)delete fp1CalibRPYC;
797 fp1CalibRPYC = (TProfile*)yc->Clone(Form("%sCalib",yc->GetName()));
798 }
799 else{
800 Printf("%s:%d centroid histogram is 0x0",(char*)__FILE__,__LINE__);
801 }
802}
803
6d597ca2 804Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
805 if(!aod)return kFALSE;
806 const AliAODVertex *vtx = aod->GetPrimaryVertex();
807 return IsVertexIn(vtx); // VertexIn calls VertexValid
808}
809
810Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
811
c0997643 812 // We can check the type of the vertex by its name and title
813 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
814 // vtx name title
815 // Tracks PrimaryVertex VertexerTracksNoConstraint
6d597ca2 816 // Tracks PrimaryVertex VertexerTracksWithConstraint
c0997643 817 // TPC TPCVertex VertexerTracksNoConstraint
6d597ca2 818 // TPC TPCVertex VertexerTracksWithConstraint
c0997643 819 // SPD SPDVertex vertexer: 3D
820 // SPD SPDVertex vertexer: Z
821
d3a3f33d 822 if(!vtx){
823 Printf("%s:%d No ESD vertex found",(char*)__FILE__,__LINE__);
824 return kFALSE;
825 }
6d597ca2 826 Int_t nCont = vtx->GetNContributors();
6d597ca2 827 if(nCont>=1){
828 fEventCutInfoESD |= kContributorsCut1;
829 if(nCont>=2){
830 fEventCutInfoESD |= kContributorsCut2;
831 if(nCont>=3){
832 fEventCutInfoESD |= kContributorsCut3;
833 }
834 }
835 }
836
837 if(nCont<3)return kFALSE;
c0997643 838
6d597ca2 839 // do not want tpc only primary vertex
840 TString vtxName(vtx->GetName());
841 if(vtxName.Contains("TPCVertex")){
842 fEventCutInfoESD |= kVertexTPC;
843 return kFALSE;
844 }
845 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
846 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
c0997643 847
c0997643 848
6d597ca2 849 TString vtxTitle(vtx->GetTitle());
850 if(vtxTitle.Contains("vertexer: Z")){
851 if(vtx->GetDispersion()>0.02)return kFALSE;
852 }
853 fEventCutInfoESD |= kSPDDispersionCut;
854 return kTRUE;
855}
856
857
858Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
859
860 // We can check the type of the vertex by its name and title
861 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
862 // vtx name title
863 // Tracks PrimaryVertex VertexerTracksNoConstraint
864 // TPC TPCVertex VertexerTracksNoConstraint
865 // SPD SPDVertex vertexer: 3D
866 // SPD SPDVertex vertexer: Z
d8f21f85 867
d3a3f33d 868 if(!vtx){
869 Printf("%s:%d No AOD vertex found",(char*)__FILE__,__LINE__);
870 return kFALSE;
871 }
872
873
d8f21f85 874 if(fDebug){
875 Printf(" n contrib %d",vtx->GetNContributors());
876 vtx->Print();
877 }
6d597ca2 878
d8f21f85 879 // if(vtx->GetNContributors()<3)return kFALSE;
c0997643 880 // do not want tpc only primary vertex
881 TString vtxName(vtx->GetName());
882 if(vtxName.Contains("TPCVertex"))return kFALSE;
6d597ca2 883
884 // no dispersion yet...
885 /*
886 TString vtxTitle(vtx->GetTitle());
887 if(vtxTitle.Contains("vertexer: Z")){
888 if(vtx->GetDispersion()>0.02)return kFALSE;
889 }
890 */
891 return kTRUE;
892}
893
894
895Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
896
897 if(!IsVertexValid(vtx))return kFALSE;
898
c0997643 899 Float_t zvtx = vtx->GetZ();
900 Float_t yvtx = vtx->GetY();
901 Float_t xvtx = vtx->GetX();
902
6d597ca2 903 xvtx -= fVtxXMean;
904 yvtx -= fVtxYMean;
905 zvtx -= fVtxZMean;
906
907
908
909 if(TMath::Abs(zvtx)>fVtxZCut){
910 return kFALSE;
911 }
912 fEventCutInfoESD |= kVertexZCut;
913 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
914 if(r2>(fVtxRCut*fVtxRCut)){
915 return kFALSE;
916 }
917 fEventCutInfoESD |= kVertexRCut;
918 return kTRUE;
919}
920
921
922Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
923
924 if(!IsVertexValid(vtx))return kFALSE;
925
926 Float_t zvtx = vtx->GetZ();
927 Float_t yvtx = vtx->GetY();
928 Float_t xvtx = vtx->GetX();
c0997643 929
930 xvtx -= fVtxXMean;
931 yvtx -= fVtxYMean;
932 zvtx -= fVtxZMean;
933
934 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
935
6d597ca2 936 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
937 return vertexIn;
c0997643 938}
939
6d597ca2 940Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
f2dd0695 941 if(!esd)return kFALSE;
942 return esd->IsPileupFromSPD();
943}
944
6d597ca2 945Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
f2dd0695 946 if(!esd)return kFALSE;
947 // add track cuts for which we look for cosmics...
df7848fc 948
949 Bool_t isCosmic = kFALSE;
950 Int_t nTracks = esd->GetNumberOfTracks();
951 Int_t nCosmicCandidates = 0;
952
953 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
954 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
955 if (!track1) continue;
75946d5d 956 UInt_t status1 = track1->GetStatus();
957 //If track is ITS stand alone track, skip the track
958 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
df7848fc 959 if(track1->Pt()<fPtMinCosmic) continue;
960 //Start 2nd track loop to look for correlations
961 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
962 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
963 if(!track2) continue;
75946d5d 964 UInt_t status2 = track2->GetStatus();
965 //If track is ITS stand alone track, skip the track
966 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
df7848fc 967 if(track2->Pt()<fPtMinCosmic) continue;
968 //Check if back-to-back
969 Double_t mom1[3],mom2[3];
970 track1->GetPxPyPz(mom1);
971 track2->GetPxPyPz(mom2);
972 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
973 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
974 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
975 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
976
977 Float_t deltaPhi = track1->Phi()-track2->Phi();
978 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
979
980 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
981 if(rIsol<fRIsolMinCosmic) continue;
982
983 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
75946d5d 984 nCosmicCandidates+=1;
df7848fc 985 isCosmic = kTRUE;
986 }
987
988 }
989 }
990
991 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
992
993 return isCosmic;
f2dd0695 994}
995
996
f4132e7d 997Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
998
aa3ba8d2 999 Float_t cent = 0;
21df3cc7 1000 if(fCollisionType==kPbPb){
1001 if(esd->GetCentrality()){
1002 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
1003 }
1004 if(cent<0)cent = 101;
1005 if(cent>80||cent<0)return 5;
1006 if(cent>50)return 4;
1007 if(cent>30)return 3;
1008 if(cent>10)return 2;
1009 return 1;
f4132e7d 1010 }
f4132e7d 1011 return 1;
f4132e7d 1012}
f2dd0695 1013
ffab794c 1014
1015Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
1016
21df3cc7 1017 if(fCollisionType==kPbPb){
1018 Float_t cent = aod->GetHeader()->GetCentrality();
1019 if(cent>80||cent<0)return 5;
1020 if(cent>50)return 4;
1021 if(cent>30)return 3;
1022 if(cent>10)return 2;
1023 return 1;
1024 }
ffab794c 1025 return 1;
1026
1027}
1028
1029
bf7b8731 1030void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
1031{
1032 // Terminate analysis
1033 //
bf7b8731 1034}
f4132e7d 1035
2eded560 1036Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngleVZERO(AliAODEvent *aod){
40445651 1037
2eded560 1038 // 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};
1039 if(!aod)return kFALSE;
1040 static bool bFirst = true;
1041 static Double_t v0phi[64] = {0,};
1042
1043 if(bFirst){
1044 int is=0;
1045 for(int iArm = 0; iArm<2; iArm++){
1046 for(int iRing = 0; iRing<4 ; iRing++){
1047 for(int iSec = 0; iSec<8 ; iSec++){
1048 v0phi[is] = 22.5 + 45. * iSec;
1049 v0phi[is] *= TMath::Pi()/180;
1050 // cout<< is <<" "<< v0phi[is]<<endl;
1051 is++;
1052 }
1053 }
40445651 1054 }
2eded560 1055 bFirst = false;
40445651 1056 }
1057
2eded560 1058 //
1059 const AliAODVZERO *aodVZERO = aod->GetVZEROData();
1060 Double_t numYZNA = 0,numXZNA = 0,sumZNA = 0;
1061 Double_t meanXA = 0,meanYA = 0;
1062
1063 Double_t numYZNC = 0,numXZNC = 0,sumZNC = 0;
1064 Double_t meanXC = 0,meanYC = 0;
1065
b360e4be 1066
1067
1068 static Int_t iOldRun = -1;
1069 static Int_t iFoundBin = -1;
2e90089b 1070
1071 if(aod->GetRunNumber()!=iOldRun&&(fp1CalibRPYA)){
b360e4be 1072 // search only or the bin in case of new runs
1073 iFoundBin = -1;
1074 Int_t ib = fp1CalibRPYA->FindBin(aod->GetRunNumber());
1075 Float_t err = fp1CalibRPYA->GetBinError(ib);
1076 if(err>0){// value can be zero...
1077 iFoundBin = ib;
1078 }
1079 else{
1080 Int_t ibLo = ib-1;
1081 Int_t ibUp = ib+1;
1082 while(iFoundBin<0&&(ibLo>0||ibUp<=fp1CalibRPYA->GetNbinsX())){
1083 err = fp1CalibRPYA->GetBinError(ibLo);
1084 if(err>0){
1085 iFoundBin = ibLo;
1086 }
1087 else{
1088 err = fp1CalibRPYA->GetBinError(ibUp);
1089 if(err>0)iFoundBin = ibUp;
1090 }
1091 ibUp++;
1092 ibLo--;
1093 }
1094 }
1095 iOldRun = aod->GetRunNumber();
1096 }
1097
2e90089b 1098 if(fDebug)Printf("%s:%d iFoundBin %d",(char*)__FILE__,__LINE__,iFoundBin);
b360e4be 1099
d889ce29 1100 if(iFoundBin>0&&(fp1CalibRPYA)){
b360e4be 1101 meanXA = fp1CalibRPXA->GetBinContent(iFoundBin);
1102 meanYA = fp1CalibRPYA->GetBinContent(iFoundBin);
1103 meanXC = fp1CalibRPXC->GetBinContent(iFoundBin);
1104 meanYC = fp1CalibRPYC->GetBinContent(iFoundBin);
1105 }
1106
2e90089b 1107 if(fDebug)Printf("%s:%d iFoundBin %1.3E %1.3E %1.3E %1.3E",(char*)__FILE__,__LINE__,meanXA,meanYA,meanXC,meanYC);
b360e4be 1108
2eded560 1109 for (int i=0; i<64; i++) {
1110 Double_t mult = aodVZERO->GetMultiplicity(i);
1111 Double_t phi = v0phi[i];
1112 if (mult>0) {
1113 if (i<32) { //C-side
1114 Double_t wZNC= mult;
1115 numYZNC += sin(2.*phi)*wZNC;
1116 numXZNC += cos(2.*phi)*wZNC;
1117 sumZNC+=wZNC;
1118 }
1119 else if(i>31){ //A-side
1120 Double_t wZNA=mult;
1121 numYZNA += sin(2.*phi)*wZNA;
1122 numXZNA += cos(2.*phi)*wZNA;
1123 sumZNA+=wZNA;
1124 }
1125 }// mult>0
1126 }// 64 sectors
1127
1128 Double_t XC = numXZNC/sumZNC;
1129 Double_t YC = numYZNC/sumZNC;
40445651 1130
2eded560 1131 Double_t XA = numXZNA/sumZNA;
1132 Double_t YA = numYZNA/sumZNA;
40445651 1133
2eded560 1134
1135 fPsiVZEROA = 0.5*TMath::ATan2(YA-meanYA, XA-meanXA);
1136 if(fPsiVZEROA>TMath::Pi()){fPsiVZEROA-=TMath::Pi();}
1137 if(fPsiVZEROA<0){fPsiVZEROA+=TMath::Pi();}
1138
1139 fPsiVZEROC = 0.5*TMath::ATan2(YC-meanYC, XA-meanXC);
1140 if(fPsiVZEROC>TMath::Pi()){fPsiVZEROC-=TMath::Pi();}
1141 if(fPsiVZEROC<0){fPsiVZEROC+=TMath::Pi();}
40445651 1142
b360e4be 1143 fh2XYA->Fill(XA-meanXA,YA-meanYA); // control
2eded560 1144 fp1RPXA->Fill(aod->GetRunNumber(),XA);
1145 fp1RPYA->Fill(aod->GetRunNumber(),YA);
b360e4be 1146 fh2XYC->Fill(XC-meanXC,YC-meanYC); // control
2eded560 1147 fp1RPXC->Fill(aod->GetRunNumber(),XC);
1148 fp1RPYC->Fill(aod->GetRunNumber(),YC);
40445651 1149 return kTRUE;
40445651 1150
40445651 1151}
1152
40445651 1153Int_t AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1154 Int_t iCount = 0;
1155 AliAODEvent *aod = 0;
1156 if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1157 else aod = AODEvent();
1158 if(!aod){
1159 return iCount;
1160 }
1161 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1162 AliAODTrack *tr = aod->GetTrack(it);
1163 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1164 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1165 if(tr->Pt()<fMinTrackPt)continue;
1166 list->Add(tr);
1167 iCount++;
1168 }
1169 list->Sort();
1170 return iCount;
1171
1172}
1173