PWG4/JetTasks/AliAnalysisTaskJetServices.cxx

   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>
  25 #include <TString.h>
  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>
  37 #include <TRandom3.h>
  38 #include <TClonesArray.h>
  39 #include "TDatabasePDG.h"
  40
  41 #include "AliAnalysisTaskJetServices.h"
  42 #include "AliCentrality.h"
  43 #include "AliAnalysisDataContainer.h"
  44 #include "AliAnalysisDataSlot.h"
  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"
  64 #include "AliPhysicsSelection.h"
  65 #include "AliTriggerAnalysis.h"
  66
  67 #include "AliAnalysisHelperJetTasks.h"
  68
  69 ClassImp(AliAnalysisTaskJetServices)
  70
  71 AliAODHeader*  AliAnalysisTaskJetServices::fgAODHeader = NULL;
  72 TClonesArray*   AliAnalysisTaskJetServices::fgAODVertices = NULL;
  73
  74 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices():
  75   AliAnalysisTaskSE(),
  76   fUseAODInput(kFALSE),
  77   fUsePhysicsSelection(kFALSE),
  78   fMC(kFALSE),
  79   fFilterAODCollisions(kFALSE),
  80   fPhysicsSelectionFlag(AliVEvent::kMB),
  81   fSelectionInfoESD(0),
  82   fEventCutInfoESD(0),
  83   fFilterMask(0),
  84   fRPSubeventMethod(0),
  85   fCollisionType(kPbPb),
  86   fAvgTrials(1),
  87   fVtxXMean(0),
  88   fVtxYMean(0),
  89   fVtxZMean(0),
  90   fVtxRCut(1.),
  91   fVtxZCut(8.),
  92   fPtMinCosmic(5.),
  93   fRIsolMinCosmic(3.),
  94   fMaxCosmicAngle(0.01),
  95   fCentrality(101),
  96   fTrackRecEtaWindow(0.9),
  97   fMinTrackPt(0.15),
  98   fRPAngle(0),
  99   fRandomizer(0),
 100   fNonStdFile(""),
 101   fh1Xsec(0x0),
 102   fh1Trials(0x0),
 103   fh1PtHard(0x0),
 104   fh1PtHardTrials(0x0),
 105   fh1SelectionInfoESD(0x0),
 106   fh1EventCutInfoESD(0),
 107   fh1CentralityESD(0),
 108   fh1Centrality(0),
 109   fh1RP(0),
 110   fh2TriggerCount(0x0),
 111   fh2ESDTriggerCount(0x0),
 112   fh2TriggerVtx(0x0),
 113   fh2ESDTriggerVtx(0x0),
 114   fh2ESDTriggerRun(0x0),
 115   fh2VtxXY(0x0),
 116   fh1NCosmicsPerEvent(0x0),
 117   fh2RPSubevents(0x0),
 118   fh2RPCentrality(0x0),
 119   fh2RPDeltaRP(0x0),
 120   fh2RPQxQy(0x0),
 121   fh2RPCosDeltaRP(0x0),
 122   fh3PhiWeights(0x0),
 123   fh3RPPhiTracks(0x0),
 124   fTriggerAnalysis(0x0),
 125   fHistList(0x0)
 126 {
 127   fRunRange[0] = fRunRange[1] = 0;
 128   fFlatA[0] =   fFlatA[1] = 0;
 129   fFlatB[0] =   fFlatB[1] = 0;
 130   fDeltaQxy[0] =   fDeltaQxy[1] = 0;
 131
 132 }
 133
 134 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
 135   AliAnalysisTaskSE(name),
 136   fUseAODInput(kFALSE),
 137   fUsePhysicsSelection(kFALSE),
 138   fMC(kFALSE),
 139   fFilterAODCollisions(kFALSE),
 140   fPhysicsSelectionFlag(AliVEvent::kMB),
 141   fSelectionInfoESD(0),
 142   fEventCutInfoESD(0),
 143   fFilterMask(0),
 144   fRPSubeventMethod(0),
 145   fCollisionType(kPbPb),
 146   fAvgTrials(1),
 147   fVtxXMean(0),
 148   fVtxYMean(0),
 149   fVtxZMean(0),
 150   fVtxRCut(1.),
 151   fVtxZCut(8.),
 152   fPtMinCosmic(5.),
 153   fRIsolMinCosmic(3.),
 154   fMaxCosmicAngle(0.01),
 155   fCentrality(101),
 156   fTrackRecEtaWindow(0.9),
 157   fMinTrackPt(0.15),
 158   fRPAngle(0),
 159   fRandomizer(0),
 160   fNonStdFile(""),
 161   fh1Xsec(0x0),
 162   fh1Trials(0x0),
 163   fh1PtHard(0x0),
 164   fh1PtHardTrials(0x0),
 165   fh1SelectionInfoESD(0x0),
 166   fh1EventCutInfoESD(0),
 167   fh1CentralityESD(0),
 168   fh1Centrality(0),
 169   fh1RP(0),
 170   fh2TriggerCount(0x0),
 171   fh2ESDTriggerCount(0x0),
 172   fh2TriggerVtx(0x0),
 173   fh2ESDTriggerVtx(0x0),
 174   fh2ESDTriggerRun(0x0),
 175   fh2VtxXY(0x0),
 176   fh1NCosmicsPerEvent(0x0),
 177   fh2RPSubevents(0x0),
 178   fh2RPCentrality(0x0),
 179   fh2RPDeltaRP(0x0),
 180   fh2RPQxQy(0x0),
 181   fh2RPCosDeltaRP(0x0),
 182   fh3PhiWeights(0x0),
 183   fh3RPPhiTracks(0x0),
 184
 185   fTriggerAnalysis(0x0),
 186   fHistList(0x0)
 187 {
 188   fRunRange[0] = fRunRange[1] = 0;
 189   fFlatA[0] =   fFlatA[1] = 0;
 190   fFlatB[0] =   fFlatB[1] = 0;
 191   fDeltaQxy[0] =   fDeltaQxy[1] = 0;
 192   DefineOutput(1,TList::Class());
 193 }
 194
 195
 196
 197 Bool_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
 228 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 229 {
 230
 231   //
 232   // Create the output container
 233   //
 234
 235   fRandomizer = new TRandom3(0);
 236   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
 237
 238   OpenFile(1);
 239   if(!fHistList)fHistList = new TList();
 240   fHistList->SetOwner();
 241   PostData(1, fHistList); // post data in any case once
 242
 243   Bool_t oldStatus = TH1::AddDirectoryStatus();
 244   TH1::AddDirectory(kFALSE);
 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
 264
 265   fh1CentralityESD = new TH1F("fh1CentralityESD","cent",103,-1,102);
 266   fHistList->Add(fh1CentralityESD);
 267
 268   fh1Centrality = new TH1F("fh1Centrality","cent",103,-1,102);
 269   fHistList->Add(fh1Centrality);
 270
 271   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
 272   fHistList->Add(fh1RP);
 273
 274   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
 275   fHistList->Add(fh2TriggerCount);
 276
 277   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",6,-0.5,5.5,kConstraints,-0.5,kConstraints-0.5);
 278   fHistList->Add(fh2ESDTriggerCount);
 279   const Int_t nBins = 6*kConstraints;
 280   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
 281   fHistList->Add(fh2TriggerVtx);
 282
 283   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.);
 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);
 291   fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
 292                                  AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
 293   fHistList->Add(fh1SelectionInfoESD);
 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
 299   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
 300   // 3 triggers BB BE/EB EE
 301
 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);
 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);
 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
 318   fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
 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
 340   fHistList->Add(fh1NCosmicsPerEvent),
 341
 342
 343   TH1::AddDirectory(oldStatus);
 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());
 350      if (fDebug > 1) AliInfo("Replicating primary vertices");
 351      fgAODVertices = new TClonesArray("AliAODVertex",3);
 352      fgAODVertices->SetName("vertices");
 353      AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
 354   }
 355 }
 356
 357 void AliAnalysisTaskJetServices::Init()
 358 {
 359   //
 360   // Initialization
 361   //
 362   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
 363
 364 }
 365
 366 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 367 {
 368
 369   //
 370   // Execute analysis for current event
 371   //
 372
 373   AliAODEvent *aod = 0;
 374   AliESDEvent *esd = 0;
 375
 376
 377
 378   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
 379   AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
 380   AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,0); // set slection to false
 381   fSelectionInfoESD = 0; // reset
 382   fEventCutInfoESD = 0; // reset
 383   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
 384
 385
 386   static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 387
 388
 389
 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   }
 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   }
 413   fSelectionInfoESD |= kNoEventCut;
 414   fEventCutInfoESD |= kNoEventCut;
 415
 416   Bool_t esdVtxValid = false;
 417   Bool_t esdVtxIn = false;
 418   Bool_t aodVtxValid = false;
 419   Bool_t aodVtxIn = false;
 420
 421   if(esd){
 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;
 438   }
 439
 440   // Apply additional constraints
 441   Bool_t esdEventSelected = IsEventSelected(esd);
 442   Bool_t esdEventPileUp = IsEventPileUp(esd);
 443   Bool_t esdEventCosmic = IsEventCosmic(esd);
 444
 445   Bool_t aodEventSelected = IsEventSelected(aod);
 446
 447   Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
 448
 449
 450   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
 451   fh1EventCutInfoESD->Fill(fEventCutInfoESD);
 452
 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   }
 463   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD);
 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();
 469       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 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();
 474     }
 475   }
 476
 477   // loop over all possible triggers for esd
 478
 479   Float_t cent = 0;
 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   }
 485   fCentrality = cent;
 486   fRPAngle = 0;
 487
 488   if(esd){
 489     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 490     esdVtxValid = IsVertexValid(vtxESD);
 491     esdVtxIn = IsVertexIn(vtxESD);
 492     if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
 493       if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
 494       if(cent<=80&&esdVtxIn){
 495         aodH->SetFillAOD(kTRUE);
 496         aodH->SetFillExtension(kTRUE);
 497       }
 498     }
 499
 500
 501     Float_t zvtx = vtxESD->GetZ();
 502     Int_t  iCl = GetEventClass(esd);
 503     AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
 504     Bool_t cand = physicsSelection;
 505
 506     if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
 507     fh2ESDTriggerCount->Fill(0.,kAllTriggered);
 508     fh2ESDTriggerCount->Fill(iCl,kAllTriggered);
 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);
 523       }
 524       if(cand){
 525         fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
 526         fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
 527         fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
 528       }
 529     }
 530
 531     if(cand&&esdVtxIn&&iCl<5){
 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
 539       if(esd->GetCentrality()){
 540         Float_t tmpCent = 100;
 541         tmpCent = esd->GetCentrality()->GetCentralityPercentile("V0M");
 542         if(tmpCent<0)tmpCent = 101;
 543         fh1CentralityESD->Fill(tmpCent);
 544       }
 545     }
 546   }
 547
 548
 549
 550   if(aod){
 551     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 552     aodVtxValid = IsVertexValid(vtxAOD);
 553     aodVtxIn = IsVertexIn(vtxAOD);
 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());
 559     fh2TriggerCount->Fill(0.,kAllTriggered);
 560     fh2TriggerCount->Fill(iCl,kAllTriggered);
 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);
 574       }
 575       if(cand){
 576         fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
 577         fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
 578         fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
 579       }
 580     }
 581
 582     if(cand&&aodVtxIn&&iCl<5){
 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);
 589       fh1Centrality->Fill(cent);
 590       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 591       if(aodH&&cand&&fFilterAODCollisions&&!esd){
 592         if(fCentrality<=80&&aodVtxIn){
 593           aodH->SetFillAOD(kTRUE);
 594           aodH->SetFillExtension(kTRUE);
 595         }
 596       }
 597
 598       TList recTracks;
 599       GetListOfTracks(&recTracks);
 600       CalculateReactionPlaneAngle(&recTracks);
 601       AliAnalysisHelperJetTasks::ReactionPlane(kTRUE,fRPAngle); // set slection to false
 602       if(fUseAODInput&&fCentrality<=80){
 603         if(fFilterAODCollisions&&aod){
 604           aodH->SetFillAOD(kTRUE);
 605         }
 606       }
 607     }
 608   }
 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
 641   // replication of
 642   if(fNonStdFile.Length()&&aod){
 643     if (fgAODHeader){
 644       *fgAODHeader =  *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
 645       Double_t q[2] = {fRPAngle,fRPAngle};
 646       fgAODHeader->SetQTheta(q,2);
 647     }
 648     if(fgAODVertices){
 649       fgAODVertices->Delete();
 650       TClonesArray &vertices = *fgAODVertices;
 651       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 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());
 665       TString vtitle = vtxAOD->GetTitle();
 666       vtx->SetNContributors(vtxAOD->GetNContributors());
 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       }
 689     }
 690   }
 691
 692   PostData(1, fHistList);
 693 }
 694
 695 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
 696   if(!esd)return kFALSE;
 697   const AliESDVertex *vtx = esd->GetPrimaryVertex();
 698   return IsVertexIn(vtx); // vertex in calls vertex valid
 699 }
 700
 701 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
 702   if(fgAODVertices){
 703     fgAODVertices->Delete();
 704     delete fgAODVertices;
 705   }
 706   delete fRandomizer;
 707   if(fgAODHeader)delete fgAODHeader;
 708 }
 709
 710
 711 Bool_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
 717 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
 718
 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
 723   // Tracks    PrimaryVertex         VertexerTracksWithConstraint
 724   // TPC       TPCVertex             VertexerTracksNoConstraint
 725   // TPC       TPCVertex             VertexerTracksWithConstraint
 726   // SPD       SPDVertex             vertexer: 3D
 727   // SPD       SPDVertex             vertexer: Z
 728
 729   Int_t nCont = vtx->GetNContributors();
 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;
 741
 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;
 750
 751
 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
 761 Bool_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
 770
 771   if(fDebug){
 772     Printf(" n contrib %d",vtx->GetNContributors());
 773     vtx->Print();
 774   }
 775
 776   //  if(vtx->GetNContributors()<3)return kFALSE;
 777   // do not want tpc only primary vertex
 778   TString vtxName(vtx->GetName());
 779   if(vtxName.Contains("TPCVertex"))return kFALSE;
 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
 792 Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
 793
 794   if(!IsVertexValid(vtx))return kFALSE;
 795
 796   Float_t zvtx = vtx->GetZ();
 797   Float_t yvtx = vtx->GetY();
 798   Float_t xvtx = vtx->GetX();
 799
 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
 819 Bool_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();
 826
 827   xvtx -= fVtxXMean;
 828   yvtx -= fVtxYMean;
 829   zvtx -= fVtxZMean;
 830
 831   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 832
 833   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 834   return vertexIn;
 835 }
 836
 837 Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
 838   if(!esd)return kFALSE;
 839   return esd->IsPileupFromSPD();
 840 }
 841
 842 Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
 843   if(!esd)return kFALSE;
 844   // add track cuts for which we look for cosmics...
 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;
 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;
 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;
 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;
 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) {
 881         nCosmicCandidates+=1;
 882         isCosmic = kTRUE;
 883       }
 884
 885     }
 886   }
 887
 888   fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
 889
 890   return isCosmic;
 891 }
 892
 893
 894 Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
 895
 896   Float_t cent = 0;
 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;
 907   }
 908   return 1;
 909 }
 910
 911
 912 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
 913
 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   }
 922   return 1;
 923
 924 }
 925
 926
 927 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 928 {
 929   // Terminate analysis
 930   //
 931 }
 932
 933 Bool_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
1054 Double_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
1061 Int_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