PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx

   1 // **************************************
   2 // Task 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 <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TFile.h>
  29 #include <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TList.h>
  35 #include <TLorentzVector.h>
  36 #include <TClonesArray.h>
  37 #include  "TDatabasePDG.h"
  38
  39 #include "AliAnalysisTaskJetSpectrum2.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliJetFinder.h"
  42 #include "AliJetHeader.h"
  43 #include "AliJetReader.h"
  44 #include "AliJetReaderHeader.h"
  45 #include "AliUA1JetHeaderV1.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODTrack.h"
  50 #include "AliAODJet.h"
  51 #include "AliAODJetEventBackground.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliJetKineReaderHeader.h"
  58 #include "AliGenCocktailEventHeader.h"
  59 #include "AliInputEventHandler.h"
  60
  61
  62 #include "AliAnalysisHelperJetTasks.h"
  63
  64 ClassImp(AliAnalysisTaskJetSpectrum2)
  65
  66 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
  67                                                             fJetHeaderRec(0x0),
  68                                                             fJetHeaderGen(0x0),
  69                                                             fAOD(0x0),
  70                                                             fhnCorrelation(0x0),
  71   fhnCorrelationPhiZRec(0x0),
  72   f1PtScale(0x0),
  73                                                             fBranchRec("jets"),
  74                                                             fBranchGen(""),
  75                                                             fBranchBkg(""),
  76                                                             fUseAODJetInput(kFALSE),
  77                                                             fUseAODTrackInput(kFALSE),
  78                                                             fUseAODMCInput(kFALSE),
  79                                                             fUseGlobalSelection(kFALSE),
  80                                                             fUseExternalWeightOnly(kFALSE),
  81                                                             fLimitGenJetEta(kFALSE),
  82                                                             fBkgSubtraction(kFALSE),
  83                                                             fFillCorrBkg(0),
  84                                                             fFilterMask(0),
  85   fEventSelectionMask(0),
  86                                                             fAnalysisType(0),
  87   fTrackTypeRec(kTrackUndef),
  88   fTrackTypeGen(kTrackUndef),
  89   fEventClass(0),
  90   fAvgTrials(1),
  91   fExternalWeight(1),
  92                                                             fRecEtaWindow(0.5),
  93   fMinJetPt(0),
  94                                                             fDeltaPhiWindow(20./180.*TMath::Pi()),
  95   fh1Xsec(0x0),
  96   fh1Trials(0x0),
  97   fh1PtHard(0x0),
  98   fh1PtHardNoW(0x0),
  99   fh1PtHardTrials(0x0),
 100   fh1ZVtx(0x0),
 101   fh1NGenJets(0x0),
 102   fh1NRecJets(0x0),
 103   fh1PtTrackRec(0x0),
 104   fh1SumPtTrackRec(0x0),
 105   fh1SumPtTrackAreaRec(0x0),
 106   fh1TmpRho(0x0),
 107   fh1PtJetsRecIn(0x0),
 108   fh1PtJetsLeadingRecIn(0x0),
 109   fh1PtTracksRecIn(0x0),
 110   fh1PtTracksLeadingRecIn(0x0),
 111   fh1PtTracksGenIn(0x0),
 112   fh2NRecJetsPt(0x0),
 113   fh2NRecTracksPt(0x0),
 114   fh2JetsLeadingPhiEta(0x0),
 115   fh2JetsLeadingPhiPt(0x0),
 116   fh2TracksLeadingPhiEta(0x0),
 117   fh2TracksLeadingPhiPt(0x0),
 118   fh2TracksLeadingJetPhiPt(0x0),
 119   fh2JetPtJetPhi(0x0),
 120   fh2TrackPtTrackPhi(0x0),
 121   fh2RelPtFGen(0x0),
 122   fh2DijetDeltaPhiPt(0x0),
 123   fh2DijetAsymPt(0x0),
 124   fh2DijetAsymPtCut(0x0),
 125   fh2DijetDeltaPhiDeltaEta(0x0),
 126   fh2DijetPt2vsPt1(0x0),
 127   fh2DijetDifvsSum(0x0),
 128   fh1DijetMinv(0x0),
 129   fh1DijetMinvCut(0x0),
 130   fh1Bkg1(0x0),
 131   fh1Bkg2(0x0),
 132   fh1Bkg3(0x0),
 133   fh1Sigma1(0x0),
 134   fh1Sigma2(0x0),
 135   fh1Sigma3(0x0),
 136   fh1Area1(0x0),
 137   fh1Area2(0x0),
 138   fh1Area3(0x0),
 139   fh1Ptjet(0x0),
 140   fh1Ptjetsub1(0x0),
 141   fh1Ptjetsub2(0x0),
 142   fh1Ptjetsub3(0x0),
 143   fh1Ptjethardest(0x0),
 144   fh1Ptjetsubhardest1(0x0),
 145   fh1Ptjetsubhardest2(0x0),
 146   fh1Ptjetsubhardest3(0x0),
 147   fh2Rhovspthardest1(0x0),
 148   fh2Rhovspthardest2(0x0),
 149   fh2Rhovspthardest3(0x0),
 150   fh2Errorvspthardest1(0x0),
 151   fh2Errorvspthardest2(0x0),
 152   fh2Errorvspthardest3(0x0),
 153   fHistList(0x0)
 154 {
 155   for(int i = 0;i < kMaxStep*2;++i){
 156     fhnJetContainer[i] = 0;
 157   }
 158   for(int i = 0;i < kMaxJets;++i){
 159     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 160
 161     fh2PhiPt[i] = 0;
 162     fh2PhiEta[i] = 0;
 163     fh2RhoPtRec[i] = 0;
 164     fh2RhoPtGen[i] = 0;
 165     fh2PsiPtGen[i] = 0;
 166     fh2PsiPtRec[i] = 0;
 167     fh2FragRec[i] = 0;
 168     fh2FragLnRec[i] = 0;
 169     fh2FragGen[i] = 0;
 170     fh2FragLnGen[i] = 0;
 171   }
 172
 173 }
 174
 175 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 176   AliAnalysisTaskSE(name),
 177   fJetHeaderRec(0x0),
 178   fJetHeaderGen(0x0),
 179   fAOD(0x0),
 180   fhnCorrelation(0x0),
 181   fhnCorrelationPhiZRec(0x0),
 182   f1PtScale(0x0),
 183   fBranchRec("jets"),
 184   fBranchGen(""),
 185   fBranchBkg(""),
 186   fUseAODJetInput(kFALSE),
 187   fUseAODTrackInput(kFALSE),
 188   fUseAODMCInput(kFALSE),
 189   fUseGlobalSelection(kFALSE),
 190   fUseExternalWeightOnly(kFALSE),
 191   fLimitGenJetEta(kFALSE),
 192   fBkgSubtraction(kFALSE),
 193   fFillCorrBkg(0),
 194   fFilterMask(0),
 195   fEventSelectionMask(0),
 196   fAnalysisType(0),
 197   fTrackTypeRec(kTrackUndef),
 198   fTrackTypeGen(kTrackUndef),
 199   fEventClass(0),
 200   fAvgTrials(1),
 201   fExternalWeight(1),
 202   fRecEtaWindow(0.5),
 203   fMinJetPt(0),
 204   fDeltaPhiWindow(20./180.*TMath::Pi()),
 205   fh1Xsec(0x0),
 206   fh1Trials(0x0),
 207   fh1PtHard(0x0),
 208   fh1PtHardNoW(0x0),
 209   fh1PtHardTrials(0x0),
 210   fh1ZVtx(0x0),
 211   fh1NGenJets(0x0),
 212   fh1NRecJets(0x0),
 213   fh1PtTrackRec(0x0),
 214   fh1SumPtTrackRec(0x0),
 215   fh1SumPtTrackAreaRec(0x0),
 216   fh1TmpRho(0x0),
 217   fh1PtJetsRecIn(0x0),
 218   fh1PtJetsLeadingRecIn(0x0),
 219   fh1PtTracksRecIn(0x0),
 220   fh1PtTracksLeadingRecIn(0x0),
 221   fh1PtTracksGenIn(0x0),
 222   fh2NRecJetsPt(0x0),
 223   fh2NRecTracksPt(0x0),
 224   fh2JetsLeadingPhiEta(0x0),
 225   fh2JetsLeadingPhiPt(0x0),
 226   fh2TracksLeadingPhiEta(0x0),
 227   fh2TracksLeadingPhiPt(0x0),
 228   fh2TracksLeadingJetPhiPt(0x0),
 229   fh2JetPtJetPhi(0x0),
 230   fh2TrackPtTrackPhi(0x0),
 231   fh2RelPtFGen(0x0),
 232   fh2DijetDeltaPhiPt(0x0),
 233   fh2DijetAsymPt(0x0),
 234   fh2DijetAsymPtCut(0x0),
 235   fh2DijetDeltaPhiDeltaEta(0x0),
 236   fh2DijetPt2vsPt1(0x0),
 237   fh2DijetDifvsSum(0x0),
 238   fh1DijetMinv(0x0),
 239   fh1DijetMinvCut(0x0),
 240   fh1Bkg1(0x0),
 241   fh1Bkg2(0x0),
 242   fh1Bkg3(0x0),
 243   fh1Sigma1(0x0),
 244   fh1Sigma2(0x0),
 245   fh1Sigma3(0x0),
 246   fh1Area1(0x0),
 247   fh1Area2(0x0),
 248   fh1Area3(0x0),
 249   fh1Ptjet(0x0),
 250   fh1Ptjetsub1(0x0),
 251   fh1Ptjetsub2(0x0),
 252   fh1Ptjetsub3(0x0),
 253   fh1Ptjethardest(0x0),
 254   fh1Ptjetsubhardest1(0x0),
 255   fh1Ptjetsubhardest2(0x0),
 256   fh1Ptjetsubhardest3(0x0),
 257   fh2Rhovspthardest1(0x0),
 258   fh2Rhovspthardest2(0x0),
 259   fh2Rhovspthardest3(0x0),
 260   fh2Errorvspthardest1(0x0),
 261   fh2Errorvspthardest2(0x0),
 262   fh2Errorvspthardest3(0x0),
 263   fHistList(0x0)
 264 {
 265
 266   for(int i = 0;i < kMaxStep*2;++i){
 267     fhnJetContainer[i] = 0;
 268   }
 269   for(int i = 0;i < kMaxJets;++i){
 270     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 271
 272     fh2PhiPt[i] = 0;
 273     fh2PhiEta[i] = 0;
 274     fh2RhoPtRec[i] = 0;
 275     fh2RhoPtGen[i] = 0;
 276     fh2PsiPtGen[i] = 0;
 277     fh2PsiPtRec[i] = 0;
 278     fh2FragRec[i] = 0;
 279     fh2FragLnRec[i] = 0;
 280     fh2FragGen[i] = 0;
 281     fh2FragLnGen[i] = 0;
 282   }
 283   DefineOutput(1, TList::Class());
 284 }
 285
 286
 287
 288 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 289 {
 290   //
 291   // Implemented Notify() to read the cross sections
 292   // and number of trials from pyxsec.root
 293   //
 294
 295   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 296   Float_t xsection = 0;
 297   Float_t ftrials  = 1;
 298
 299   fAvgTrials = 1;
 300   if(tree){
 301     TFile *curfile = tree->GetCurrentFile();
 302     if (!curfile) {
 303       Error("Notify","No current file");
 304       return kFALSE;
 305     }
 306     if(!fh1Xsec||!fh1Trials){
 307       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 308       return kFALSE;
 309     }
 310     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 311     fh1Xsec->Fill("<#sigma>",xsection);
 312     // construct a poor man average trials
 313     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 314     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 315   }
 316   return kTRUE;
 317 }
 318
 319 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 320 {
 321
 322   //
 323   // Create the output container
 324   //
 325
 326
 327   // Connect the AOD
 328
 329
 330   MakeJetContainer();
 331
 332
 333   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 334
 335   OpenFile(1);
 336   if(!fHistList)fHistList = new TList();
 337   fHistList->SetOwner(kTRUE);
 338   Bool_t oldStatus = TH1::AddDirectoryStatus();
 339   TH1::AddDirectory(kFALSE);
 340
 341   //
 342   //  Histogram
 343
 344   const Int_t nBinPt = 320;
 345   Double_t binLimitsPt[nBinPt+1];
 346   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 347     if(iPt == 0){
 348       binLimitsPt[iPt] = 0.0;
 349     }
 350     else {// 1.0
 351       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 352     }
 353   }
 354
 355   const Int_t nBinPhi = 90;
 356   Double_t binLimitsPhi[nBinPhi+1];
 357   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 358     if(iPhi==0){
 359       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 360     }
 361     else{
 362       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 363     }
 364   }
 365
 366
 367   const Int_t nBinPhi2 = 360;
 368   Double_t binLimitsPhi2[nBinPhi2+1];
 369   for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
 370     if(iPhi2==0){
 371       binLimitsPhi2[iPhi2] = 0.;
 372     }
 373     else{
 374       binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
 375     }
 376   }
 377
 378
 379
 380   const Int_t nBinEta = 40;
 381   Double_t binLimitsEta[nBinEta+1];
 382   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 383     if(iEta==0){
 384       binLimitsEta[iEta] = -2.0;
 385     }
 386     else{
 387       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 388     }
 389   }
 390
 391   const Int_t nBinFrag = 25;
 392
 393
 394   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 395   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 396
 397   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 398   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 399
 400   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 401   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 402   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 403
 404   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 405
 406   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 407   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 408
 409   fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 410   fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
 411   fh1SumPtTrackAreaRec = new TH1F("fh1SumPtTrackAreaRec","Sum Rec track P_T #eta <0.9 / 1.8 * 2 * 0.4*0.4;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
 412
 413   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 414   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 415   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 416   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 417   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 418
 419   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 420   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 421   //
 422
 423   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 424                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 425   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 426                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 427   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 428                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 429   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 430                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 431
 432   fh2JetPtJetPhi = new TH2F("fh2JetPtJetPhi","Reconstructed jet phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
 433   fh2TrackPtTrackPhi = new TH2F("fh2TrackPtTrackPhi","Reconstructed track phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
 434
 435   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-1.2,1.2);
 436
 437   for(int ij = 0;ij < kMaxJets;++ij){
 438     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 439     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 440
 441     fh2PhiPt[ij] =  new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
 442                              nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 443
 444     fh2PhiEta[ij] =  new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
 445                               nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 446
 447     fh2RhoPtRec[ij] =  new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
 448                                 40,0.,1.,nBinPt,binLimitsPt);
 449     fh2PsiPtRec[ij] =  new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
 450                                 40,0.,2.,nBinPt,binLimitsPt);
 451
 452     fh2RhoPtGen[ij] =  new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
 453                                40,0.,2.,nBinPt,binLimitsPt);
 454     fh2PsiPtGen[ij] =  new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
 455                                 40,0.,2.,nBinPt,binLimitsPt);
 456     if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 457
 458
 459     fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
 460                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 461     fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 462                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 463
 464     fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
 465                               nBinFrag,0.,1.,nBinPt,binLimitsPt);
 466     fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
 467                                 nBinFrag,0.,10.,nBinPt,binLimitsPt);
 468   }
 469
 470   // Dijet histograms
 471
 472   fh2DijetDeltaPhiPt       = new TH2F("fh2DijetDeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
 473   fh2DijetAsymPt            = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
 474   fh2DijetAsymPtCut         = new TH2F("fh2DijetAsymCut","Pt asymmetry after delta phi cut;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
 475   fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
 476   fh2DijetPt2vsPt1          = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
 477   fh2DijetDifvsSum          = new TH2F("fh2DijetDifvsSum","Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
 478   fh1DijetMinv               = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 479   fh1DijetMinvCut           = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 480   //background histograms
 481   if(fBkgSubtraction){
 482     fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
 483     fh1Bkg2 = new TH1F("fh1Bkg2","Background estimate 2",100,0.,10.);
 484     fh1Bkg3 = new TH1F("fh1Bkg3","Background estimate 3",100,0.,10.);
 485     fh1Sigma1 = new TH1F("fh1Sigma1","Background fluctuations 1",100,0.,10.);
 486     fh1Sigma2 = new TH1F("fh1Sigma2","Background fluctuations 2",100,0.,10.);
 487     fh1Sigma3 = new TH1F("fh1Sigma3","Background fluctuations 3",100,0.,10.);
 488     fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
 489     fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
 490     fh1Area3 = new TH1F("fh1Area3","Background mean area 3",50,0.,5.);
 491     fh1Ptjet = new TH1F("fh1Ptjet","Jet spectrum",100,0.,200.);
 492     fh1Ptjetsub1 = new TH1F("fh1Ptjetsub1","Subtracted spectrum 1",50,0.,200.);
 493     fh1Ptjetsub2 = new TH1F("fh1Ptjetsub2","Subtracted spectrum 2",50,0.,200.);
 494     fh1Ptjetsub3 = new TH1F("fh1Ptjetsub3","Subtracted spectrum 3",50,0.,200.);
 495     fh1Ptjethardest = new TH1F("fh1Ptjethardest","Hardest jet spectrum",50,0.,200.);
 496     fh1Ptjetsubhardest1 = new TH1F("fh1Pthardestsub1","Subtracted hardest jet spectrum 1",100,0.,200.);
 497     fh1Ptjetsubhardest2 = new TH1F("fh1Pthardestsub2","Subtracted hardest jet spectrum 2",100,0.,200.);
 498     fh1Ptjetsubhardest3 = new TH1F("fh1Pthardestsub3","Subtracted hardest jet spectrum 3",100,0.,200.);
 499     fh2Rhovspthardest1 = new TH2F("fh2Rhovspthardest1","Background vs pTjet 1",100,0.,200.,50,0.,5.);
 500     fh2Rhovspthardest2 = new TH2F("fh2Rhovspthardest2","Background vs pTjet 2",100,0.,200.,50,0.,5.);
 501     fh2Rhovspthardest3 = new TH2F("fh2Rhovspthardest3","Background vs pTjet 3",100,0.,200.,50,0.,5.);
 502     fh2Errorvspthardest1 = new TH2F("fh2Errorvspthardest1","Relative error 1",100,0.,200.,50,0.,5.);
 503     fh2Errorvspthardest2 = new TH2F("fh2Errorvspthardest2","Relative error 2",100,0.,200.,50,0.,5.);
 504     fh2Errorvspthardest3 = new TH2F("fh2Errorvspthardest3","Relative error 3",100,0.,200.,50,0.,5.);
 505   }
 506
 507
 508
 509
 510   const Int_t saveLevel = 3; // large save level more histos
 511   if(saveLevel>0){
 512     fHistList->Add(fh1Xsec);
 513     fHistList->Add(fh1Trials);
 514     fHistList->Add(fh1PtHard);
 515     fHistList->Add(fh1PtHardNoW);
 516     fHistList->Add(fh1PtHardTrials);
 517     fHistList->Add(fh1ZVtx);
 518     if(fBranchGen.Length()>0||fBkgSubtraction){
 519       fHistList->Add(fh1NGenJets);
 520       fHistList->Add(fh1PtTracksGenIn);
 521     }
 522     fHistList->Add(fh1PtJetsRecIn);
 523     fHistList->Add(fh1PtJetsLeadingRecIn);
 524     fHistList->Add(fh1PtTracksRecIn);
 525     fHistList->Add(fh1PtTracksLeadingRecIn);
 526     fHistList->Add(fh1NRecJets);
 527     fHistList->Add(fh1PtTrackRec);
 528     fHistList->Add(fh1SumPtTrackRec);
 529     fHistList->Add(fh1SumPtTrackAreaRec);
 530     fHistList->Add(fh2NRecJetsPt);
 531     fHistList->Add(fh2NRecTracksPt);
 532     fHistList->Add(fh2JetsLeadingPhiEta );
 533     fHistList->Add(fh2JetsLeadingPhiPt );
 534     fHistList->Add(fh2TracksLeadingPhiEta);
 535     fHistList->Add(fh2TracksLeadingPhiPt);
 536     for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 537     for(int ij = 0;ij<kMaxJets;++ij){
 538       fHistList->Add( fh1PtRecIn[ij]);
 539
 540       if(fBranchGen.Length()>0||fBkgSubtraction){
 541         fHistList->Add(fh1PtGenIn[ij]);
 542         fHistList->Add(fh2FragGen[ij]);
 543         fHistList->Add(fh2FragLnGen[ij]);
 544         fHistList->Add(fh2RhoPtGen[ij]);
 545         fHistList->Add(fh2PsiPtGen[ij]);
 546       }
 547       fHistList->Add( fh2PhiPt[ij]);
 548       fHistList->Add( fh2PhiEta[ij]);
 549       fHistList->Add( fh2RhoPtRec[ij]);
 550       fHistList->Add( fh2PsiPtRec[ij]);
 551       fHistList->Add( fh2FragRec[ij]);
 552       fHistList->Add( fh2FragLnRec[ij]);
 553     }
 554     fHistList->Add(fhnCorrelation);
 555     fHistList->Add(fhnCorrelationPhiZRec);
 556     fHistList->Add(fh2JetPtJetPhi);
 557     fHistList->Add(fh2TrackPtTrackPhi);
 558     fHistList->Add(fh2RelPtFGen);
 559
 560     fHistList->Add(fh2DijetDeltaPhiPt);
 561     fHistList->Add(fh2DijetAsymPt);
 562     fHistList->Add(fh2DijetAsymPtCut);
 563     fHistList->Add(fh2DijetDeltaPhiDeltaEta);
 564     fHistList->Add(fh2DijetPt2vsPt1);
 565     fHistList->Add(fh2DijetDifvsSum);
 566     fHistList->Add(fh1DijetMinv);
 567     fHistList->Add(fh1DijetMinvCut);
 568     if(fBkgSubtraction){
 569       fHistList->Add(fh1Bkg1);
 570       fHistList->Add(fh1Bkg2);
 571       fHistList->Add(fh1Bkg3);
 572       fHistList->Add(fh1Sigma1);
 573       fHistList->Add(fh1Sigma2);
 574       fHistList->Add(fh1Sigma3);
 575       fHistList->Add(fh1Area1);
 576       fHistList->Add(fh1Area2);
 577       fHistList->Add(fh1Area3);
 578       fHistList->Add(fh1Ptjet);
 579       fHistList->Add(fh1Ptjethardest);
 580       fHistList->Add(fh1Ptjetsub1);
 581       fHistList->Add(fh1Ptjetsub2);
 582       fHistList->Add(fh1Ptjetsub3);
 583       fHistList->Add(fh1Ptjetsubhardest1);
 584       fHistList->Add(fh1Ptjetsubhardest2);
 585       fHistList->Add(fh1Ptjetsubhardest3);
 586       fHistList->Add(fh2Rhovspthardest1);
 587       fHistList->Add(fh2Rhovspthardest2);
 588       fHistList->Add(fh2Rhovspthardest3);
 589       fHistList->Add(fh2Errorvspthardest1);
 590       fHistList->Add(fh2Errorvspthardest2);
 591      fHistList->Add(fh2Errorvspthardest3);
 592     }
 593   }
 594
 595   // =========== Switch on Sumw2 for all histos ===========
 596   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 597     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 598     if (h1){
 599       h1->Sumw2();
 600       continue;
 601     }
 602     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 603     if(hn)hn->Sumw2();
 604   }
 605   TH1::AddDirectory(oldStatus);
 606 }
 607
 608 void AliAnalysisTaskJetSpectrum2::Init()
 609 {
 610   //
 611   // Initialization
 612   //
 613
 614   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 615
 616 }
 617
 618 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 619 {
 620
 621   Bool_t selected = kTRUE;
 622
 623   if(fUseGlobalSelection&&fEventSelectionMask==0){
 624     selected = AliAnalysisHelperJetTasks::Selected();
 625   }
 626   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 627     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 628   }
 629
 630   if(fEventClass>0){
 631     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 632   }
 633
 634   if(!selected){
 635     // no selection by the service task, we continue
 636     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 637     PostData(1, fHistList);
 638     return;
 639   }
 640
 641
 642   //
 643   // Execute analysis for current event
 644   //
 645   AliESDEvent *fESD = 0;
 646   if(fUseAODJetInput){
 647     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 648     if(!fAOD){
 649       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
 650       return;
 651     }
 652     // fethc the header
 653   }
 654   else{
 655     //  assume that the AOD is in the general output...
 656     fAOD  = AODEvent();
 657     if(!fAOD){
 658       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 659       return;
 660     }
 661     fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 662   }
 663
 664
 665
 666
 667   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 668
 669
 670   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 671
 672   if(!aodH){
 673     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 674     return;
 675   }
 676
 677   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 678   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 679   if(!aodRecJets){
 680     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 681     return;
 682   }
 683
 684   Int_t nJets = aodRecJets->GetEntriesFast();
 685
 686   // ==== General variables needed
 687   // We use statice array, not to fragment the memory
 688   AliAODJet genJets[kMaxJets];
 689   Int_t nGenJets = 0;
 690   AliAODJet recJets[kMaxJets];
 691   Int_t nRecJets = 0;
 692   ///////////////////////////
 693
 694
 695
 696
 697   if(fBkgSubtraction){
 698      AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
 699
 700
 701      if(!evBkg){
 702        Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
 703        return;
 704      }
 705
 706
 707      ///just to start: some very simple plots containing rho, sigma and area of the background.
 708
 709
 710      Float_t pthardest=0.;
 711
 712      if(nJets!=0){
 713
 714        Float_t bkg1=evBkg->GetBackground(0);
 715        Float_t bkg2=evBkg->GetBackground(1);
 716        Float_t bkg3=evBkg->GetBackground(2);
 717        Float_t sigma1=evBkg->GetSigma(0);
 718        Float_t sigma2=evBkg->GetSigma(1);
 719        Float_t sigma3=evBkg->GetSigma(2);
 720        Float_t area1=evBkg->GetMeanarea(0);
 721        Float_t area2=evBkg->GetMeanarea(1);
 722        Float_t area3=evBkg->GetMeanarea(2);
 723        fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
 724        fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
 725        fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
 726        fh1Sigma1->Fill(sigma1);
 727        fh1Sigma2->Fill(sigma2);
 728        fh1Sigma3->Fill(sigma3);
 729        fh1Area1->Fill(area1);
 730        fh1Area2->Fill(area2);
 731        fh1Area3->Fill(area3);
 732
 733        Int_t iSubJetCounter = 0;
 734        for(Int_t k=0;k<nJets;k++){
 735          AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
 736          fh1Ptjet->Fill(jet->Pt());
 737          Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
 738          Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
 739          Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
 740          if(ptsub2<0.) ptsub2=0.;
 741          if(ptsub1<0.) ptsub1=0.;
 742          if(ptsub3<0.) ptsub3=0.;
 743          Float_t err1=sigma1*sqrt(area1);
 744          Float_t err2=sigma2*sqrt(area2);
 745          Float_t err3=sigma3*sqrt(area3);
 746          fh1Ptjetsub1->Fill(ptsub1);
 747          fh1Ptjetsub2->Fill(ptsub2);
 748          fh1Ptjetsub3->Fill(ptsub3);
 749          if(k==0) {
 750            pthardest=jet->Pt();
 751            fh1Ptjethardest->Fill(pthardest);
 752            fh1Ptjetsubhardest1->Fill(ptsub1);
 753            fh1Ptjetsubhardest2->Fill(ptsub2);
 754            fh1Ptjetsubhardest3->Fill(ptsub3);
 755            if (ptsub1 > 0)
 756              fh2Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
 757            if (ptsub2 > 0)
 758              fh2Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
 759            if (ptsub3 > 0)
 760              fh2Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
 761          }
 762
 763          Float_t ptsub=0.;
 764          if(fFillCorrBkg==1) ptsub=ptsub1;
 765          if(fFillCorrBkg==2) ptsub=ptsub2;
 766          if(fFillCorrBkg==3) ptsub=ptsub3;
 767          if(ptsub>0){// avoid unphysical jets pT
 768            Float_t subphi=jet->Phi();
 769            Float_t subtheta=jet->Theta();
 770            Float_t subpz = ptsub/TMath::Tan(subtheta);
 771            Float_t subpx=ptsub*TMath::Cos(subphi);
 772            Float_t subpy=ptsub * TMath::Sin(subphi);
 773            Float_t subp  = TMath::Sqrt(ptsub*ptsub+subpz*subpz);
 774            if(k<kMaxJets){// only store the jets which are assoicated to anohter one
 775              genJets[iSubJetCounter].SetPxPyPzE(subpx,subpy,subpz,subp);
 776              iSubJetCounter++;
 777            }
 778          }
 779        }
 780        nGenJets = iSubJetCounter;
 781        fh2Rhovspthardest1->Fill(pthardest,bkg1);
 782        fh2Rhovspthardest2->Fill(pthardest,bkg2);
 783        fh2Rhovspthardest3->Fill(pthardest,bkg3);
 784      }
 785   }// background subtraction
 786
 787
 788   Double_t eventW = 1;
 789   Double_t ptHard = 0;
 790   Double_t nTrials = 1; // Trials for MC trigger
 791
 792   if(fUseExternalWeightOnly){
 793     eventW = fExternalWeight;
 794   }
 795
 796   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 797   //  if(fDebug>0)aodH->SetFillAOD(kFALSE);
 798   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 799   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 800     // this is the part we only use when we have MC information
 801     AliMCEvent* mcEvent = MCEvent();
 802     //    AliStack *pStack = 0;
 803     if(!mcEvent){
 804       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 805       return;
 806     }
 807     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 808     Int_t iCount = 0;
 809     if(pythiaGenHeader){
 810       nTrials = pythiaGenHeader->Trials();
 811       ptHard  = pythiaGenHeader->GetPtHard();
 812       int iProcessType = pythiaGenHeader->ProcessType();
 813       // 11 f+f -> f+f
 814       // 12 f+barf -> f+barf
 815       // 13 f+barf -> g+g
 816       // 28 f+g -> f+g
 817       // 53 g+g -> f+barf
 818       // 68 g+g -> g+g
 819       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 820       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 821
 822       // fetch the pythia generated jets only to be used here
 823       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 824       AliAODJet pythiaGenJets[kMaxJets];
 825       for(int ip = 0;ip < nPythiaGenJets;++ip){
 826         if(iCount>=kMaxJets)continue;
 827         Float_t p[4];
 828         pythiaGenHeader->TriggerJet(ip,p);
 829         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 830
 831         if(fBranchGen.Length()==0){
 832           /*
 833           if(fLimitGenJetEta){
 834             if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 835                pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 836           }
 837           */
 838           if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
 839           // if we have MC particles and we do not read from the aod branch
 840           // use the pythia jets
 841           if(!fBkgSubtraction){
 842             genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 843             iCount++;
 844           }
 845         }
 846       }
 847     }
 848     if(fBranchGen.Length()==0&&!fBkgSubtraction)nGenJets = iCount;
 849   }// (fAnalysisType&kMCESD)==kMCESD)
 850
 851
 852   // we simply fetch the tracks/mc particles as a list of AliVParticles
 853
 854   TList recParticles;
 855   TList genParticles;
 856
 857
 858
 859
 860   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 861   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 862   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 863   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 864
 865
 866   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 867   fh1PtHard->Fill(ptHard,eventW);
 868   fh1PtHardNoW->Fill(ptHard,1);
 869   fh1PtHardTrials->Fill(ptHard,nTrials);
 870   fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
 871
 872
 873   // If we set a second branch for the input jets fetch this
 874   if(fBranchGen.Length()>0&&!fBkgSubtraction){
 875     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 876     if(aodGenJets){
 877       Int_t iCount = 0;
 878       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 879         if(iCount>=kMaxJets)continue;
 880         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 881         if(!tmp)continue;
 882         /*
 883         if(fLimitGenJetEta){
 884           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 885              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 886         }
 887         */
 888         if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
 889         genJets[iCount] = *tmp;
 890         iCount++;
 891       }
 892       nGenJets = iCount;
 893     }
 894     else{
 895       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 896       if(fDebug>2)fAOD->Print();
 897     }
 898   }
 899
 900   fh1NGenJets->Fill(nGenJets);
 901   // We do not want to exceed the maximum jet number
 902   nGenJets = TMath::Min(nGenJets,kMaxJets);
 903
 904   // Fetch the reconstructed jets...
 905
 906   nRecJets = aodRecJets->GetEntries();
 907
 908   nRecJets = aodRecJets->GetEntries();
 909   fh1NRecJets->Fill(nRecJets);
 910
 911   // Do something with the all rec jets
 912   Int_t nRecOver = nRecJets;
 913
 914   // check that the jets are sorted
 915   Float_t ptOld = 999999;
 916   for(int ir = 0;ir < nRecJets;ir++){
 917     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
 918     Float_t tmpPt = tmp->Pt();
 919     if(tmpPt>ptOld){
 920       Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
 921     }
 922     ptOld = tmpPt;
 923   }
 924
 925
 926   if(nRecOver>0){
 927     TIterator *recIter = aodRecJets->MakeIterator();
 928     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
 929     Float_t pt = tmpRec->Pt();
 930     if(tmpRec){
 931       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 932         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 933         while(pt<ptCut&&tmpRec){
 934           nRecOver--;
 935           tmpRec = (AliAODJet*)(recIter->Next());
 936           if(tmpRec){
 937             pt = tmpRec->Pt();
 938           }
 939         }
 940         if(nRecOver<=0)break;
 941         fh2NRecJetsPt->Fill(ptCut,nRecOver);
 942       }
 943     }
 944     recIter->Reset();
 945     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
 946     Float_t phi = leading->Phi();
 947     if(phi<0)phi+=TMath::Pi()*2.;
 948     Float_t eta = leading->Eta();
 949     pt = leading->Pt();
 950     while((tmpRec = (AliAODJet*)(recIter->Next()))){
 951       Float_t tmpPt = tmpRec->Pt();
 952       fh1PtJetsRecIn->Fill(tmpPt);
 953       if(tmpRec==leading){
 954         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 955         continue;
 956       }
 957       // correlation
 958       Float_t tmpPhi =  tmpRec->Phi();
 959
 960       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 961       Float_t dPhi = phi - tmpRec->Phi();
 962       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 963       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 964       Float_t dEta = eta - tmpRec->Eta();
 965       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 966       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 967     }
 968     delete recIter;
 969   }
 970
 971   Int_t nTrackOver = recParticles.GetSize();
 972   // do the same for tracks and jets
 973   if(nTrackOver>0){
 974     TIterator *recIter = recParticles.MakeIterator();
 975     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 976     Float_t pt = tmpRec->Pt();
 977     //    Printf("Leading track p_t %3.3E",pt);
 978     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 979       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 980       while(pt<ptCut&&tmpRec){
 981         nTrackOver--;
 982         tmpRec = (AliVParticle*)(recIter->Next());
 983         if(tmpRec){
 984           pt = tmpRec->Pt();
 985         }
 986       }
 987       if(nTrackOver<=0)break;
 988       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 989     }
 990
 991     recIter->Reset();
 992     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 993     Float_t phi = leading->Phi();
 994     if(phi<0)phi+=TMath::Pi()*2.;
 995     Float_t eta = leading->Eta();
 996     pt  = leading->Pt();
 997     while((tmpRec = (AliVParticle*)(recIter->Next()))){
 998       Float_t tmpPt = tmpRec->Pt();
 999       fh1PtTracksRecIn->Fill(tmpPt);
1000       if(tmpRec==leading){
1001         fh1PtTracksLeadingRecIn->Fill(tmpPt);
1002         continue;
1003       }
1004       // correlation
1005       Float_t tmpPhi =  tmpRec->Phi();
1006
1007       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1008       Float_t dPhi = phi - tmpRec->Phi();
1009       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1010       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1011       Float_t dEta = eta - tmpRec->Eta();
1012       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1013       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1014     }
1015     delete recIter;
1016   }
1017
1018   if(genParticles.GetSize()){
1019     TIterator *genIter = genParticles.MakeIterator();
1020     AliVParticle *tmpGen = 0;
1021     while((tmpGen = (AliVParticle*)(genIter->Next()))){
1022       if(TMath::Abs(tmpGen->Eta())<0.9){
1023         Float_t tmpPt = tmpGen->Pt();
1024         fh1PtTracksGenIn->Fill(tmpPt);
1025       }
1026     }
1027     delete genIter;
1028   }
1029
1030   nRecJets = TMath::Min(nRecJets,kMaxJets);
1031   nJets = TMath::Min(nJets,kMaxJets);
1032   Int_t iCountRec = 0;
1033   for(int ir = 0;ir < nRecJets;++ir){
1034     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1035     if(!tmp)continue;
1036     if(tmp->Pt()<fMinJetPt)continue;
1037     recJets[ir] = *tmp;
1038     iCountRec++;
1039   }
1040   nRecJets = iCountRec;
1041
1042
1043   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1044   // Relate the jets
1045   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
1046   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
1047
1048
1049   for(int i = 0;i<kMaxJets;++i){
1050     iGenIndex[i] = iRecIndex[i] = -1;
1051   }
1052
1053
1054   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1055                                             iGenIndex,iRecIndex,fDebug);
1056
1057
1058   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1059
1060   if(fDebug){
1061     for(int i = 0;i<kMaxJets;++i){
1062       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1063       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1064     }
1065   }
1066
1067
1068
1069
1070   Double_t container[6];
1071   Double_t containerPhiZ[6];
1072
1073   // loop over generated jets
1074
1075   // radius; tmp, get from the jet header itself
1076   // at some point, how todeal woht FastJet on MC?
1077   Float_t radiusGen =0.4;
1078   Float_t radiusRec =0.4;
1079
1080   for(int ig = 0;ig < nGenJets;++ig){
1081     Double_t ptGen = genJets[ig].Pt();
1082     Double_t phiGen = genJets[ig].Phi();
1083     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1084     Double_t etaGen = genJets[ig].Eta();
1085
1086     container[3] = ptGen;
1087     container[4] = etaGen;
1088     container[5] = phiGen;
1089     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1090     Int_t ir = iRecIndex[ig];
1091
1092     if(TMath::Abs(etaGen)<fRecEtaWindow){
1093       fh1TmpRho->Reset();
1094
1095       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1096       fh1PtGenIn[ig]->Fill(ptGen,eventW);
1097       // fill the fragmentation function
1098       for(int it = 0;it<genParticles.GetEntries();++it){
1099         AliVParticle *part = (AliVParticle*)genParticles.At(it);
1100         Float_t deltaR = genJets[ig].DeltaR(part);
1101         if(ptGen>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1102         if(deltaR<radiusGen&&ptGen>0){
1103           Float_t z = part->Pt()/ptGen;
1104           Float_t lnz =  -1.*TMath::Log(z);
1105           fh2FragGen[ig]->Fill(z,ptGen,eventW);
1106           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1107         }
1108
1109       }
1110       Float_t rhoSum = 0;
1111       for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
1112         Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1113         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1114         rhoSum += rho;
1115         fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1116         fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
1117       }
1118     }
1119     if(ir>=0&&ir<nRecJets){
1120       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1121       Double_t etaRec = recJets[ir].Eta();
1122       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1123       if(TMath::Abs(etaRec)<fRecEtaWindow&&TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
1124     }
1125   }// loop over generated jets
1126
1127   Float_t sumPt = 0;
1128   for(int it = 0;it<recParticles.GetEntries();++it){
1129     AliVParticle *part = (AliVParticle*)recParticles.At(it);
1130     // fill sum pt and P_t of all paritles
1131     if(TMath::Abs(part->Eta())<0.9){
1132       Float_t pt = part->Pt();
1133       fh1PtTrackRec->Fill(pt,eventW);
1134       fh2TrackPtTrackPhi->Fill(pt,part->Phi());
1135       sumPt += pt;
1136     }
1137   }
1138   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1139   fh1SumPtTrackRec->Fill(sumPt,eventW);
1140
1141
1142   // loop over reconstructed jets
1143   for(int ir = 0;ir < nRecJets;++ir){
1144     Double_t etaRec = recJets[ir].Eta();
1145     Double_t ptRec = recJets[ir].Pt();
1146     Double_t phiRec = recJets[ir].Phi();
1147     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1148
1149
1150   // do something with dijets...
1151   if(ir==1){
1152     Double_t etaRec1 = recJets[0].Eta();
1153     Double_t ptRec1 = recJets[0].Pt();
1154     Double_t phiRec1 = recJets[0].Phi();
1155     if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
1156
1157
1158     if(TMath::Abs(etaRec1)<fRecEtaWindow
1159       &&TMath::Abs(etaRec)<fRecEtaWindow){
1160
1161       Float_t deltaPhi = phiRec1 - phiRec;
1162
1163       if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1164       if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1165       deltaPhi = TMath::Abs(deltaPhi);
1166       fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
1167       Float_t asym = 9999;
1168       if((ptRec1+ptRec)>0)asym = (ptRec1-ptRec)/(ptRec1+ptRec);
1169       fh2DijetAsymPt->Fill(asym,ptRec1);
1170       fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
1171       fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
1172       fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
1173       Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
1174                          recJets[0].Px()*recJets[1].Px()-
1175                          recJets[0].Py()*recJets[1].Py()-
1176                          recJets[0].Pz()*recJets[1].Pz());
1177       if(minv<0)minv=0; // prevent numerical instabilities
1178       minv = TMath::Sqrt(minv);
1179       // with mass == 0;
1180
1181       fh1DijetMinv->Fill(minv);
1182       if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
1183         fh1DijetMinvCut->Fill(minv);
1184         fh2DijetAsymPtCut->Fill(asym,ptRec1);
1185       }
1186     }
1187   }
1188
1189
1190     container[0] = ptRec;
1191     container[1] = etaRec;
1192     container[2] = phiRec;
1193     containerPhiZ[0] = ptRec;
1194     containerPhiZ[1] = phiRec;
1195     if(ptRec>30.&&fDebug>2){
1196       // need to cast to int, otherwise the printf overwrites
1197       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1198       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1199       if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
1200       //  aodH->SetFillAOD(kTRUE);
1201       fAOD->GetHeader()->Print();
1202       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
1203       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1204         AliAODTrack *tr = fAOD->GetTrack(it);
1205         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1206         tr->Print();
1207         //      tr->Dump();
1208         if(fESD){
1209           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1210           esdTr->Print("");
1211           // esdTr->Dump();
1212         }
1213       }
1214     }
1215
1216
1217     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1218     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1219
1220     Float_t zLeading = -1;
1221     if(TMath::Abs(etaRec)<fRecEtaWindow){
1222       fh2JetPtJetPhi->Fill(ptRec,phiRec);
1223       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1224       fh1PtRecIn[ir]->Fill(ptRec,eventW);
1225       // fill the fragmentation function
1226
1227       fh1TmpRho->Reset();
1228
1229       for(int it = 0;it<recParticles.GetEntries();++it){
1230         AliVParticle *part = (AliVParticle*)recParticles.At(it);
1231         Float_t eta = part->Eta();
1232         if(TMath::Abs(eta)<0.9){
1233           Float_t phi = part->Phi();
1234           if(phi<0)phi+=TMath::Pi()*2.;
1235           Float_t dPhi = phi - phiRec;
1236           Float_t dEta = eta - etaRec;
1237           if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1238           if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1239           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1240           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1241         }
1242
1243         Float_t deltaR = recJets[ir].DeltaR(part);
1244         if(ptRec>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1245
1246
1247         if(deltaR<radiusRec&&ptRec>0){
1248           Float_t z = part->Pt()/ptRec;
1249           if(z>zLeading)zLeading=z;
1250           Float_t lnz =  -1.*TMath::Log(z);
1251           fh2FragRec[ir]->Fill(z,ptRec,eventW);
1252           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1253         }
1254       }
1255       // fill the jet shapes
1256       Float_t rhoSum = 0;
1257       for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1258         Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1259         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1260         rhoSum += rho;
1261         fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1262         fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1263       }
1264     }
1265
1266     // Fill Correlation
1267     Int_t ig = iGenIndex[ir];
1268     if(ig>=0 && ig<nGenJets){
1269       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1270       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1271       Double_t ptGen  = genJets[ig].Pt();
1272       Double_t phiGen = genJets[ig].Phi();
1273       if(phiGen<0)phiGen+=TMath::Pi()*2.;
1274       Double_t etaGen = genJets[ig].Eta();
1275
1276       container[3] = ptGen;
1277       container[4] = etaGen;
1278       container[5] = phiGen;
1279       containerPhiZ[3] = ptGen;
1280       //
1281       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1282       //
1283
1284       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1285       if(TMath::Abs(etaRec)<fRecEtaWindow){
1286         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1287         fhnCorrelation->Fill(container);
1288         if(ptGen>0){
1289           Float_t delta = (ptRec-ptGen)/ptGen;
1290           fh2RelPtFGen->Fill(ptGen,delta,eventW);
1291         }
1292         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1293       }// if etarec in window
1294     }
1295     else{
1296       containerPhiZ[3] = 0;
1297       if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1298     }
1299   }// loop over reconstructed jets
1300
1301
1302   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1303   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1304   PostData(1, fHistList);
1305 }
1306
1307
1308 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1309   //
1310   // Create the particle container for the correction framework manager and
1311   // link it
1312   //
1313   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1314   const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1315   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1316   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1317   const Double_t kZmin = 0., kZmax = 1;
1318
1319   // can we neglect migration in eta and phi?
1320   // phi should be no problem since we cover full phi and are phi symmetric
1321   // eta migration is more difficult  due to needed acceptance correction
1322   // in limited eta range
1323
1324   //arrays for the number of bins in each dimension
1325   Int_t iBin[kNvar];
1326   iBin[0] = 320; //bins in pt
1327   iBin[1] =  1; //bins in eta
1328   iBin[2] = 1; // bins in phi
1329
1330   //arrays for lower bounds :
1331   Double_t* binEdges[kNvar];
1332   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1333     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1334
1335   //values for bin lower bounds
1336   //  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);
1337   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1338   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1339   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1340
1341
1342   for(int i = 0;i < kMaxStep*2;++i){
1343     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1344     for (int k=0; k<kNvar; k++) {
1345       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1346     }
1347   }
1348   //create correlation matrix for unfolding
1349   Int_t thnDim[2*kNvar];
1350   for (int k=0; k<kNvar; k++) {
1351     //first half  : reconstructed
1352     //second half : MC
1353     thnDim[k]      = iBin[k];
1354     thnDim[k+kNvar] = iBin[k];
1355   }
1356
1357   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1358   for (int k=0; k<kNvar; k++) {
1359     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1360     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1361   }
1362   fhnCorrelation->Sumw2();
1363
1364   // for second correlation histogram
1365
1366
1367   const Int_t kNvarPhiZ   = 4;
1368   //arrays for the number of bins in each dimension
1369   Int_t iBinPhiZ[kNvarPhiZ];
1370   iBinPhiZ[0] = 80; //bins in pt
1371   iBinPhiZ[1] = 72; //bins in phi
1372   iBinPhiZ[2] = 20; // bins in Z
1373   iBinPhiZ[3] = 80; //bins in ptgen
1374
1375   //arrays for lower bounds :
1376   Double_t* binEdgesPhiZ[kNvarPhiZ];
1377   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1378     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1379
1380   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1381   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1382   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1383   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1384
1385   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1386   for (int k=0; k<kNvarPhiZ; k++) {
1387     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1388   }
1389   fhnCorrelationPhiZRec->Sumw2();
1390
1391
1392   // Add a histogram for Fake jets
1393
1394   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1395     delete [] binEdges[ivar];
1396
1397   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1398     delete [] binEdgesPhiZ[ivar];
1399
1400 }
1401
1402 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1403 {
1404 // Terminate analysis
1405 //
1406     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1407 }
1408
1409
1410 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1411
1412   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1413
1414   Int_t iCount = 0;
1415   if(type==kTrackAOD){
1416     AliAODEvent *aod = 0;
1417     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1418     else aod = AODEvent();
1419     if(!aod){
1420       return iCount;
1421     }
1422     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1423       AliAODTrack *tr = aod->GetTrack(it);
1424       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1425       if(TMath::Abs(tr->Eta())>0.9)continue;
1426       if(fDebug>0){
1427         if(tr->Pt()>20){
1428           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1429           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1430           tr->Print();
1431           //    tr->Dump();
1432           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1433           if(fESD){
1434             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1435             esdTr->Print("");
1436             // esdTr->Dump();
1437           }
1438         }
1439       }
1440       list->Add(tr);
1441       iCount++;
1442     }
1443   }
1444   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1445     AliMCEvent* mcEvent = MCEvent();
1446     if(!mcEvent)return iCount;
1447     // we want to have alivpartilces so use get track
1448     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1449       if(!mcEvent->IsPhysicalPrimary(it))continue;
1450       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1451       if(type == kTrackKineAll){
1452         list->Add(part);
1453         iCount++;
1454       }
1455       else if(type == kTrackKineCharged){
1456         if(part->Particle()->GetPDG()->Charge()==0)continue;
1457         list->Add(part);
1458         iCount++;
1459       }
1460     }
1461   }
1462   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1463     AliAODEvent *aod = 0;
1464     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1465     else aod = AODEvent();
1466     if(!aod)return iCount;
1467     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1468     if(!tca)return iCount;
1469     for(int it = 0;it < tca->GetEntriesFast();++it){
1470       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1471       if(!part->IsPhysicalPrimary())continue;
1472       if(type == kTrackAODMCAll){
1473         list->Add(part);
1474         iCount++;
1475       }
1476       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1477         if(part->Charge()==0)continue;
1478         if(kTrackAODMCCharged){
1479           list->Add(part);
1480         }
1481         else {
1482           if(TMath::Abs(part->Eta())>0.9)continue;
1483           list->Add(part);
1484         }
1485         iCount++;
1486       }
1487     }
1488   }// AODMCparticle
1489   list->Sort();
1490   return iCount;
1491
1492 }