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