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