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