PWG4/JetTasks/AliAnalysisTaskJetServices.cxx

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