]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetServices.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added switch for event slection in cluster task, updated plot notand JetSpectrum...
[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 }
90a006c2 576
c08c3ad2 577 if(cand&&aodVtxIn&&iCl<5){
ffab794c 578 fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
579 fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
580 fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
581 fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
582 fh2TriggerCount->Fill(iCl,kSelected);
583 fh2TriggerCount->Fill(0.,kSelected);
742ee86c 584 fh1Centrality->Fill(cent);
585 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
90a006c2 586 if(aodH&&cand&&fFilterAODCollisions&&!esd){
587 if(cent<=80&&aodVtxIn){
588 aodH->SetFillAOD(kTRUE);
589 aodH->SetFillExtension(kTRUE);
590 }
591 }
592
40445651 593 TList recTracks;
594 GetListOfTracks(&recTracks);
595 CalculateReactionPlaneAngle(&recTracks);
596 AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
742ee86c 597 if(fUseAODInput&&cent<=80){
221d8aa9 598 if(fFilterAODCollisions&&aod){
742ee86c 599 aodH->SetFillAOD(kTRUE);
221d8aa9 600 }
601 }
221d8aa9 602 }
c0997643 603 }
bf7b8731 604
605 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
606
607
608 Double_t ptHard = 0;
609 Double_t nTrials = 1; // Trials for MC trigger
610
611 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
612 AliMCEvent* mcEvent = MCEvent();
613 // AliStack *pStack = 0;
614 if(mcEvent){
615 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
616 if(pythiaGenHeader){
617 nTrials = pythiaGenHeader->Trials();
618 ptHard = pythiaGenHeader->GetPtHard();
619 int iProcessType = pythiaGenHeader->ProcessType();
620 // 11 f+f -> f+f
621 // 12 f+barf -> f+barf
622 // 13 f+barf -> g+g
623 // 28 f+g -> f+g
624 // 53 g+g -> f+barf
625 // 68 g+g -> g+g
626 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
627 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
628 fh1PtHard->Fill(ptHard);
629 fh1PtHardTrials->Fill(ptHard,nTrials);
630
631 }// if pythia gen header
632 }
633
634 // trigger selection
635
f4132e7d 636 // replication of
637 if(fNonStdFile.Length()&&aod){
638 if (fgAODHeader){
639 *fgAODHeader = *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
40445651 640 Double_t q[2] = {fRPAngle,fRPAngle};
641 fgAODHeader->SetQTheta(q,2);
f4132e7d 642 }
80a790fd 643 if(fgAODVertices){
80a790fd 644 fgAODVertices->Delete();
645 TClonesArray &vertices = *fgAODVertices;
80a790fd 646 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
070f06d2 647 // we only use some basic information,
648
649
650 Double_t pos[3];
651 Double_t covVtx[6];
652
653 vtxAOD->GetXYZ(pos); // position
654 vtxAOD->GetCovMatrix(covVtx); //covariance matrix
655 Int_t jVertices = 0;
656 AliAODVertex * vtx = new(vertices[jVertices++])
657 AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
658 vtx->SetName(vtxAOD->GetName());
659 vtx->SetTitle(vtxAOD->GetTitle());
070f06d2 660 TString vtitle = vtxAOD->GetTitle();
d8f21f85 661 vtx->SetNContributors(vtxAOD->GetNContributors());
070f06d2 662
663 // Add SPD "main" vertex
664 const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
665 vtxS->GetXYZ(pos); // position
666 vtxS->GetCovMatrix(covVtx); //covariance matrix
667 AliAODVertex * mVSPD = new(vertices[jVertices++])
668 AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
669 mVSPD->SetName(vtxS->GetName());
670 mVSPD->SetTitle(vtxS->GetTitle());
671 mVSPD->SetNContributors(vtxS->GetNContributors());
672
673 // Add tpc only vertex
674 if(esd){
675 const AliESDVertex *vtxT = esd->GetPrimaryVertexTPC();
676 vtxT->GetXYZ(pos); // position
677 vtxT->GetCovMatrix(covVtx); //covariance matrix
678 AliAODVertex * mVTPC = new(vertices[jVertices++])
679 AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
680 mVTPC->SetName(vtxT->GetName());
681 mVTPC->SetTitle(vtxT->GetTitle());
682 mVTPC->SetNContributors(vtxT->GetNContributors());
683 }
80a790fd 684 }
f4132e7d 685 }
686
bf7b8731 687 PostData(1, fHistList);
80a790fd 688}
bf7b8731 689
6d597ca2 690Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
c0997643 691 if(!esd)return kFALSE;
692 const AliESDVertex *vtx = esd->GetPrimaryVertex();
6d597ca2 693 return IsVertexIn(vtx); // vertex in calls vertex valid
694}
695
80a790fd 696AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
697 if(fgAODVertices){
698 fgAODVertices->Delete();
699 delete fgAODVertices;
700 }
40445651 701 delete fRandomizer;
070f06d2 702 if(fgAODHeader)delete fgAODHeader;
80a790fd 703}
704
6d597ca2 705
706Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
707 if(!aod)return kFALSE;
708 const AliAODVertex *vtx = aod->GetPrimaryVertex();
709 return IsVertexIn(vtx); // VertexIn calls VertexValid
710}
711
712Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
713
c0997643 714 // We can check the type of the vertex by its name and title
715 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
716 // vtx name title
717 // Tracks PrimaryVertex VertexerTracksNoConstraint
6d597ca2 718 // Tracks PrimaryVertex VertexerTracksWithConstraint
c0997643 719 // TPC TPCVertex VertexerTracksNoConstraint
6d597ca2 720 // TPC TPCVertex VertexerTracksWithConstraint
c0997643 721 // SPD SPDVertex vertexer: 3D
722 // SPD SPDVertex vertexer: Z
723
6d597ca2 724 Int_t nCont = vtx->GetNContributors();
6d597ca2 725 if(nCont>=1){
726 fEventCutInfoESD |= kContributorsCut1;
727 if(nCont>=2){
728 fEventCutInfoESD |= kContributorsCut2;
729 if(nCont>=3){
730 fEventCutInfoESD |= kContributorsCut3;
731 }
732 }
733 }
734
735 if(nCont<3)return kFALSE;
c0997643 736
6d597ca2 737 // do not want tpc only primary vertex
738 TString vtxName(vtx->GetName());
739 if(vtxName.Contains("TPCVertex")){
740 fEventCutInfoESD |= kVertexTPC;
741 return kFALSE;
742 }
743 if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
744 if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
c0997643 745
c0997643 746
6d597ca2 747 TString vtxTitle(vtx->GetTitle());
748 if(vtxTitle.Contains("vertexer: Z")){
749 if(vtx->GetDispersion()>0.02)return kFALSE;
750 }
751 fEventCutInfoESD |= kSPDDispersionCut;
752 return kTRUE;
753}
754
755
756Bool_t AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
757
758 // We can check the type of the vertex by its name and title
759 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
760 // vtx name title
761 // Tracks PrimaryVertex VertexerTracksNoConstraint
762 // TPC TPCVertex VertexerTracksNoConstraint
763 // SPD SPDVertex vertexer: 3D
764 // SPD SPDVertex vertexer: Z
d8f21f85 765
766 if(fDebug){
767 Printf(" n contrib %d",vtx->GetNContributors());
768 vtx->Print();
769 }
6d597ca2 770
d8f21f85 771 // if(vtx->GetNContributors()<3)return kFALSE;
c0997643 772 // do not want tpc only primary vertex
773 TString vtxName(vtx->GetName());
774 if(vtxName.Contains("TPCVertex"))return kFALSE;
6d597ca2 775
776 // no dispersion yet...
777 /*
778 TString vtxTitle(vtx->GetTitle());
779 if(vtxTitle.Contains("vertexer: Z")){
780 if(vtx->GetDispersion()>0.02)return kFALSE;
781 }
782 */
783 return kTRUE;
784}
785
786
787Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
788
789 if(!IsVertexValid(vtx))return kFALSE;
790
c0997643 791 Float_t zvtx = vtx->GetZ();
792 Float_t yvtx = vtx->GetY();
793 Float_t xvtx = vtx->GetX();
794
6d597ca2 795 xvtx -= fVtxXMean;
796 yvtx -= fVtxYMean;
797 zvtx -= fVtxZMean;
798
799
800
801 if(TMath::Abs(zvtx)>fVtxZCut){
802 return kFALSE;
803 }
804 fEventCutInfoESD |= kVertexZCut;
805 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
806 if(r2>(fVtxRCut*fVtxRCut)){
807 return kFALSE;
808 }
809 fEventCutInfoESD |= kVertexRCut;
810 return kTRUE;
811}
812
813
814Bool_t AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
815
816 if(!IsVertexValid(vtx))return kFALSE;
817
818 Float_t zvtx = vtx->GetZ();
819 Float_t yvtx = vtx->GetY();
820 Float_t xvtx = vtx->GetX();
c0997643 821
822 xvtx -= fVtxXMean;
823 yvtx -= fVtxYMean;
824 zvtx -= fVtxZMean;
825
826 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
827
6d597ca2 828 Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
829 return vertexIn;
c0997643 830}
831
6d597ca2 832Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
f2dd0695 833 if(!esd)return kFALSE;
834 return esd->IsPileupFromSPD();
835}
836
6d597ca2 837Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
f2dd0695 838 if(!esd)return kFALSE;
839 // add track cuts for which we look for cosmics...
df7848fc 840
841 Bool_t isCosmic = kFALSE;
842 Int_t nTracks = esd->GetNumberOfTracks();
843 Int_t nCosmicCandidates = 0;
844
845 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
846 AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
847 if (!track1) continue;
75946d5d 848 UInt_t status1 = track1->GetStatus();
849 //If track is ITS stand alone track, skip the track
850 if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
df7848fc 851 if(track1->Pt()<fPtMinCosmic) continue;
852 //Start 2nd track loop to look for correlations
853 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
854 AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
855 if(!track2) continue;
75946d5d 856 UInt_t status2 = track2->GetStatus();
857 //If track is ITS stand alone track, skip the track
858 if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
df7848fc 859 if(track2->Pt()<fPtMinCosmic) continue;
860 //Check if back-to-back
861 Double_t mom1[3],mom2[3];
862 track1->GetPxPyPz(mom1);
863 track2->GetPxPyPz(mom2);
864 TVector3 momv1(mom1[0],mom1[1],mom1[2]);
865 TVector3 momv2(mom2[0],mom2[1],mom2[2]);
866 Float_t theta = (float)(momv1.Phi()-momv2.Phi());
867 if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
868
869 Float_t deltaPhi = track1->Phi()-track2->Phi();
870 if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
871
872 Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
873 if(rIsol<fRIsolMinCosmic) continue;
874
875 if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
75946d5d 876 nCosmicCandidates+=1;
df7848fc 877 isCosmic = kTRUE;
878 }
879
880 }
881 }
882
883 fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
884
885 return isCosmic;
f2dd0695 886}
887
888
f4132e7d 889Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
890
891 Float_t cent = 999;
892 if(esd->GetCentrality()){
893 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
894 }
742ee86c 895 if(cent>80||cent<0)return 5;
f4132e7d 896 if(cent>50)return 4;
897 if(cent>30)return 3;
898 if(cent>10)return 2;
899 return 1;
bf7b8731 900
f4132e7d 901}
f2dd0695 902
ffab794c 903
904Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
905
906 Float_t cent = aod->GetHeader()->GetCentrality();
742ee86c 907 if(cent>80||cent<0)return 5;
ffab794c 908 if(cent>50)return 4;
909 if(cent>30)return 3;
910 if(cent>10)return 2;
911 return 1;
912
913}
914
915
bf7b8731 916void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
917{
918 // Terminate analysis
919 //
bf7b8731 920}
f4132e7d 921
40445651 922Bool_t AliAnalysisTaskJetServices::CalculateReactionPlaneAngle(const TList *trackList)
923{
924
925 if(!trackList)return kFALSE;
926 fRPAngle=0;
927
928 // need to get this info from elsewhere??
929
930 Double_t fPsiRP =0,fDeltaPsiRP = 0;
931
932
933
934 TVector2 mQ,mQ1,mQ2;
935 Float_t mQx= fDeltaQxy[0], mQy=fDeltaQxy[1];
936
937 Float_t mQx1=fDeltaQxy[0], mQy1=fDeltaQxy[1];
938 Float_t mQx2=fDeltaQxy[0], mQy2=fDeltaQxy[1];
939
940 AliVParticle *track=0x0;
941 Int_t count[3]={0,0,0};
942
943
944 for (Int_t iter=0;iter<trackList->GetEntries();iter++){
945
946 track=(AliVParticle*)trackList->At(iter);
947
948 //cuts already applied before
949 // Comment DCA not correctly implemented yet for AOD tracks
950
951 Double_t momentum;
952 if(track->Charge()>0){momentum=track->Pt();}
953 else{momentum=-track->Pt();}
954
955
956
957 // For Weighting
958 fh3RPPhiTracks->Fill(fCentrality,momentum,track->Phi());
959 count[0]++;
960
961 Double_t phiweight=GetPhiWeight(track->Phi(),momentum);
962 // Double_t phiweight=1;
963 Double_t weight=2;
964 if(track->Pt()<2){weight=track->Pt();}
965
966
967 mQx += (cos(2*track->Phi()))*weight*phiweight;
968 mQy += (sin(2*track->Phi()))*weight*phiweight;
969
970 // Make random Subevents
971
972 if(fRPSubeventMethod==0){
973 if(fRandomizer->Binomial(1,0.5)){
974 mQx1 += (cos(2*track->Phi()))*weight*phiweight;
975 mQy1 += (sin(2*track->Phi()))*weight*phiweight;
976 count[1]++;}
977 else{
978 mQx2 += (cos(2*track->Phi()))*weight*phiweight;
979 mQy2 += (sin(2*track->Phi()))*weight*phiweight;
980 count[2]++;}
981 }
982 else if(fRPSubeventMethod==1){
983 // Make eta dependent subevents
984 if(track->Eta()>0){
985 mQx1 += (cos(2*track->Phi()))*weight*phiweight;
986 mQy1 += (sin(2*track->Phi()))*weight*phiweight;
987 count[1]++;}
988 else{
989 mQx2 += (cos(2*track->Phi()))*weight*phiweight;
990 mQy2 += (sin(2*track->Phi()))*weight*phiweight;
991 count[2]++;}
992 }
993
994 }
995
996
997
998 //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
999 if(count[0]==0||count[1]==0||count[2]==0){
1000 return kFALSE;
1001 }
1002
1003 mQ.Set(mQx,mQy);
1004 mQ1.Set(mQx1,mQy1);
1005 mQ2.Set(mQx2,mQy2);
1006
1007 // cout<<"MQ"<<mQx<<" " <<mQy<<" psi"<<endl;
1008
1009 fPsiRP=mQ.Phi()/2;
1010
1011 //Correction
1012 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);
1013
1014 Double_t fPsiRP1=mQ1.Phi()/2;
1015 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);
1016 Double_t fPsiRP2=mQ2.Phi()/2;
1017 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);
1018 fDeltaPsiRP=fPsiRP1-fPsiRP2;
1019
1020 if(fPsiRP>TMath::Pi()){fPsiRP-=TMath::Pi();}
1021 if(fPsiRP<0){fPsiRP+=TMath::Pi();}
1022
1023 // reactionplaneangle + Pi() is the same angle
1024 if(TMath::Abs(fDeltaPsiRP)>TMath::Pi()/2){
1025 if(fDeltaPsiRP>0)fDeltaPsiRP-=TMath::Pi();
1026 else fDeltaPsiRP+=TMath::Pi();
1027 }
1028
1029 Double_t cos2deltaRP=TMath::Cos(2*fDeltaPsiRP);
1030
1031 // FillHistograms
1032 fh2RPSubevents->Fill(fPsiRP1,fPsiRP2);
1033 fh1RP->Fill(fPsiRP);
1034 fh2RPCentrality->Fill(fCentrality,fPsiRP);
1035 fh2RPDeltaRP->Fill(fDeltaPsiRP,fCentrality);
1036 fh2RPQxQy->Fill(mQx,mQy);
1037 fh2RPCosDeltaRP->Fill(fCentrality,cos2deltaRP);
1038
1039 fRPAngle=fPsiRP;
1040 return kTRUE;
1041}
1042
1043Double_t AliAnalysisTaskJetServices::GetPhiWeight(Double_t phi,Double_t signedpt){
1044 if(!fh3PhiWeights)return 1;
1045 else return fh3PhiWeights->GetBinContent(fh3PhiWeights->GetXaxis()->FindBin(fCentrality),fh3PhiWeights->GetYaxis()->FindBin(signedpt),fh3PhiWeights->GetZaxis()->FindBin(phi));
1046}
1047
1048 //________________________________________________________________________
1049
1050Int_t AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
1051 Int_t iCount = 0;
1052 AliAODEvent *aod = 0;
1053 if(fUseAODInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1054 else aod = AODEvent();
1055 if(!aod){
1056 return iCount;
1057 }
1058 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1059 AliAODTrack *tr = aod->GetTrack(it);
1060 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1061 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1062 if(tr->Pt()<fMinTrackPt)continue;
1063 list->Add(tr);
1064 iCount++;
1065 }
1066 list->Sort();
1067 return iCount;
1068
1069}
1070