PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx

   1 // **************************************
   2 //  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 <TRandom3.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 <TClonesArray.h>
  38 #include <TRefArray.h>
  39 #include  "TDatabasePDG.h"
  40
  41 #include "AliAnalysisTaskJetSpectrum2.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliJetFinder.h"
  44 #include "AliJetHeader.h"
  45 #include "AliJetReader.h"
  46 #include "AliJetReaderHeader.h"
  47 #include "AliUA1JetHeaderV1.h"
  48 #include "AliESDEvent.h"
  49 #include "AliAODEvent.h"
  50 #include "AliAODHandler.h"
  51 #include "AliAODTrack.h"
  52 #include "AliAODJet.h"
  53 #include "AliAODJetEventBackground.h"
  54 #include "AliAODMCParticle.h"
  55 #include "AliMCEventHandler.h"
  56 #include "AliMCEvent.h"
  57 #include "AliStack.h"
  58 #include "AliGenPythiaEventHeader.h"
  59 #include "AliJetKineReaderHeader.h"
  60 #include "AliGenCocktailEventHeader.h"
  61 #include "AliInputEventHandler.h"
  62
  63
  64 #include "AliAnalysisHelperJetTasks.h"
  65
  66 ClassImp(AliAnalysisTaskJetSpectrum2)
  67
  68 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
  69   AliAnalysisTaskSE(),
  70   fJetHeaderRec(0x0),
  71   fJetHeaderGen(0x0),
  72   fAODIn(0x0),
  73   fAODOut(0x0),
  74   fAODExtension(0x0),
  75   fhnCorrelation(0x0),
  76   fhnCorrelationPhiZRec(0x0),
  77   f1PtScale(0x0),
  78   fBranchRec("jets"),
  79   fBranchGen(""),
  80   fBranchBkgRec(""),
  81   fBranchBkgGen(""),
  82   fNonStdFile(""),
  83   fRandomizer(0x0),
  84   fUseAODJetInput(kFALSE),
  85   fUseAODTrackInput(kFALSE),
  86   fUseAODMCInput(kFALSE),
  87   fUseGlobalSelection(kFALSE),
  88   fUseExternalWeightOnly(kFALSE),
  89   fLimitGenJetEta(kFALSE),
  90   fNMatchJets(5),
  91   fNRPBins(3),
  92   fFilterMask(0),
  93   fEventSelectionMask(0),
  94   fAnalysisType(0),
  95   fTrackTypeRec(kTrackUndef),
  96   fTrackTypeGen(kTrackUndef),
  97   fEventClass(0),
  98   fRPSubeventMethod(0),
  99   fAvgTrials(1),
 100   fExternalWeight(1),
 101   fJetRecEtaWindow(0.5),
 102   fTrackRecEtaWindow(0.5),
 103   fMinJetPt(0),
 104   fMinTrackPt(0.15),
 105   fDeltaPhiWindow(90./180.*TMath::Pi()),
 106   fCentrality(100),
 107   fRPAngle(0),
 108   fMultRec(0),
 109   fMultGen(0),
 110   fh1Xsec(0x0),
 111   fh1Trials(0x0),
 112   fh1PtHard(0x0),
 113   fh1PtHardNoW(0x0),
 114   fh1PtHardTrials(0x0),
 115   fh1ZVtx(0x0),
 116   fh1RP(0x0),
 117   fh1Centrality(0x0),
 118   fh1TmpRho(0x0),
 119   fh2MultRec(0x0),
 120   fh2MultGen(0x0),
 121   fh2RPSubevents(0x0),
 122   fh2RPCentrality(0x0),
 123   fh2RPDeltaRP(0x0),
 124   fh2RPQxQy(0x0),
 125   fh2RPCosDeltaRP(0x0),
 126   fh2PtFGen(0x0),
 127   fh2RelPtFGen(0x0),
 128   fh3RPPhiTracks(0x0),
 129   fHistList(0x0)
 130 {
 131   for(int i = 0;i < kMaxStep*2;++i){
 132     fhnJetContainer[i] = 0;
 133   }
 134
 135   for(int ij = 0;ij <kJetTypes;++ij){
 136     fFlagJetType[ij] = 1; // default = on
 137     fh1NJets[ij] = 0;
 138     fh1SumPtTrack[ij] = 0;
 139     fh1PtJetsIn[ij] = 0;
 140     fh1PtTracksIn[ij] = 0;
 141     fh1PtTracksInLow[ij] = 0;
 142     fh1PtTracksLeadingIn[ij] = 0;
 143     fh2NJetsPt[ij]  = 0;
 144     fh2NTracksPt[ij]  = 0;
 145     fh3MultTrackPtRP[ij] = 0;
 146     fh2LeadingTrackPtTrackPhi[ij] = 0;
 147     for(int i = 0;i <= kMaxJets;++i){
 148       fh3MultPtRP[ij][i] = 0;
 149       fh2PhiPt[ij][i] = 0;
 150       fh2EtaPt[ij][i] = 0;
 151       fh2AreaPt[ij][i] = 0;
 152       fh2EtaArea[ij][i] = 0;
 153       fh2PhiEta[ij][i] = 0;
 154       fh2LTrackPtJetPt[ij][i] = 0;
 155       fh1PtIn[ij][i] = 0;
 156       fh2RhoPt[ij][i] = 0;
 157       fh2PsiPt[ij][i] = 0;
 158     }
 159
 160     fh1DijetMinv[ij] = 0;
 161     fh2DijetDeltaPhiPt[ij] = 0;
 162     fh2DijetAsymPt[ij] = 0;
 163     fh2DijetPt2vsPt1[ij] = 0;
 164     fh2DijetDifvsSum[ij] = 0;
 165
 166   }
 167 }
 168
 169 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 170   AliAnalysisTaskSE(name),
 171   fJetHeaderRec(0x0),
 172   fJetHeaderGen(0x0),
 173   fAODIn(0x0),
 174   fAODOut(0x0),
 175   fAODExtension(0x0),
 176   fhnCorrelation(0x0),
 177   fhnCorrelationPhiZRec(0x0),
 178   f1PtScale(0x0),
 179   fBranchRec("jets"),
 180   fBranchGen(""),
 181   fBranchBkgRec(""),
 182   fBranchBkgGen(""),
 183   fNonStdFile(""),
 184   fRandomizer(0x0),
 185   fUseAODJetInput(kFALSE),
 186   fUseAODTrackInput(kFALSE),
 187   fUseAODMCInput(kFALSE),
 188   fUseGlobalSelection(kFALSE),
 189   fUseExternalWeightOnly(kFALSE),
 190   fLimitGenJetEta(kFALSE),
 191   fNMatchJets(5),
 192   fNRPBins(3),
 193   fFilterMask(0),
 194   fEventSelectionMask(0),
 195   fAnalysisType(0),
 196   fTrackTypeRec(kTrackUndef),
 197   fTrackTypeGen(kTrackUndef),
 198   fEventClass(0),
 199   fRPSubeventMethod(0),
 200   fAvgTrials(1),
 201   fExternalWeight(1),
 202   fJetRecEtaWindow(0.5),
 203   fTrackRecEtaWindow(0.5),
 204   fMinJetPt(0),
 205   fMinTrackPt(0.15),
 206   fDeltaPhiWindow(90./180.*TMath::Pi()),
 207   fCentrality(100),
 208   fRPAngle(0),
 209   fMultRec(0),
 210   fMultGen(0),
 211   fh1Xsec(0x0),
 212   fh1Trials(0x0),
 213   fh1PtHard(0x0),
 214   fh1PtHardNoW(0x0),
 215   fh1PtHardTrials(0x0),
 216   fh1ZVtx(0x0),
 217   fh1RP(0x0),
 218   fh1Centrality(0x0),
 219   fh1TmpRho(0x0),
 220   fh2MultRec(0x0),
 221   fh2MultGen(0x0),
 222   fh2RPSubevents(0x0),
 223   fh2RPCentrality(0x0),
 224   fh2RPDeltaRP(0x0),
 225   fh2RPQxQy(0x0),
 226   fh2RPCosDeltaRP(0x0),
 227   fh2PtFGen(0x0),
 228   fh2RelPtFGen(0x0),
 229   fh3RPPhiTracks(0x0),
 230   fHistList(0x0)
 231 {
 232
 233   for(int i = 0;i < kMaxStep*2;++i){
 234     fhnJetContainer[i] = 0;
 235   }
 236
 237   for(int ij = 0;ij <kJetTypes;++ij){
 238     fFlagJetType[ij] = 1; // default = on
 239     fh1NJets[ij] = 0;
 240     fh1SumPtTrack[ij] = 0;
 241     fh1PtJetsIn[ij] = 0;
 242     fh1PtTracksIn[ij] = 0;
 243     fh1PtTracksInLow[ij] = 0;
 244     fh1PtTracksLeadingIn[ij] = 0;
 245     fh2NJetsPt[ij]  = 0;
 246     fh2NTracksPt[ij]  = 0;
 247     fh3MultTrackPtRP[ij] = 0;
 248     fh2LeadingTrackPtTrackPhi[ij] = 0;
 249     for(int i = 0;i <= kMaxJets;++i){
 250       fh3MultPtRP[ij][i] = 0;
 251       fh2PhiPt[ij][i] = 0;
 252       fh2EtaPt[ij][i] = 0;
 253       fh2AreaPt[ij][i] = 0;
 254       fh2EtaArea[ij][i] = 0;
 255       fh2PhiEta[ij][i] = 0;
 256       fh2LTrackPtJetPt[ij][i] = 0;
 257       fh1PtIn[ij][i] = 0;
 258       fh2RhoPt[ij][i] = 0;
 259       fh2PsiPt[ij][i] = 0;
 260     }
 261
 262     fh1DijetMinv[ij] = 0;
 263     fh2DijetDeltaPhiPt[ij] = 0;
 264     fh2DijetAsymPt[ij] = 0;
 265     fh2DijetPt2vsPt1[ij] = 0;
 266     fh2DijetDifvsSum[ij] = 0;
 267   }
 268
 269   DefineOutput(1, TList::Class());
 270 }
 271
 272
 273
 274 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 275 {
 276
 277
 278
 279   //
 280   // Implemented Notify() to read the cross sections
 281   // and number of trials from pyxsec.root
 282   //
 283
 284   // Fetch the aod also from the input in,
 285   // have todo it in notify
 286
 287
 288   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 289   //  assume that the AOD is in the general output...
 290   fAODOut  = AODEvent();
 291
 292   if(fNonStdFile.Length()!=0){
 293     // case that we have an AOD extension we need can fetch the jets from the extended output
 294     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 295     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 296     if(!fAODExtension){
 297       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 298     }
 299   }
 300
 301
 302   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 303   Float_t xsection = 0;
 304   Float_t ftrials  = 1;
 305
 306   fAvgTrials = 1;
 307   if(tree){
 308     TFile *curfile = tree->GetCurrentFile();
 309     if (!curfile) {
 310       Error("Notify","No current file");
 311       return kFALSE;
 312     }
 313     if(!fh1Xsec||!fh1Trials){
 314       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 315       return kFALSE;
 316     }
 317     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 318     fh1Xsec->Fill("<#sigma>",xsection);
 319     // construct a poor man average trials
 320     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 321     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 322   }
 323
 324   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 325
 326   return kTRUE;
 327 }
 328
 329 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 330 {
 331
 332
 333   // Connect the AOD
 334
 335   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 336   OpenFile(1);
 337   if(!fHistList)fHistList = new TList();
 338   PostData(1, fHistList); // post data in any case once
 339
 340   if(!fRandomizer)fRandomizer = new TRandom3(0);
 341
 342   fHistList->SetOwner(kTRUE);
 343   Bool_t oldStatus = TH1::AddDirectoryStatus();
 344   TH1::AddDirectory(kFALSE);
 345
 346   MakeJetContainer();
 347   fHistList->Add(fhnCorrelation);
 348   if(fhnCorrelationPhiZRec)fHistList->Add(fhnCorrelationPhiZRec);
 349   for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 350
 351   //
 352   //  Histogram
 353
 354   const Int_t nBinPt = 160;
 355   Double_t binLimitsPt[nBinPt+1];
 356   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 357     if(iPt == 0){
 358       binLimitsPt[iPt] = 0.0;
 359     }
 360     else {// 1.0
 361       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
 362     }
 363   }
 364   const Int_t nBinPhi = 90;
 365   Double_t binLimitsPhi[nBinPhi+1];
 366   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 367     if(iPhi==0){
 368       binLimitsPhi[iPhi] = 0;
 369     }
 370     else{
 371       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 372     }
 373   }
 374
 375   const Int_t nBinEta = 40;
 376   Double_t binLimitsEta[nBinEta+1];
 377   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 378     if(iEta==0){
 379       binLimitsEta[iEta] = -2.0;
 380     }
 381     else{
 382       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 383     }
 384   }
 385
 386
 387   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 388   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 389    fHistList->Add(fh1Xsec);
 390   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 391   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 392   fHistList->Add(fh1Trials);
 393   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 394   fHistList->Add(fh1PtHard);
 395   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 396   fHistList->Add(fh1PtHardNoW);
 397   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 398   fHistList->Add(fh1PtHardTrials);
 399
 400   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 401   fHistList->Add(fh1ZVtx);
 402
 403   fh1RP = new TH1F("fh1RP","RP;#Psi",480,-180,360);
 404   fHistList->Add(fh1RP);
 405
 406   fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",101,-0.5,100.5);
 407   fHistList->Add(fh1Centrality);
 408
 409   fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
 410   fHistList->Add(fh2MultRec);
 411   fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
 412   fHistList->Add(fh2MultGen);
 413
 414   fh2RPSubevents = new TH2F("fh2RPSubevents" ,"Reaction Plane Angle" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
 415   fHistList->Add( fh2RPSubevents);
 416
 417   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
 418   fHistList->Add(fh2RPCentrality);
 419
 420   fh2RPDeltaRP   = new TH2F("fh2DeltaRP" ,"Delta Reaction Plane Angle" , 100, -TMath::Pi()/2, TMath::Pi()/2,20,0.,100.0);
 421   fHistList->Add(fh2RPDeltaRP);
 422
 423   fh2RPQxQy      = new TH2F("fh2RPQxQy" ,"" , 100, -100,100,100,-100,100);
 424   fHistList->Add(fh2RPQxQy);
 425
 426   fh2RPCosDeltaRP = new TH2F("fh2RPCosDeltaRP" ,"" , 20, 0.001,100.001,100,-1,1);
 427   fHistList->Add(fh2RPCosDeltaRP);
 428
 429
 430   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 431   fHistList->Add(fh2PtFGen);
 432
 433   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
 434   fHistList->Add(fh2RelPtFGen);
 435
 436   fh3RPPhiTracks = new TH3F("fh3RPPhiTracks","Phi Tracks Pt Centrality", 10, 0.,100.,20,-5,5,180, 0, 2*TMath::Pi());
 437   fHistList->Add(fh3RPPhiTracks);
 438
 439   for(int ij = 0;ij <kJetTypes;++ij){
 440     TString cAdd = "";
 441     TString cJetBranch = "";
 442     if(ij==kJetRec){
 443       cAdd = "Rec";
 444       cJetBranch = fBranchRec.Data();
 445     }
 446     else if (ij==kJetGen){
 447       cAdd = "Gen";
 448       cJetBranch = fBranchGen.Data();
 449     }
 450     else if (ij==kJetRecFull){
 451       cAdd = "RecFull";
 452       cJetBranch = fBranchRec.Data();
 453     }
 454     else if (ij==kJetGenFull){
 455       cAdd = "GenFull";
 456       cJetBranch = fBranchGen.Data();
 457     }
 458
 459     if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
 460     if(!fFlagJetType[ij])continue;
 461
 462     fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
 463     fHistList->Add(fh1NJets[ij]);
 464
 465     fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 466     fHistList->Add(fh1PtJetsIn[ij]);
 467
 468     fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 469     fHistList->Add(fh1PtTracksIn[ij]);
 470
 471     fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,10.);
 472     fHistList->Add(fh1PtTracksInLow[ij]);
 473
 474     fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 475     fHistList->Add(fh1PtTracksLeadingIn[ij]);
 476
 477     fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),900,0.,900.);
 478     fHistList->Add(fh1SumPtTrack[ij]);
 479
 480     fh2NJetsPt[ij]  = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
 481     fHistList->Add(fh2NJetsPt[ij]);
 482
 483     fh2NTracksPt[ij]  = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,4000);
 484     fHistList->Add(fh2NTracksPt[ij]);
 485
 486     fh3MultTrackPtRP[ij]  = new TH3F(Form("fh3MultTrackPtRP%s",
 487                                           cAdd.Data()),Form("%s track p_T;# tracks;;p_{T} (GeV/c)",cAdd.Data()),400,0,4000,nBinPt,0,300,(Int_t)fNRPBins,-0.5,fNRPBins-0.5);
 488     fHistList->Add(fh3MultTrackPtRP[ij]);
 489
 490     fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
 491                                                nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
 492       fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
 493
 494       for(int i = 0;i <= kMaxJets;++i){
 495         fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
 496         fHistList->Add(fh1PtIn[ij][i]);
 497
 498         fh2RhoPt[ij][i] = new TH2F(Form("fh2RhoPt%s_j%d",cAdd.Data(),i),Form("jet shape rho for %s jets;r;p_{T};",cAdd.Data()),
 499                                    40,0.,2.,nBinPt,binLimitsPt);
 500         fHistList->Add(fh2RhoPt[ij][i]);
 501         if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 502         fh2PsiPt[ij][i] = new TH2F(Form("fh2PsiPt%s_j%d",cAdd.Data(),i),Form("jet shape #psi for %s jets;r;p_{T};",cAdd.Data()),
 503                                      40,0.,2.,nBinPt,binLimitsPt);
 504         fHistList->Add(fh2PsiPt[ij][i]);
 505         fh2PhiPt[ij][i] = new TH2F(Form("fh2PhiPt%s_j%d",cAdd.Data(),i),Form("pt vs phi %s;#phi;p_{T};",cAdd.Data()),
 506                                    nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 507
 508         fh3MultPtRP[ij][i]  = new TH3F(Form("fh3MultPtRP%s_j%d",cAdd.Data(),i),Form("%s jets p_T;# tracks;;p_{T} (GeV/c)",cAdd.Data()),400,0,4000,nBinPt,0,300,(Int_t)fNRPBins,-0.5,fNRPBins-0.5);
 509       fHistList->Add(fh3MultPtRP[ij][i]);
 510
 511
 512
 513         fHistList->Add(fh2PhiPt[ij][i]);
 514         fh2EtaPt[ij][i] = new TH2F(Form("fh2EtaPt%s_j%d",cAdd.Data(),i),Form("pt vs eta %s;#eta;p_{T};",cAdd.Data()),
 515                                    50,-1.,1.,nBinPt,binLimitsPt);
 516         fHistList->Add(fh2EtaPt[ij][i]);
 517         fh2AreaPt[ij][i] = new TH2F(Form("fh2AreaPt%s_j%d",cAdd.Data(),i),
 518                                     Form("pt vs area %s;area;p_{T};",
 519                                          cAdd.Data()),
 520                                     50,0.,1.,nBinPt,binLimitsPt);
 521         fHistList->Add(fh2AreaPt[ij][i]);
 522         fh2EtaArea[ij][i] = new TH2F(Form("fh2EtaArea%s_j%d",cAdd.Data(),i),
 523                                      Form("area vs eta %s;#eta;area;",
 524                                           cAdd.Data()),
 525                                      50,-1.,1.,50,0,1.);
 526         fHistList->Add(fh2EtaArea[ij][i]);
 527
 528         fh2PhiEta[ij][i] =  new TH2F(Form("fh2PhiEta%s_j%d",cAdd.Data(),i),Form("phi vs eta %s ;phi;p_{T};",cAdd.Data()),
 529                                      nBinPhi,binLimitsPhi,20,-1.,1.);
 530         fHistList->Add(fh2PhiEta[ij][i]);
 531
 532         fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
 533                                            Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
 534                                                 cAdd.Data()),
 535                                            200,0.,200.,nBinPt,binLimitsPt);
 536         fHistList->Add(fh2LTrackPtJetPt[ij][i]);
 537       }
 538
 539
 540       fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 541       fHistList->Add(fh1DijetMinv[ij]);
 542
 543       fh2DijetDeltaPhiPt[ij]       = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
 544       fHistList->Add(fh2DijetDeltaPhiPt[ij]);
 545
 546       fh2DijetAsymPt[ij]            = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
 547       fHistList->Add(fh2DijetAsymPt[ij]);
 548
 549       fh2DijetPt2vsPt1[ij]          = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
 550       fHistList->Add(fh2DijetPt2vsPt1[ij]);
 551       fh2DijetDifvsSum[ij]         = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
 552       fHistList->Add( fh2DijetDifvsSum[ij]);
 553     }
 554   // =========== Switch on Sumw2 for all histos ===========
 555   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 556     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 557     if (h1){
 558       h1->Sumw2();
 559       continue;
 560     }
 561     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 562     if(hn)hn->Sumw2();
 563   }
 564   TH1::AddDirectory(oldStatus);
 565 }
 566
 567 void AliAnalysisTaskJetSpectrum2::Init()
 568 {
 569   //
 570   // Initialization
 571   //
 572
 573   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 574
 575 }
 576
 577 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
 578
 579
 580   Bool_t selected = kTRUE;
 581
 582   if(fUseGlobalSelection&&fEventSelectionMask==0){
 583     selected = AliAnalysisHelperJetTasks::Selected();
 584   }
 585   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 586     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 587   }
 588
 589   if(fEventClass>0){
 590     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 591   }
 592
 593   if(!selected){
 594     // no selection by the service task, we continue
 595     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 596     PostData(1, fHistList);
 597     return;
 598   }
 599
 600
 601   static AliAODEvent* aod = 0;
 602
 603   // take all other information from the aod we take the tracks from
 604   if(!aod){
 605    if(fUseAODTrackInput)aod = fAODIn;
 606    else aod = fAODOut;
 607   }
 608
 609
 610   //
 611   // Execute analysis for current event
 612   //
 613   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 614   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 615   if(!aodH){
 616     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 617     return;
 618   }
 619
 620   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 621   TClonesArray *aodRecJets = 0;
 622
 623   if(fAODOut&&!aodRecJets){
 624     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
 625   }
 626   if(fAODExtension&&!aodRecJets){
 627     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
 628   }
 629   if(fAODIn&&!aodRecJets){
 630     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
 631   }
 632
 633
 634
 635   if(!aodRecJets){
 636     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 637     return;
 638   }
 639
 640   TClonesArray *aodGenJets = 0;
 641   if(fBranchGen.Length()>0){
 642     if(fAODOut&&!aodGenJets){
 643       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
 644     }
 645     if(fAODExtension&&!aodGenJets){
 646       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
 647     }
 648     if(fAODIn&&!aodGenJets){
 649       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
 650     }
 651
 652     if(!aodGenJets){
 653       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
 654       return;
 655     }
 656   }
 657
 658   TClonesArray *aodBackRecJets = 0;
 659   if(fBranchBkgRec.Length()>0){
 660     if(fAODOut&&!aodBackRecJets){
 661       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
 662     }
 663     if(fAODExtension&&!aodBackRecJets){
 664       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
 665     }
 666     if(fAODIn&&!aodBackRecJets){
 667       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
 668     }
 669
 670     if(!aodBackRecJets){
 671       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
 672       return;
 673     }
 674   }
 675
 676
 677   TClonesArray *aodBackGenJets = 0;
 678
 679   if(fBranchBkgGen.Length()>0){
 680     if(fAODOut&&!aodBackGenJets){
 681       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
 682     }
 683     if(fAODExtension&&!aodBackGenJets){
 684       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
 685     }
 686     if(fAODIn&&!aodBackGenJets){
 687       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
 688     }
 689
 690     if(!aodBackGenJets){
 691       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
 692       return;
 693     }
 694   }
 695
 696
 697   // new Scheme
 698   // first fill all the pure  histograms, i.e. full jets
 699   // in case of the correaltion limit to the n-leading jets
 700
 701   // reconstructed
 702
 703
 704   // generated
 705
 706
 707   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
 708
 709
 710
 711   TList genJetsList;         // full acceptance
 712   TList genJetsListCut;      // acceptance cut
 713   TList recJetsList;         // full acceptance
 714   TList recJetsListCut;      // acceptance cut
 715
 716
 717   GetListOfJets(&recJetsList,aodRecJets,0);
 718   GetListOfJets(&recJetsListCut,aodRecJets,1);
 719
 720   if(fBranchGen.Length()>0){
 721     GetListOfJets(&genJetsList,aodGenJets,0);
 722     GetListOfJets(&genJetsListCut,aodGenJets,1);
 723   }
 724
 725   Double_t eventW = 1;
 726   Double_t ptHard = 0;
 727   Double_t nTrials = 1; // Trials for MC trigger
 728   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 729
 730   // Getting some global properties
 731   fCentrality = GetCentrality();
 732   fh1Centrality->Fill(fCentrality);
 733
 734
 735   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 736     // this is the part we only use when we have MC information
 737     AliMCEvent* mcEvent = MCEvent();
 738     //    AliStack *pStack = 0;
 739     if(!mcEvent){
 740       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 741       return;
 742     }
 743     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 744     if(pythiaGenHeader){
 745       nTrials = pythiaGenHeader->Trials();
 746       ptHard  = pythiaGenHeader->GetPtHard();
 747       int iProcessType = pythiaGenHeader->ProcessType();
 748       // 11 f+f -> f+f
 749       // 12 f+barf -> f+barf
 750       // 13 f+barf -> g+g
 751       // 28 f+g -> f+g
 752       // 53 g+g -> f+barf
 753       // 68 g+g -> g+g
 754       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 755       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 756     }
 757   }// (fAnalysisType&kMCESD)==kMCESD)
 758
 759
 760   // we simply fetch the tracks/mc particles as a list of AliVParticles
 761
 762   TList recParticles;
 763   TList genParticles;
 764
 765   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 766   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 767   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 768   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 769
 770   CalculateReactionPlaneAngle(&recParticles);
 771
 772   // Event control and counting ...
 773   // MC
 774   fh1PtHard->Fill(ptHard,eventW);
 775   fh1PtHardNoW->Fill(ptHard,1);
 776   fh1PtHardTrials->Fill(ptHard,nTrials);
 777
 778   // Real
 779   if(aod->GetPrimaryVertex()){// No vtx for pure MC
 780     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
 781   }
 782
 783
 784   Int_t recMult1 = recParticles.GetEntries();
 785   Int_t genMult1 = genParticles.GetEntries();
 786
 787   Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
 788   Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
 789
 790   fh2MultRec->Fill(recMult1,recMult2);
 791   fh2MultGen->Fill(genMult1,genMult2);
 792   fMultRec = recMult1;
 793   if(fMultRec<=0)fMultRec = recMult2;
 794   fMultGen = genMult1;
 795   if(fMultGen<=0)fMultGen = genMult2;
 796
 797   // the loops for rec and gen should be indentical... pass it to a separate
 798   // function ...
 799   // Jet Loop
 800   // Track Loop
 801   // Jet Jet Loop
 802   // Jet Track loop
 803
 804   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 805   FillJetHistos(recJetsList,recParticles,kJetRecFull);
 806   FillTrackHistos(recParticles,kJetRec);
 807
 808   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 809   FillJetHistos(genJetsList,genParticles,kJetGenFull);
 810   FillTrackHistos(genParticles,kJetGen);
 811
 812   // Here follows the jet matching part
 813   // delegate to a separated method?
 814
 815   FillMatchHistos(recJetsList,genJetsList);
 816
 817   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 818   PostData(1, fHistList);
 819 }
 820
 821 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 822
 823   if(iType>=kJetTypes){
 824     return;
 825   }
 826   if(!fFlagJetType[iType])return;
 827
 828   Int_t refMult = fMultRec;
 829   if(iType==kJetGen||iType==kJetGenFull){
 830     refMult = fMultGen;
 831   }
 832
 833   Int_t nJets = jetsList.GetEntries();
 834   fh1NJets[iType]->Fill(nJets);
 835
 836   if(nJets<=0)return;
 837
 838   Float_t ptOld = 1.E+32;
 839   Float_t pT0 = 0;
 840   Float_t pT1 = 0;
 841   Float_t phi0 = 0;
 842   Float_t phi1 = 0;
 843   Int_t ij0 = -1;
 844   Int_t ij1 = -1;
 845
 846
 847   for(int ij = 0;ij < nJets;ij++){
 848     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
 849     Float_t ptJet = jet->Pt();
 850     fh1PtJetsIn[iType]->Fill(ptJet);
 851     if(ptJet>ptOld){
 852       Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
 853     }
 854     ptOld = ptJet;
 855
 856     // find the dijets assume sorting and acceptance cut...
 857     if(ij==0){
 858       ij0 = ij;
 859       pT0 = ptJet;
 860       phi0 = jet->Phi();
 861       if(phi0<0)phi0+=TMath::Pi()*2.;
 862     }
 863     else if(ptJet>pT1){
 864       // take only the backward hemisphere??
 865       phi1 = jet->Phi();
 866       if(phi1<0)phi1+=TMath::Pi()*2.;
 867       Float_t dPhi = phi1 - phi0;
 868       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 869       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 870       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
 871         ij1 = ij;
 872         pT1 = ptJet;
 873       }
 874     }
 875     // fill jet histos for kmax jets
 876
 877       Float_t phiJet = jet->Phi();
 878       Float_t etaJet = jet->Eta();
 879       if(phiJet<0)phiJet+=TMath::Pi()*2.;
 880       fh1TmpRho->Reset();
 881       if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
 882       fh1PtIn[iType][kMaxJets]->Fill(ptJet);
 883       // fill leading jets...
 884       AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
 885       //      AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
 886       Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
 887       if(ptJet>10){
 888         if(ij<kMaxJets){
 889           fh2PhiEta[iType][ij]->Fill(phiJet,etaJet);
 890           fh2AreaPt[iType][ij]->Fill(jet->EffectiveAreaCharged(),ptJet);
 891           fh2EtaArea[iType][ij]->Fill(etaJet,jet->EffectiveAreaCharged());
 892         }
 893         fh2PhiEta[iType][kMaxJets]->Fill(phiJet,etaJet);
 894         fh2AreaPt[iType][kMaxJets]->Fill(jet->EffectiveAreaCharged(),ptJet);
 895         fh2EtaArea[iType][kMaxJets]->Fill(etaJet,jet->EffectiveAreaCharged());
 896       }
 897       if(ij<kMaxJets){
 898         fh3MultPtRP[iType][ij]->Fill(refMult,ptJet,phiBin);
 899         fh2PhiPt[iType][ij]->Fill(phiJet,ptJet);
 900         fh2EtaPt[iType][ij]->Fill(etaJet,ptJet);
 901         if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
 902       }
 903       fh3MultPtRP[iType][kMaxJets]->Fill(refMult,ptJet,phiBin);
 904       fh2PhiPt[iType][kMaxJets]->Fill(phiJet,ptJet);
 905       fh2EtaPt[iType][kMaxJets]->Fill(etaJet,ptJet);
 906       if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
 907
 908       if(particlesList.GetSize()&&ij<kMaxJets){
 909         // Particles... correlated with jets...
 910         for(int it = 0;it<particlesList.GetEntries();++it){
 911           AliVParticle *part = (AliVParticle*)particlesList.At(it);
 912           Float_t deltaR = jet->DeltaR(part);
 913           if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
 914         }
 915         // fill the jet shapes
 916         Float_t rhoSum = 0;
 917         for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
 918           Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
 919           Float_t rho = fh1TmpRho->GetBinContent(ibx);
 920           rhoSum += rho;
 921           fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
 922           fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
 923         }
 924       }// if we have particles
 925   }// Jet Loop
 926
 927
 928   // do something with dijets...
 929   if(ij0>=0&&ij1>0){
 930     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
 931     Double_t ptJet0 = jet0->Pt();
 932     Double_t phiJet0 = jet0->Phi();
 933     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
 934
 935     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
 936     Double_t ptJet1 = jet1->Pt();
 937     Double_t phiJet1 = jet1->Phi();
 938     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
 939
 940     Float_t deltaPhi = phiJet0 - phiJet1;
 941     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
 942     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
 943     deltaPhi = TMath::Abs(deltaPhi);
 944     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
 945
 946     Float_t asym = 9999;
 947     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
 948       fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
 949       fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
 950       fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
 951       Float_t minv = 2.*(jet0->P()*jet1->P()-
 952                          jet0->Px()*jet1->Px()-
 953                          jet0->Py()*jet1->Py()-
 954                          jet0->Pz()*jet1->Pz());    // assume mass == 0;
 955       if(minv<0)minv=0; // prevent numerical instabilities
 956       minv = TMath::Sqrt(minv);
 957       fh1DijetMinv[iType]->Fill(minv);
 958   }
 959
 960
 961
 962   // count down the jets above thrueshold
 963   Int_t nOver = nJets;
 964   if(nOver>0){
 965     TIterator *jetIter = jetsList.MakeIterator();
 966     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
 967     if(tmpJet){
 968       Float_t pt = tmpJet->Pt();
 969       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
 970         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
 971         while(pt<ptCut&&tmpJet){
 972           nOver--;
 973           tmpJet = (AliAODJet*)(jetIter->Next());
 974           if(tmpJet){
 975             pt = tmpJet->Pt();
 976           }
 977         }
 978         if(nOver<=0)break;
 979         fh2NJetsPt[iType]->Fill(ptCut,nOver);
 980       }
 981     }
 982     delete jetIter;
 983   }
 984 }
 985
 986 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
 987
 988   if(!fFlagJetType[iType])return;
 989   Int_t refMult = fMultRec;
 990   if(iType==kJetGen||iType==kJetGenFull){
 991     refMult = fMultGen;
 992
 993   }
 994
 995   Int_t nTrackOver = particlesList.GetSize();
 996   // do the same for tracks and jets
 997   if(nTrackOver>0){
 998     TIterator *trackIter = particlesList.MakeIterator();
 999     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1000     Float_t pt = tmpTrack->Pt();
1001     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1002       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1003       while(pt<ptCut&&tmpTrack){
1004         nTrackOver--;
1005         tmpTrack = (AliVParticle*)(trackIter->Next());
1006         if(tmpTrack){
1007           pt = tmpTrack->Pt();
1008         }
1009       }
1010       if(nTrackOver<=0)break;
1011       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1012     }
1013
1014
1015     trackIter->Reset();
1016     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1017     Float_t sumPt = 0;
1018
1019     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1020       Float_t tmpPt = tmpTrack->Pt();
1021       fh1PtTracksIn[iType]->Fill(tmpPt);
1022       fh1PtTracksInLow[iType]->Fill(tmpPt);
1023
1024       sumPt += tmpPt;
1025       Float_t tmpPhi = tmpTrack->Phi();
1026       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1027       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1028       fh3MultTrackPtRP[iType]->Fill(refMult,tmpPt,phiBin);
1029
1030       if(tmpTrack==leading){
1031         fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
1032         fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
1033         continue;
1034       }
1035     }
1036     fh1SumPtTrack[iType]->Fill(sumPt);
1037     delete trackIter;
1038   }
1039
1040 }
1041
1042
1043 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1044
1045
1046   // Fill al the matching histos
1047   // It is important that the acceptances for the mathing are as large as possible
1048   // to avoid false matches on the edge of acceptance
1049   // therefore we add some extra matching jets as overhead
1050
1051   static TArrayI aGenIndex(recJetsList.GetEntries());
1052   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1053
1054   static TArrayI aRecIndex(genJetsList.GetEntries());
1055   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1056
1057   if(fDebug){
1058     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1059     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1060   }
1061   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1062                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1063                                             aGenIndex,aRecIndex,fDebug);
1064
1065   if(fDebug){
1066     for(int i = 0;i< aGenIndex.GetSize();++i){
1067       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1068     }
1069     for(int i = 0;i< aRecIndex.GetSize();++i){
1070       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1071     }
1072   }
1073
1074   Double_t container[6];
1075   Double_t containerPhiZ[6];
1076
1077   // loop over generated jets
1078   // consider the
1079   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1080     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1081     Double_t ptGen = genJet->Pt();
1082     Double_t phiGen = genJet->Phi();
1083     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1084     Double_t etaGen = genJet->Eta();
1085     container[3] = ptGen;
1086     container[4] = etaGen;
1087     container[5] = phiGen;
1088     fhnJetContainer[kStep0]->Fill(&container[3]);
1089     if(JetSelected(genJet)){
1090       fhnJetContainer[kStep1]->Fill(&container[3]);
1091       Int_t ir = aRecIndex[ig];
1092       if(ir>=0&&ir<recJetsList.GetEntries()){
1093         fhnJetContainer[kStep2]->Fill(&container[3]);
1094         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1095         if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
1096         if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
1097       }
1098     }
1099   }// loop over generated jets used for matching...
1100
1101
1102
1103   // loop over reconstructed jets
1104   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1105     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1106     Double_t etaRec = recJet->Eta();
1107     Double_t ptRec = recJet->Pt();
1108     Double_t phiRec = recJet->Phi();
1109     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1110     // do something with dijets...
1111
1112     container[0] = ptRec;
1113     container[1] = etaRec;
1114     container[2] = phiRec;
1115     containerPhiZ[0] = ptRec;
1116     containerPhiZ[1] = phiRec;
1117
1118     fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1119     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1120
1121     if(JetSelected(recJet)){
1122       fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1123       // Fill Correlation
1124       Int_t ig = aGenIndex[ir];
1125       if(ig>=0 && ig<genJetsList.GetEntries()){
1126         fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1127         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1128         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1129         Double_t ptGen  = genJet->Pt();
1130         Double_t phiGen = genJet->Phi();
1131         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1132         Double_t etaGen = genJet->Eta();
1133
1134         container[3] = ptGen;
1135         container[4] = etaGen;
1136         container[5] = phiGen;
1137         containerPhiZ[3] = ptGen;
1138         //
1139         // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1140         //
1141         if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1142         fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1143         fhnCorrelation->Fill(container);
1144         if(ptGen>0){
1145           Float_t delta = (ptRec-ptGen)/ptGen;
1146           fh2RelPtFGen->Fill(ptGen,delta);
1147           fh2PtFGen->Fill(ptGen,ptRec);
1148         }
1149         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1150       }
1151       else{
1152         containerPhiZ[3] = 0;
1153         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1154       }
1155     }// loop over reconstructed jets
1156   }
1157   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1158 }
1159
1160
1161 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1162   //
1163   // Create the particle container for the correction framework manager and
1164   // link it
1165   //
1166   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1167   const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1168   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1169   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1170   const Double_t kZmin = 0., kZmax = 1;
1171
1172   // can we neglect migration in eta and phi?
1173   // phi should be no problem since we cover full phi and are phi symmetric
1174   // eta migration is more difficult  due to needed acceptance correction
1175   // in limited eta range
1176
1177   //arrays for the number of bins in each dimension
1178   Int_t iBin[kNvar];
1179   iBin[0] = 160; //bins in pt
1180   iBin[1] =  1; //bins in eta
1181   iBin[2] = 1; // bins in phi
1182
1183   //arrays for lower bounds :
1184   Double_t* binEdges[kNvar];
1185   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1186     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1187
1188   //values for bin lower bounds
1189   //  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
1190   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1191   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1192   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1193
1194
1195   for(int i = 0;i < kMaxStep*2;++i){
1196     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1197     for (int k=0; k<kNvar; k++) {
1198       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1199     }
1200   }
1201   //create correlation matrix for unfolding
1202   Int_t thnDim[2*kNvar];
1203   for (int k=0; k<kNvar; k++) {
1204     //first half  : reconstructed
1205     //second half : MC
1206     thnDim[k]      = iBin[k];
1207     thnDim[k+kNvar] = iBin[k];
1208   }
1209
1210   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1211   for (int k=0; k<kNvar; k++) {
1212     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1213     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1214   }
1215   fhnCorrelation->Sumw2();
1216
1217
1218   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1219     delete [] binEdges[ivar];
1220
1221
1222
1223   // for second correlation histogram
1224
1225
1226   const Int_t kNvarPhiZ   = 4;
1227   //arrays for the number of bins in each dimension
1228   Int_t iBinPhiZ[kNvarPhiZ];
1229   iBinPhiZ[0] = 80; //bins in pt
1230   iBinPhiZ[1] = 72; //bins in phi
1231   iBinPhiZ[2] = 20; // bins in Z
1232   iBinPhiZ[3] = 80; //bins in ptgen
1233
1234
1235   return;
1236   //arrays for lower bounds :
1237   Double_t* binEdgesPhiZ[kNvarPhiZ];
1238   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1239     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1240
1241   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1242   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1243   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1244   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1245
1246   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1247   for (int k=0; k<kNvarPhiZ; k++) {
1248     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1249   }
1250   fhnCorrelationPhiZRec->Sumw2();
1251
1252   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1253     delete [] binEdgesPhiZ[ivar];
1254
1255 }
1256
1257 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1258 {
1259 // Terminate analysis
1260 //
1261     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1262 }
1263
1264
1265 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1266
1267   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1268
1269   Int_t iCount = 0;
1270   if(type==kTrackAOD){
1271     AliAODEvent *aod = 0;
1272     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1273     else aod = AODEvent();
1274     if(!aod){
1275       return iCount;
1276     }
1277     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1278       AliAODTrack *tr = aod->GetTrack(it);
1279       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1280       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1281       if(tr->Pt()<fMinTrackPt)continue;
1282       if(fDebug>0){
1283         if(tr->Pt()>20){
1284           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1285           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1286           tr->Print();
1287           //    tr->Dump();
1288           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1289           if(fESD){
1290             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1291             esdTr->Print("");
1292             // esdTr->Dump();
1293           }
1294         }
1295       }
1296       list->Add(tr);
1297       iCount++;
1298     }
1299   }
1300   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1301     AliMCEvent* mcEvent = MCEvent();
1302     if(!mcEvent)return iCount;
1303     // we want to have alivpartilces so use get track
1304     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1305       if(!mcEvent->IsPhysicalPrimary(it))continue;
1306       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1307       if(part->Pt()<fMinTrackPt)continue;
1308       if(type == kTrackKineAll){
1309         list->Add(part);
1310         iCount++;
1311       }
1312       else if(type == kTrackKineCharged){
1313         if(part->Particle()->GetPDG()->Charge()==0)continue;
1314         list->Add(part);
1315         iCount++;
1316       }
1317     }
1318   }
1319   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1320     AliAODEvent *aod = 0;
1321     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1322     else aod = AODEvent();
1323     if(!aod)return iCount;
1324     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1325     if(!tca)return iCount;
1326     for(int it = 0;it < tca->GetEntriesFast();++it){
1327       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1328       if(!part)continue;
1329       if(part->Pt()<fMinTrackPt)continue;
1330       if(!part->IsPhysicalPrimary())continue;
1331       if(type == kTrackAODMCAll){
1332         list->Add(part);
1333         iCount++;
1334       }
1335       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1336         if(part->Charge()==0)continue;
1337         if(kTrackAODMCCharged){
1338           list->Add(part);
1339         }
1340         else {
1341           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1342           list->Add(part);
1343         }
1344         iCount++;
1345       }
1346     }
1347   }// AODMCparticle
1348   list->Sort();
1349   return iCount;
1350
1351 }
1352
1353
1354 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1355     AliAODEvent *aod = 0;
1356     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1357     else aod = AODEvent();
1358     if(!aod){
1359       return 100;
1360     }
1361     return aod->GetHeader()->GetCentrality();
1362 }
1363
1364
1365
1366 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1367   Bool_t selected = false;
1368
1369   if(!jet)return selected;
1370
1371   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1372     selected = kTRUE;
1373   }
1374   return selected;
1375
1376 }
1377
1378 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1379
1380   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1381   Int_t iCount = 0;
1382
1383   if(!jarray){
1384     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1385     return iCount;
1386   }
1387
1388
1389   for(int ij=0;ij<jarray->GetEntries();ij++){
1390     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1391     if(!jet)continue;
1392     if(type==0){
1393       // No cut at all, main purpose here is sorting
1394       list->Add(jet);
1395       iCount++;
1396     }
1397     else if(type == 1){
1398       // eta cut
1399       if(JetSelected(jet)){
1400         list->Add(jet);
1401         iCount++;
1402       }
1403     }
1404   }
1405
1406   list->Sort();
1407   return iCount;
1408
1409 }
1410
1411
1412 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1413   if(!jets)return 0;
1414
1415   Int_t refMult = 0;
1416   for(int ij = 0;ij < jets->GetEntries();++ij){
1417     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1418     if(!jet)continue;
1419     TRefArray *refs = jet->GetRefTracks();
1420     if(!refs)continue;
1421     refMult += refs->GetEntries();
1422   }
1423   return refMult;
1424
1425 }
1426
1427
1428 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1429   if(!jet)return 0;
1430   TRefArray *refs = jet->GetRefTracks();
1431   if(!refs) return 0;
1432   AliVParticle *leading = 0;
1433   Float_t fMaxPt = 0;
1434   for(int i = 0;i<refs->GetEntries();i++){
1435     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1436     if(!tmp)continue;
1437     if(tmp->Pt()>fMaxPt){
1438       leading = tmp;
1439       fMaxPt = tmp->Pt();
1440     }
1441   }
1442   return leading;
1443 }
1444
1445
1446 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1447   if(!jet)return 0;
1448   if(!list) return 0;
1449   AliVParticle *leading = 0;
1450   Float_t fMaxPt = 0;
1451   for(int i = 0;i<list->GetEntries();i++){
1452     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1453     if(!tmp)continue;
1454     if(jet->DeltaR(tmp)>r)continue;
1455     if(tmp->Pt()>fMaxPt){
1456       leading = tmp;
1457       fMaxPt = tmp->Pt();
1458     }
1459   }
1460   return leading;
1461 }
1462
1463 Bool_t AliAnalysisTaskJetSpectrum2::CalculateReactionPlaneAngle(const TList *trackList)
1464 {
1465
1466   if(!trackList)return kFALSE;
1467   fRPAngle=0;
1468
1469   // need to get this info from elsewhere??
1470   Double_t fFlatA[2] = {1,1};
1471   Double_t fFlatB[2] = {1,1};
1472
1473
1474   Double_t fPsiRP =0,fDeltaPsiRP = 0;
1475
1476
1477
1478   TVector2 mQ,mQ1,mQ2;
1479   Float_t mQx=0, mQy=0;
1480
1481   Float_t mQx1=0, mQy1=0;
1482   Float_t mQx2=0, mQy2=0;
1483
1484   AliVParticle *track=0x0;
1485   Int_t count[3]={0,0,0};
1486
1487
1488   for (Int_t iter=0;iter<trackList->GetEntries();iter++){
1489
1490     track=(AliVParticle*)trackList->At(iter);
1491
1492     //cuts already applied before
1493     // Comment DCA not correctly implemented yet for AOD tracks
1494
1495     Double_t momentum;
1496     if(track->Charge()>0){momentum=track->Pt();}
1497     else{momentum=-track->Pt();}
1498
1499
1500
1501     // For Weighting
1502     fh3RPPhiTracks->Fill(fCentrality,momentum,track->Phi());
1503     count[0]++;
1504
1505     //    Double_t phiweight=GetPhiWeight(track->Phi(),momentum);
1506     Double_t phiweight=1;
1507     Double_t weight=2;
1508     if(track->Pt()<2){weight=track->Pt();}
1509
1510
1511     mQx += (cos(2*track->Phi()))*weight*phiweight;
1512     mQy += (sin(2*track->Phi()))*weight*phiweight;
1513
1514     // Make random Subevents
1515
1516     if(fRPSubeventMethod==0){
1517       if(fRandomizer->Binomial(1,0.5)){
1518         mQx1 += (cos(2*track->Phi()))*weight*phiweight;
1519         mQy1 += (sin(2*track->Phi()))*weight*phiweight;
1520         count[1]++;}
1521       else{
1522         mQx2 += (cos(2*track->Phi()))*weight*phiweight;
1523         mQy2 += (sin(2*track->Phi()))*weight*phiweight;
1524         count[2]++;}
1525     }
1526     else if(fRPSubeventMethod==1){
1527       // Make eta dependent subevents
1528       if(track->Eta()>0){
1529         mQx1 += (cos(2*track->Phi()))*weight*phiweight;
1530         mQy1 += (sin(2*track->Phi()))*weight*phiweight;
1531         count[1]++;}
1532       else{
1533         mQx2 += (cos(2*track->Phi()))*weight*phiweight;
1534         mQy2 += (sin(2*track->Phi()))*weight*phiweight;
1535         count[2]++;}
1536     }
1537
1538   }
1539
1540
1541
1542   //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
1543   if(count[0]==0||count[1]==0||count[2]==0){
1544     return kFALSE;
1545   }
1546
1547   mQ.Set(mQx,mQy);
1548   mQ1.Set(mQx1,mQy1);
1549   mQ2.Set(mQx2,mQy2);
1550
1551   // cout<<"MQ"<<mQx<<" " <<mQy<<" psi"<<endl;
1552
1553   fPsiRP=mQ.Phi()/2;
1554
1555   //Correction
1556   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);
1557
1558   Double_t fPsiRP1=mQ1.Phi()/2;
1559   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);
1560   Double_t fPsiRP2=mQ2.Phi()/2;
1561   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);
1562   fDeltaPsiRP=fPsiRP1-fPsiRP2;
1563
1564   if(fPsiRP>TMath::Pi()){fPsiRP-=TMath::Pi();}
1565   if(fPsiRP<0){fPsiRP+=TMath::Pi();}
1566
1567   // reactionplaneangle + Pi() is the same angle
1568   if(TMath::Abs(fDeltaPsiRP)>TMath::Pi()/2){
1569     if(fDeltaPsiRP>0)fDeltaPsiRP-=TMath::Pi();
1570     else fDeltaPsiRP+=TMath::Pi();
1571   }
1572
1573   Double_t cos2deltaRP=TMath::Cos(2*fDeltaPsiRP);
1574
1575   // FillHistograms
1576   fh2RPSubevents->Fill(fPsiRP1,fPsiRP2);
1577   fh1RP->Fill(fPsiRP);
1578   fh2RPCentrality->Fill(fCentrality,fPsiRP);
1579   fh2RPDeltaRP->Fill(fDeltaPsiRP,fCentrality);
1580   fh2RPQxQy->Fill(mQx,mQy);
1581   fh2RPCosDeltaRP->Fill(fCentrality,cos2deltaRP);
1582
1583   fRPAngle=fPsiRP;
1584
1585   return kTRUE;
1586 }
1587
1588 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1589 {
1590     Int_t phibin=-1;
1591     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1592     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1593     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1594     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1595     return phibin;
1596 }
1597
1598  //________________________________________________________________________