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