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