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