PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx

   1
   2 // **************************************
   3 // Task used for the correction of determiantion of reconstructed jet spectra
   4 // Compares input (gen) and output (rec) jets
   5 // *******************************************
   6
   7
   8 /**************************************************************************
   9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  10  *                                                                        *
  11  * Author: The ALICE Off-line Project.                                    *
  12  * Contributors are mentioned in the code where appropriate.              *
  13  *                                                                        *
  14  * Permission to use, copy, modify and distribute this software and its   *
  15  * documentation strictly for non-commercial purposes is hereby granted   *
  16  * without fee, provided that the above copyright notice appears in all   *
  17  * copies and that both the copyright notice and this permission notice   *
  18  * appear in the supporting documentation. The authors make no claims     *
  19  * about the suitability of this software for any purpose. It is          *
  20  * provided "as is" without express or implied warranty.                  *
  21  **************************************************************************/
  22
  23
  24 #include <TROOT.h>
  25 #include <TRandom.h>
  26 #include <TSystem.h>
  27 #include <TInterpreter.h>
  28 #include <TChain.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetSpectrum2.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliJetFinder.h"
  43 #include "AliJetHeader.h"
  44 #include "AliJetReader.h"
  45 #include "AliJetReaderHeader.h"
  46 #include "AliUA1JetHeaderV1.h"
  47 #include "AliESDEvent.h"
  48 #include "AliAODEvent.h"
  49 #include "AliAODHandler.h"
  50 #include "AliAODTrack.h"
  51 #include "AliAODJet.h"
  52 #include "AliAODJetEventBackground.h"
  53 #include "AliAODMCParticle.h"
  54 #include "AliMCEventHandler.h"
  55 #include "AliMCEvent.h"
  56 #include "AliStack.h"
  57 #include "AliGenPythiaEventHeader.h"
  58 #include "AliJetKineReaderHeader.h"
  59 #include "AliGenCocktailEventHeader.h"
  60 #include "AliInputEventHandler.h"
  61
  62
  63 #include "AliAnalysisHelperJetTasks.h"
  64
  65 ClassImp(AliAnalysisTaskJetSpectrum2)
  66
  67 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
  68   fJetHeaderRec(0x0),
  69   fJetHeaderGen(0x0),
  70   fAODIn(0x0),
  71   fAODOut(0x0),
  72   fAODExtension(0x0),
  73   fhnCorrelation(0x0),
  74   fhnCorrelationPhiZRec(0x0),
  75   f1PtScale(0x0),
  76   fBranchRec("jets"),
  77   fBranchGen(""),
  78   fBranchBkg(""),
  79   fNonStdFile(""),
  80   fUseAODJetInput(kFALSE),
  81   fUseAODTrackInput(kFALSE),
  82   fUseAODMCInput(kFALSE),
  83   fUseGlobalSelection(kFALSE),
  84   fUseExternalWeightOnly(kFALSE),
  85   fLimitGenJetEta(kFALSE),
  86   fBkgSubtraction(kFALSE),
  87   fUseOldFill(kFALSE),
  88   fNMatchJets(5),
  89   fFillCorrBkg(0),
  90   fFilterMask(0),
  91   fEventSelectionMask(0),
  92   fAnalysisType(0),
  93   fTrackTypeRec(kTrackUndef),
  94   fTrackTypeGen(kTrackUndef),
  95   fEventClass(0),
  96   fAvgTrials(1),
  97   fExternalWeight(1),
  98   fJetRecEtaWindow(0.5),
  99   fTrackRecEtaWindow(0.5),
 100   fMinJetPt(0),
 101   fMinTrackPt(0.15),
 102   fDeltaPhiWindow(90./180.*TMath::Pi()),
 103   fh1Xsec(0x0),
 104   fh1Trials(0x0),
 105   fh1PtHard(0x0),
 106   fh1PtHardNoW(0x0),
 107   fh1PtHardTrials(0x0),
 108   fh1ZVtx(0x0),
 109   fh1NGenJets(0x0),
 110   fh1NRecJets(0x0),
 111   fh1PtTrackRec(0x0),
 112   fh1SumPtTrackRec(0x0),
 113   fh1SumPtTrackAreaRec(0x0),
 114   fh1TmpRho(0x0),
 115   fh1PtJetsRecIn(0x0),
 116   fh1PtJetsGenIn(0x0),
 117   fh1PtJetsLeadingRecIn(0x0),
 118   fh1PtTracksRecIn(0x0),
 119   fh1PtTracksLeadingRecIn(0x0),
 120   fh1PtTracksGenIn(0x0),
 121   fh2NRecJetsPt(0x0),
 122   fh2NRecTracksPt(0x0),
 123   fh2NGenJetsPt(0x0),
 124   fh2NGenTracksPt(0x0),
 125   fh2PtFGen(0x0),
 126   fh2RelPtFGen(0x0),
 127   fh1Bkg1(0x0),
 128   fh1Bkg2(0x0),
 129   fh1Bkg3(0x0),
 130   fh1Sigma1(0x0),
 131   fh1Sigma2(0x0),
 132   fh1Sigma3(0x0),
 133   fh1Area1(0x0),
 134   fh1Area2(0x0),
 135   fh1Area3(0x0),
 136   fh1Ptjet(0x0),
 137   fh1Ptjetsub1(0x0),
 138   fh1Ptjetsub2(0x0),
 139   fh1Ptjetsub3(0x0),
 140   fh1Ptjethardest(0x0),
 141   fh1Ptjetsubhardest1(0x0),
 142   fh1Ptjetsubhardest2(0x0),
 143   fh1Ptjetsubhardest3(0x0),
 144   fh2Rhovspthardest1(0x0),
 145   fh2Rhovspthardest2(0x0),
 146   fh2Rhovspthardest3(0x0),
 147   fh2Errorvspthardest1(0x0),
 148   fh2Errorvspthardest2(0x0),
 149   fh2Errorvspthardest3(0x0),
 150   fHistList(0x0)
 151 {
 152   for(int i = 0;i < kMaxStep*2;++i){
 153     fhnJetContainer[i] = 0;
 154   }
 155   for(int i = 0;i < kMaxJets;++i){
 156     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 157
 158     fh2PhiPt[i] = 0;
 159     fh2PhiEta[i] = 0;
 160     fh2RhoPtRec[i] = 0;
 161     fh2RhoPtGen[i] = 0;
 162     fh2PsiPtGen[i] = 0;
 163     fh2PsiPtRec[i] = 0;
 164     fh2FragRec[i] = 0;
 165     fh2FragLnRec[i] = 0;
 166     fh2FragGen[i] = 0;
 167     fh2FragLnGen[i] = 0;
 168   }
 169
 170
 171   for(int ij = 0;ij <kJetTypes;++ij){
 172     fh1NJets[ij] = 0;
 173     fh1SumPtTrack[ij] = 0;
 174     fh1PtJetsIn[ij] = 0;
 175     fh1PtTracksIn[ij] = 0;
 176     fh1PtTracksLeadingIn[ij] = 0;
 177     fh2NJetsPt[ij]  = 0;
 178     fh2NTracksPt[ij]  = 0;
 179     fh2LeadingJetPtJetPhi[ij] = 0;
 180     fh2LeadingTrackPtTrackPhi[ij] = 0;
 181     for(int i = 0;i < kMaxJets;++i){
 182       fh1PtIn[ij][i] = 0;
 183       fh2RhoPt[ij][i] = 0;
 184       fh2PsiPt[ij][i] = 0;
 185     }
 186
 187     fh1DijetMinv[ij] = 0;
 188     fh2DijetDeltaPhiPt[ij] = 0;
 189     fh2DijetAsymPt[ij] = 0;
 190     fh2DijetPt2vsPt1[ij] = 0;
 191     fh2DijetDifvsSum[ij] = 0;
 192
 193   }
 194 }
 195
 196 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 197   AliAnalysisTaskSE(name),
 198   fJetHeaderRec(0x0),
 199   fJetHeaderGen(0x0),
 200   fAODIn(0x0),
 201   fAODOut(0x0),
 202   fAODExtension(0x0),
 203   fhnCorrelation(0x0),
 204   fhnCorrelationPhiZRec(0x0),
 205   f1PtScale(0x0),
 206   fBranchRec("jets"),
 207   fBranchGen(""),
 208   fBranchBkg(""),
 209   fNonStdFile(""),
 210   fUseAODJetInput(kFALSE),
 211   fUseAODTrackInput(kFALSE),
 212   fUseAODMCInput(kFALSE),
 213   fUseGlobalSelection(kFALSE),
 214   fUseExternalWeightOnly(kFALSE),
 215   fLimitGenJetEta(kFALSE),
 216   fBkgSubtraction(kFALSE),
 217   fUseOldFill(kFALSE),
 218   fNMatchJets(5),
 219   fFillCorrBkg(0),
 220   fFilterMask(0),
 221   fEventSelectionMask(0),
 222   fAnalysisType(0),
 223   fTrackTypeRec(kTrackUndef),
 224   fTrackTypeGen(kTrackUndef),
 225   fEventClass(0),
 226   fAvgTrials(1),
 227   fExternalWeight(1),
 228   fJetRecEtaWindow(0.5),
 229   fTrackRecEtaWindow(0.5),
 230   fMinJetPt(0),
 231   fMinTrackPt(0.15),
 232   fDeltaPhiWindow(90./180.*TMath::Pi()),
 233   fh1Xsec(0x0),
 234   fh1Trials(0x0),
 235   fh1PtHard(0x0),
 236   fh1PtHardNoW(0x0),
 237   fh1PtHardTrials(0x0),
 238   fh1ZVtx(0x0),
 239   fh1NGenJets(0x0),
 240   fh1NRecJets(0x0),
 241   fh1PtTrackRec(0x0),
 242   fh1SumPtTrackRec(0x0),
 243   fh1SumPtTrackAreaRec(0x0),
 244   fh1TmpRho(0x0),
 245   fh1PtJetsRecIn(0x0),
 246   fh1PtJetsGenIn(0x0),
 247   fh1PtJetsLeadingRecIn(0x0),
 248   fh1PtTracksRecIn(0x0),
 249   fh1PtTracksLeadingRecIn(0x0),
 250   fh1PtTracksGenIn(0x0),
 251   fh2NRecJetsPt(0x0),
 252   fh2NRecTracksPt(0x0),
 253   fh2NGenJetsPt(0x0),
 254   fh2NGenTracksPt(0x0),
 255   fh2PtFGen(0x0),
 256   fh2RelPtFGen(0x0),
 257   fh1Bkg1(0x0),
 258   fh1Bkg2(0x0),
 259   fh1Bkg3(0x0),
 260   fh1Sigma1(0x0),
 261   fh1Sigma2(0x0),
 262   fh1Sigma3(0x0),
 263   fh1Area1(0x0),
 264   fh1Area2(0x0),
 265   fh1Area3(0x0),
 266   fh1Ptjet(0x0),
 267   fh1Ptjetsub1(0x0),
 268   fh1Ptjetsub2(0x0),
 269   fh1Ptjetsub3(0x0),
 270   fh1Ptjethardest(0x0),
 271   fh1Ptjetsubhardest1(0x0),
 272   fh1Ptjetsubhardest2(0x0),
 273   fh1Ptjetsubhardest3(0x0),
 274   fh2Rhovspthardest1(0x0),
 275   fh2Rhovspthardest2(0x0),
 276   fh2Rhovspthardest3(0x0),
 277   fh2Errorvspthardest1(0x0),
 278   fh2Errorvspthardest2(0x0),
 279   fh2Errorvspthardest3(0x0),
 280   fHistList(0x0)
 281 {
 282
 283   for(int i = 0;i < kMaxStep*2;++i){
 284     fhnJetContainer[i] = 0;
 285   }
 286   for(int i = 0;i < kMaxJets;++i){
 287     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 288
 289     fh2PhiPt[i] = 0;
 290     fh2PhiEta[i] = 0;
 291     fh2RhoPtRec[i] = 0;
 292     fh2RhoPtGen[i] = 0;
 293     fh2PsiPtGen[i] = 0;
 294     fh2PsiPtRec[i] = 0;
 295     fh2FragRec[i] = 0;
 296     fh2FragLnRec[i] = 0;
 297     fh2FragGen[i] = 0;
 298     fh2FragLnGen[i] = 0;
 299   }
 300
 301   for(int ij = 0;ij <kJetTypes;++ij){
 302     fh1NJets[ij] = 0;
 303     fh1SumPtTrack[ij] = 0;
 304     fh1PtJetsIn[ij] = 0;
 305     fh1PtTracksIn[ij] = 0;
 306     fh1PtTracksLeadingIn[ij] = 0;
 307     fh2NJetsPt[ij]  = 0;
 308     fh2NTracksPt[ij]  = 0;
 309     fh2LeadingJetPtJetPhi[ij] = 0;
 310     fh2LeadingTrackPtTrackPhi[ij] = 0;
 311     for(int i = 0;i < kMaxJets;++i){
 312       fh1PtIn[ij][i] = 0;
 313       fh2RhoPt[ij][i] = 0;
 314       fh2PsiPt[ij][i] = 0;
 315     }
 316
 317     fh1DijetMinv[ij] = 0;
 318     fh2DijetDeltaPhiPt[ij] = 0;
 319     fh2DijetAsymPt[ij] = 0;
 320     fh2DijetPt2vsPt1[ij] = 0;
 321     fh2DijetDifvsSum[ij] = 0;
 322   }
 323
 324  DefineOutput(1, TList::Class());
 325 }
 326
 327
 328
 329 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 330 {
 331
 332
 333
 334   //
 335   // Implemented Notify() to read the cross sections
 336   // and number of trials from pyxsec.root
 337   //
 338
 339   // Fetch the aod also from the input in,
 340   // have todo it in notify
 341
 342
 343   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 344   //  assume that the AOD is in the general output...
 345   fAODOut  = AODEvent();
 346
 347   if(fNonStdFile.Length()!=0){
 348     // case that we have an AOD extension we need can fetch the jets from the extended output
 349     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 350     fAODExtension = aodH->GetExtension(fNonStdFile.Data());
 351     if(!fAODExtension){
 352       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 353     }
 354   }
 355
 356
 357   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 358   Float_t xsection = 0;
 359   Float_t ftrials  = 1;
 360
 361   fAvgTrials = 1;
 362   if(tree){
 363     TFile *curfile = tree->GetCurrentFile();
 364     if (!curfile) {
 365       Error("Notify","No current file");
 366       return kFALSE;
 367     }
 368     if(!fh1Xsec||!fh1Trials){
 369       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 370       return kFALSE;
 371     }
 372     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 373     fh1Xsec->Fill("<#sigma>",xsection);
 374     // construct a poor man average trials
 375     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 376     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 377   }
 378
 379
 380   return kTRUE;
 381 }
 382
 383 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 384 {
 385
 386
 387   // Connect the AOD
 388
 389
 390
 391
 392
 393   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 394
 395   OpenFile(1);
 396   if(!fHistList)fHistList = new TList();
 397   fHistList->SetOwner(kTRUE);
 398   Bool_t oldStatus = TH1::AddDirectoryStatus();
 399   TH1::AddDirectory(kFALSE);
 400
 401
 402   MakeJetContainer();
 403   fHistList->Add(fhnCorrelation);
 404   fHistList->Add(fhnCorrelationPhiZRec);
 405   for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 406
 407   //
 408   //  Histogram
 409
 410   const Int_t nBinPt = 320;
 411   Double_t binLimitsPt[nBinPt+1];
 412   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 413     if(iPt == 0){
 414       binLimitsPt[iPt] = 0.0;
 415     }
 416     else {// 1.0
 417       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 418     }
 419   }
 420
 421   const Int_t nBinPhi = 90;
 422   Double_t binLimitsPhi[nBinPhi+1];
 423   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 424     if(iPhi==0){
 425       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 426     }
 427     else{
 428       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 429     }
 430   }
 431
 432
 433   const Int_t nBinPhi2 = 360;
 434   Double_t binLimitsPhi2[nBinPhi2+1];
 435   for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
 436     if(iPhi2==0){
 437       binLimitsPhi2[iPhi2] = 0.;
 438     }
 439     else{
 440       binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
 441     }
 442   }
 443
 444   const Int_t nBinEta = 40;
 445   Double_t binLimitsEta[nBinEta+1];
 446   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 447     if(iEta==0){
 448       binLimitsEta[iEta] = -2.0;
 449     }
 450     else{
 451       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 452     }
 453   }
 454
 455   const Int_t nBinFrag = 25;
 456
 457
 458
 459
 460   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 461   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 462    fHistList->Add(fh1Xsec);
 463
 464   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 465   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 466   fHistList->Add(fh1Trials);
 467
 468   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 469   fHistList->Add(fh1PtHard);
 470
 471   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 472   fHistList->Add(fh1PtHardNoW);
 473
 474   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 475   fHistList->Add(fh1PtHardTrials);
 476
 477
 478   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 479   fHistList->Add(fh1ZVtx);
 480
 481   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 482   fHistList->Add(fh2PtFGen);
 483
 484   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-2.4,2.4);
 485   fHistList->Add(fh2RelPtFGen);
 486
 487
 488   if(fUseOldFill){
 489
 490     fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 491     fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
 492
 493   fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 494   fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
 495   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);
 496
 497   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 498   fh1PtJetsGenIn  = new TH1F("fh1PtJetsGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 499   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 500   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 501   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 502   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 503
 504   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 505   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 506
 507   fh2NGenJetsPt = new TH2F("fh2NGenJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 508   fh2NGenTracksPt = new TH2F("fh2NGenTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 509   //
 510
 511   for(int ij = 0;ij < kMaxJets;++ij){
 512     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
 513     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
 514
 515     fh2PhiPt[ij] =  new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
 516                              nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 517
 518     fh2PhiEta[ij] =  new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
 519                               nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 520
 521     fh2RhoPtRec[ij] =  new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
 522                                 40,0.,1.,nBinPt,binLimitsPt);
 523     fh2PsiPtRec[ij] =  new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
 524                                 40,0.,2.,nBinPt,binLimitsPt);
 525
 526     fh2RhoPtGen[ij] =  new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
 527                                40,0.,2.,nBinPt,binLimitsPt);
 528     fh2PsiPtGen[ij] =  new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
 529                                 40,0.,2.,nBinPt,binLimitsPt);
 530     if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 531
 532     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",
 533                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
 534     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",
 535                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
 536
 537     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",
 538                               nBinFrag,0.,1.,nBinPt,binLimitsPt);
 539     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",
 540                                 nBinFrag,0.,10.,nBinPt,binLimitsPt);
 541   }
 542
 543   // Dijet histograms
 544   //background histograms
 545   if(fBkgSubtraction){
 546     fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
 547     fh1Bkg2 = new TH1F("fh1Bkg2","Background estimate 2",100,0.,10.);
 548     fh1Bkg3 = new TH1F("fh1Bkg3","Background estimate 3",100,0.,10.);
 549     fh1Sigma1 = new TH1F("fh1Sigma1","Background fluctuations 1",100,0.,10.);
 550     fh1Sigma2 = new TH1F("fh1Sigma2","Background fluctuations 2",100,0.,10.);
 551     fh1Sigma3 = new TH1F("fh1Sigma3","Background fluctuations 3",100,0.,10.);
 552     fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
 553     fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
 554     fh1Area3 = new TH1F("fh1Area3","Background mean area 3",50,0.,5.);
 555     fh1Ptjet = new TH1F("fh1Ptjet","Jet spectrum",100,0.,200.);
 556     fh1Ptjetsub1 = new TH1F("fh1Ptjetsub1","Subtracted spectrum 1",50,0.,200.);
 557     fh1Ptjetsub2 = new TH1F("fh1Ptjetsub2","Subtracted spectrum 2",50,0.,200.);
 558     fh1Ptjetsub3 = new TH1F("fh1Ptjetsub3","Subtracted spectrum 3",50,0.,200.);
 559     fh1Ptjethardest = new TH1F("fh1Ptjethardest","Hardest jet spectrum",50,0.,200.);
 560     fh1Ptjetsubhardest1 = new TH1F("fh1Pthardestsub1","Subtracted hardest jet spectrum 1",100,0.,200.);
 561     fh1Ptjetsubhardest2 = new TH1F("fh1Pthardestsub2","Subtracted hardest jet spectrum 2",100,0.,200.);
 562     fh1Ptjetsubhardest3 = new TH1F("fh1Pthardestsub3","Subtracted hardest jet spectrum 3",100,0.,200.);
 563     fh2Rhovspthardest1 = new TH2F("fh2Rhovspthardest1","Background vs pTjet 1",100,0.,200.,50,0.,5.);
 564     fh2Rhovspthardest2 = new TH2F("fh2Rhovspthardest2","Background vs pTjet 2",100,0.,200.,50,0.,5.);
 565     fh2Rhovspthardest3 = new TH2F("fh2Rhovspthardest3","Background vs pTjet 3",100,0.,200.,50,0.,5.);
 566     fh2Errorvspthardest1 = new TH2F("fh2Errorvspthardest1","Relative error 1",100,0.,200.,50,0.,5.);
 567     fh2Errorvspthardest2 = new TH2F("fh2Errorvspthardest2","Relative error 2",100,0.,200.,50,0.,5.);
 568     fh2Errorvspthardest3 = new TH2F("fh2Errorvspthardest3","Relative error 3",100,0.,200.,50,0.,5.);
 569   }
 570
 571
 572   const Int_t saveLevel = 3; // large save level more histos
 573   if(saveLevel>0){
 574      if(fBranchGen.Length()>0||fBkgSubtraction){
 575        fHistList->Add(fh1NGenJets);
 576        fHistList->Add(fh1PtTracksGenIn);
 577     }
 578     fHistList->Add(fh1PtJetsRecIn);
 579     fHistList->Add(fh1PtJetsGenIn);
 580     fHistList->Add(fh1PtJetsLeadingRecIn);
 581     fHistList->Add(fh1PtTracksRecIn);
 582     fHistList->Add(fh1PtTracksLeadingRecIn);
 583     fHistList->Add(fh1NRecJets);
 584     fHistList->Add(fh1PtTrackRec);
 585     fHistList->Add(fh1SumPtTrackRec);
 586     fHistList->Add(fh1SumPtTrackAreaRec);
 587     fHistList->Add(fh2NRecJetsPt);
 588     fHistList->Add(fh2NRecTracksPt);
 589     fHistList->Add(fh2NGenJetsPt);
 590     fHistList->Add(fh2NGenTracksPt);
 591     for(int ij = 0;ij<kMaxJets;++ij){
 592       fHistList->Add( fh1PtRecIn[ij]);
 593
 594       if(fBranchGen.Length()>0||fBkgSubtraction){
 595         fHistList->Add(fh1PtGenIn[ij]);
 596         fHistList->Add(fh2FragGen[ij]);
 597         fHistList->Add(fh2FragLnGen[ij]);
 598         fHistList->Add(fh2RhoPtGen[ij]);
 599         fHistList->Add(fh2PsiPtGen[ij]);
 600       }
 601       fHistList->Add( fh2PhiPt[ij]);
 602       fHistList->Add( fh2PhiEta[ij]);
 603       fHistList->Add( fh2RhoPtRec[ij]);
 604       fHistList->Add( fh2PsiPtRec[ij]);
 605       fHistList->Add( fh2FragRec[ij]);
 606       fHistList->Add( fh2FragLnRec[ij]);
 607     }
 608     //    fHistList->Add(fh2TrackPtTrackPhi);
 609     if(fBkgSubtraction){
 610       fHistList->Add(fh1Bkg1);
 611       fHistList->Add(fh1Bkg2);
 612       fHistList->Add(fh1Bkg3);
 613       fHistList->Add(fh1Sigma1);
 614       fHistList->Add(fh1Sigma2);
 615       fHistList->Add(fh1Sigma3);
 616       fHistList->Add(fh1Area1);
 617       fHistList->Add(fh1Area2);
 618       fHistList->Add(fh1Area3);
 619       fHistList->Add(fh1Ptjet);
 620       fHistList->Add(fh1Ptjethardest);
 621       fHistList->Add(fh1Ptjetsub1);
 622       fHistList->Add(fh1Ptjetsub2);
 623       fHistList->Add(fh1Ptjetsub3);
 624       fHistList->Add(fh1Ptjetsubhardest1);
 625       fHistList->Add(fh1Ptjetsubhardest2);
 626       fHistList->Add(fh1Ptjetsubhardest3);
 627       fHistList->Add(fh2Rhovspthardest1);
 628       fHistList->Add(fh2Rhovspthardest2);
 629       fHistList->Add(fh2Rhovspthardest3);
 630       fHistList->Add(fh2Errorvspthardest1);
 631       fHistList->Add(fh2Errorvspthardest2);
 632      fHistList->Add(fh2Errorvspthardest3);
 633     }
 634   }
 635
 636
 637   }// ue old fill
 638   else{
 639
 640     for(int ij = 0;ij <kJetTypes;++ij){
 641       TString cAdd = "";
 642       TString cJetBranch = "";
 643       if(ij==kJetRec){
 644         cAdd = "Rec";
 645         cJetBranch = fBranchRec.Data();
 646       }
 647       else if (ij==kJetGen){
 648         cAdd = "Gen";
 649         cJetBranch = fBranchGen.Data();
 650       }
 651
 652       fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
 653       fHistList->Add(fh1NJets[ij]);
 654
 655       fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 656       fHistList->Add(fh1PtJetsIn[ij]);
 657
 658       fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 659       fHistList->Add(fh1PtTracksIn[ij]);
 660
 661       fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 662       fHistList->Add(fh1PtTracksLeadingIn[ij]);
 663
 664       fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 665       fHistList->Add(fh1SumPtTrack[ij]);
 666
 667       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);
 668       fHistList->Add(fh2NJetsPt[ij]);
 669
 670       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,50,-0.5,49.5);
 671       fHistList->Add(fh2NTracksPt[ij]);
 672
 673       fh2LeadingJetPtJetPhi[ij] = new TH2F(Form("fh2Leading%sJetPtJetPhi",cAdd.Data()),Form("phi of leading %s jet;p_{T};#phi;",cAdd.Data()),
 674                                            nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
 675       fHistList->Add(fh2LeadingJetPtJetPhi[ij]);
 676
 677       fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
 678                                                nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
 679       fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
 680
 681       for(int i = 0;i < kMaxJets;++i){
 682
 683         fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
 684         fHistList->Add(fh1PtIn[ij][i]);
 685
 686         fh2RhoPt[ij][i] = new TH2F(Form("fh2RhoPt%s_j%d",cAdd.Data(),i),Form("jet shape rho for %s jets;r;p_{T};",cAdd.Data()),
 687                                    40,0.,2.,nBinPt,binLimitsPt);
 688         fHistList->Add(fh2RhoPt[ij][i]);
 689         if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 690
 691
 692         fh2PsiPt[ij][i] = new TH2F(Form("fh2PsiPt%s_j%d",cAdd.Data(),i),Form("jet shape #psi for %s jets;r;p_{T};",cAdd.Data()),
 693                                      40,0.,2.,nBinPt,binLimitsPt);
 694         fHistList->Add(fh2PsiPt[ij][i]);
 695       }
 696
 697
 698       fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 699       fHistList->Add(fh1DijetMinv[ij]);
 700
 701       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);
 702       fHistList->Add(fh2DijetDeltaPhiPt[ij]);
 703
 704       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);
 705       fHistList->Add(fh2DijetAsymPt[ij]);
 706
 707       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.);
 708       fHistList->Add(fh2DijetPt2vsPt1[ij]);
 709
 710       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.);
 711       fHistList->Add( fh2DijetDifvsSum[ij]);
 712     }
 713   }// new filling
 714
 715   // =========== Switch on Sumw2 for all histos ===========
 716   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 717     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 718     if (h1){
 719       h1->Sumw2();
 720       continue;
 721     }
 722     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 723     if(hn)hn->Sumw2();
 724   }
 725   TH1::AddDirectory(oldStatus);
 726 }
 727
 728 void AliAnalysisTaskJetSpectrum2::Init()
 729 {
 730   //
 731   // Initialization
 732   //
 733
 734   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 735
 736 }
 737
 738 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
 739
 740
 741   Bool_t selected = kTRUE;
 742
 743   if(fUseGlobalSelection&&fEventSelectionMask==0){
 744     selected = AliAnalysisHelperJetTasks::Selected();
 745   }
 746   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 747     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 748   }
 749
 750   if(fEventClass>0){
 751     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 752   }
 753
 754   if(!selected){
 755     // no selection by the service task, we continue
 756     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 757     PostData(1, fHistList);
 758     return;
 759   }
 760
 761
 762   static AliAODEvent* aod = 0;
 763
 764   // take all other information from the aod we take the tracks from
 765   if(!aod){
 766    if(fUseAODTrackInput)aod = fAODIn;
 767    else aod = fAODOut;
 768   }
 769
 770
 771   //
 772   // Execute analysis for current event
 773   //
 774   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 775   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 776   if(!aodH){
 777     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 778     return;
 779   }
 780
 781   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 782   TClonesArray *aodRecJets = 0;
 783
 784   if(fAODOut&&!aodRecJets){
 785     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
 786   }
 787   if(fAODExtension&&!aodRecJets){
 788     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
 789   }
 790   if(fAODIn&&!aodRecJets){
 791     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
 792   }
 793
 794
 795
 796   if(!aodRecJets){
 797     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 798     return;
 799   }
 800
 801   TClonesArray *aodGenJets = 0;
 802   if(fBranchGen.Length()>0){
 803     if(fAODOut&&!aodGenJets){
 804       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
 805     }
 806     if(fAODExtension&&!aodGenJets){
 807       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
 808     }
 809     if(fAODIn&&!aodGenJets){
 810       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
 811     }
 812
 813     if(!aodGenJets){
 814       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
 815       return;
 816     }
 817   }
 818
 819
 820   if(false){
 821     UserExecOld();
 822     return;
 823   }
 824
 825
 826   // new Scheme
 827   // first fill all the pure  histograms, i.e. full jets
 828   // in case of the correaltion limit to the n-leading jets
 829
 830   // reconstructed
 831
 832
 833   // generated
 834
 835
 836   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
 837
 838
 839   TList genJetsList;         // full acceptance
 840   TList genJetsListCut;      // acceptance cut
 841   TList recJetsList;         // full acceptance
 842   TList recJetsListCut;      // acceptance cut
 843
 844
 845   GetListOfJets(&recJetsList,aodRecJets,0);
 846   GetListOfJets(&recJetsListCut,aodRecJets,1);
 847
 848   if(fBranchGen.Length()>0){
 849     GetListOfJets(&genJetsList,aodGenJets,0);
 850     GetListOfJets(&genJetsListCut,aodGenJets,1);
 851   }
 852
 853   Double_t eventW = 1;
 854   Double_t ptHard = 0;
 855   Double_t nTrials = 1; // Trials for MC trigger
 856   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 857
 858   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 859     // this is the part we only use when we have MC information
 860     AliMCEvent* mcEvent = MCEvent();
 861     //    AliStack *pStack = 0;
 862     if(!mcEvent){
 863       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 864       return;
 865     }
 866     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 867     if(pythiaGenHeader){
 868       nTrials = pythiaGenHeader->Trials();
 869       ptHard  = pythiaGenHeader->GetPtHard();
 870       int iProcessType = pythiaGenHeader->ProcessType();
 871       // 11 f+f -> f+f
 872       // 12 f+barf -> f+barf
 873       // 13 f+barf -> g+g
 874       // 28 f+g -> f+g
 875       // 53 g+g -> f+barf
 876       // 68 g+g -> g+g
 877       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 878       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 879     }
 880   }// (fAnalysisType&kMCESD)==kMCESD)
 881
 882
 883   // we simply fetch the tracks/mc particles as a list of AliVParticles
 884
 885   TList recParticles;
 886   TList genParticles;
 887
 888   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 889   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 890   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 891   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 892
 893   // Event control and counting ...
 894   // MC
 895   fh1PtHard->Fill(ptHard,eventW);
 896   fh1PtHardNoW->Fill(ptHard,1);
 897   fh1PtHardTrials->Fill(ptHard,nTrials);
 898
 899   // Real
 900   if(aod->GetPrimaryVertex()){// No vtx for pure MC
 901     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
 902   }
 903
 904   // the loops for rec and gen should be indentical... pass it to a separate
 905   // function ...
 906   // Jet Loop
 907   // Track Loop
 908   // Jet Jet Loop
 909   // Jet Track loop
 910
 911   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 912   FillTrackHistos(recParticles,kJetRec);
 913
 914   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 915   FillTrackHistos(genParticles,kJetGen);
 916
 917   // Here follows the jet matching part
 918   // delegate to a separated method?
 919
 920   FillMatchHistos(recJetsList,genJetsList);
 921
 922   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 923   PostData(1, fHistList);
 924 }
 925
 926 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 927
 928   if(iType>=kJetTypes){
 929     return;
 930   }
 931
 932   Int_t nJets = jetsList.GetEntries();
 933   fh1NJets[iType]->Fill(nJets);
 934
 935   if(nJets<=0)return;
 936
 937   Float_t ptOld = 1.E+32;
 938   Float_t pT0 = 0;
 939   Float_t pT1 = 0;
 940   Float_t phi0 = 0;
 941   Float_t phi1 = 0;
 942   Int_t ij0 = -1;
 943   Int_t ij1 = -1;
 944
 945
 946   for(int ij = 0;ij < nJets;ij++){
 947     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
 948     Float_t ptJet = jet->Pt();
 949     fh1PtJetsIn[iType]->Fill(ptJet);
 950     if(ptJet>ptOld){
 951       Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
 952     }
 953     ptOld = ptJet;
 954
 955     // find the dijets assume sorting and acceptance cut...
 956     if(ij==0){
 957       ij0 = ij;
 958       pT0 = ptJet;
 959       phi0 = jet->Phi();
 960       if(phi0<0)phi0+=TMath::Pi()*2.;
 961     }
 962     else if(ptJet>pT1){
 963       // take only the backward hemisphere??
 964       phi1 = jet->Phi();
 965       if(phi1<0)phi1+=TMath::Pi()*2.;
 966       Float_t dPhi = phi1 - phi0;
 967       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 968       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 969       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
 970         ij1 = ij;
 971         pT1 = ptJet;
 972       }
 973     }
 974     // fill jet histos for kmax jets
 975     if(ij<kMaxJets){
 976       Float_t phiJet = jet->Phi();
 977       if(phiJet<0)phiJet+=TMath::Pi()*2.;
 978       fh1TmpRho->Reset();
 979       fh1PtIn[iType][ij]->Fill(ptJet);
 980       // fill leading jets...
 981       if(ij==0){
 982         fh2LeadingJetPtJetPhi[iType]->Fill(ptJet,ptJet);
 983         if(ij==0&&iType==0&&fDebug>1){
 984           Printf("%d %3.3f %p %s",iType,ptJet,jet,fBranchRec.Data());
 985         }
 986       }
 987       if(particlesList.GetSize()){
 988
 989         // Particles... correlated with jets...
 990         for(int it = 0;it<particlesList.GetEntries();++it){
 991           AliVParticle *part = (AliVParticle*)particlesList.At(it);
 992           Float_t deltaR = jet->DeltaR(part);
 993           if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
 994         }
 995         // fill the jet shapes
 996         Float_t rhoSum = 0;
 997         for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
 998           Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
 999           Float_t rho = fh1TmpRho->GetBinContent(ibx);
1000           rhoSum += rho;
1001           fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
1002           fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
1003         }
1004       }// if we have particles
1005     }// ij < kMaxJets
1006   }// Jet Loop
1007
1008
1009   // do something with dijets...
1010   if(ij0>=0&&ij1>0){
1011     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1012     Double_t ptJet0 = jet0->Pt();
1013     Double_t phiJet0 = jet0->Phi();
1014     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1015
1016     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1017     Double_t ptJet1 = jet1->Pt();
1018     Double_t phiJet1 = jet1->Phi();
1019     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1020
1021     Float_t deltaPhi = phiJet0 - phiJet1;
1022     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1023     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1024     deltaPhi = TMath::Abs(deltaPhi);
1025     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1026
1027     Float_t asym = 9999;
1028     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1029       fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1030       fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1031       fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1032       Float_t minv = 2.*(jet0->P()*jet1->P()-
1033                          jet0->Px()*jet1->Px()-
1034                          jet0->Py()*jet1->Py()-
1035                          jet0->Pz()*jet1->Pz());    // assume mass == 0;
1036       if(minv<0)minv=0; // prevent numerical instabilities
1037       minv = TMath::Sqrt(minv);
1038       fh1DijetMinv[iType]->Fill(minv);
1039   }
1040
1041
1042
1043   // count down the jets above thrueshold
1044   Int_t nOver = nJets;
1045   if(nOver>0){
1046     TIterator *jetIter = jetsList.MakeIterator();
1047     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1048     Float_t pt = tmpJet->Pt();
1049     if(tmpJet){
1050       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1051         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1052         while(pt<ptCut&&tmpJet){
1053           nOver--;
1054           tmpJet = (AliAODJet*)(jetIter->Next());
1055           if(tmpJet){
1056             pt = tmpJet->Pt();
1057           }
1058         }
1059         if(nOver<=0)break;
1060         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1061       }
1062     }
1063     delete jetIter;
1064   }
1065 }
1066
1067 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1068
1069   Int_t nTrackOver = particlesList.GetSize();
1070   // do the same for tracks and jets
1071   if(nTrackOver>0){
1072     TIterator *trackIter = particlesList.MakeIterator();
1073     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1074     Float_t pt = tmpTrack->Pt();
1075     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1076       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1077       while(pt<ptCut&&tmpTrack){
1078         nTrackOver--;
1079         tmpTrack = (AliVParticle*)(trackIter->Next());
1080         if(tmpTrack){
1081           pt = tmpTrack->Pt();
1082         }
1083       }
1084       if(nTrackOver<=0)break;
1085       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1086     }
1087
1088
1089     trackIter->Reset();
1090     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1091     Float_t sumPt = 0;
1092
1093     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1094       Float_t tmpPt = tmpTrack->Pt();
1095       fh1PtTracksIn[iType]->Fill(tmpPt);
1096       sumPt += tmpPt;
1097       Float_t tmpPhi = tmpTrack->Phi();
1098       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1099       if(tmpTrack==leading){
1100         fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
1101         fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
1102         continue;
1103       }
1104     }
1105     fh1SumPtTrack[iType]->Fill(sumPt);
1106     delete trackIter;
1107   }
1108
1109 }
1110
1111
1112 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1113
1114
1115   // Fill al the matching histos
1116   // It is important that the acceptances for the mathing are as large as possible
1117   // to avoid false matches on the edge of acceptance
1118   // therefore we add some extra matching jets as overhead
1119
1120   static TArrayI aGenIndex(recJetsList.GetEntries());
1121   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1122
1123   static TArrayI aRecIndex(genJetsList.GetEntries());
1124   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1125
1126   if(fDebug){
1127     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1128     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1129   }
1130   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1131                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1132                                             aGenIndex,aRecIndex,fDebug);
1133
1134   if(fDebug){
1135     for(int i = 0;i< aGenIndex.GetSize();++i){
1136       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1137     }
1138     for(int i = 0;i< aRecIndex.GetSize();++i){
1139       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1140     }
1141   }
1142
1143   Double_t container[6];
1144   Double_t containerPhiZ[6];
1145
1146   // loop over generated jets
1147   // consider the
1148   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1149     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1150     Double_t ptGen = genJet->Pt();
1151     Double_t phiGen = genJet->Phi();
1152     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1153     Double_t etaGen = genJet->Eta();
1154     container[3] = ptGen;
1155     container[4] = etaGen;
1156     container[5] = phiGen;
1157     fhnJetContainer[kStep0]->Fill(&container[3]);
1158     if(JetSelected(genJet)){
1159       fhnJetContainer[kStep1]->Fill(&container[3]);
1160       Int_t ir = aRecIndex[ig];
1161       if(ir>=0&&ir<recJetsList.GetEntries()){
1162         fhnJetContainer[kStep2]->Fill(&container[3]);
1163         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1164         if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
1165         if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
1166       }
1167     }
1168   }// loop over generated jets used for matching...
1169
1170
1171
1172   // loop over reconstructed jets
1173   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1174     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1175     Double_t etaRec = recJet->Eta();
1176     Double_t ptRec = recJet->Pt();
1177     Double_t phiRec = recJet->Phi();
1178     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1179     // do something with dijets...
1180
1181     container[0] = ptRec;
1182     container[1] = etaRec;
1183     container[2] = phiRec;
1184     containerPhiZ[0] = ptRec;
1185     containerPhiZ[1] = phiRec;
1186
1187     fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1188     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1189
1190     if(JetSelected(recJet)){
1191       fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1192       // Fill Correlation
1193       Int_t ig = aGenIndex[ir];
1194       if(ig>=0 && ig<genJetsList.GetEntries()){
1195         fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1196         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1197         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1198         Double_t ptGen  = genJet->Pt();
1199         Double_t phiGen = genJet->Phi();
1200         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1201         Double_t etaGen = genJet->Eta();
1202
1203         container[3] = ptGen;
1204         container[4] = etaGen;
1205         container[5] = phiGen;
1206         containerPhiZ[3] = ptGen;
1207         //
1208         // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1209         //
1210         if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1211         fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1212         fhnCorrelation->Fill(container);
1213         if(ptGen>0){
1214           Float_t delta = (ptRec-ptGen)/ptGen;
1215           fh2RelPtFGen->Fill(ptGen,delta);
1216           fh2PtFGen->Fill(ptGen,ptRec);
1217         }
1218         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1219       }
1220       else{
1221         containerPhiZ[3] = 0;
1222         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1223       }
1224     }// loop over reconstructed jets
1225   }
1226   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1227 }
1228
1229
1230 void AliAnalysisTaskJetSpectrum2::UserExecOld(Option_t */*option*/)
1231 {
1232
1233
1234   Bool_t selected = kTRUE;
1235
1236
1237
1238
1239   if(fUseGlobalSelection&&fEventSelectionMask==0){
1240     selected = AliAnalysisHelperJetTasks::Selected();
1241   }
1242   else if(fUseGlobalSelection&&fEventSelectionMask>0){
1243     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
1244   }
1245
1246   if(fEventClass>0){
1247     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
1248   }
1249
1250   if(!selected){
1251     // no selection by the service task, we continue
1252     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
1253     PostData(1, fHistList);
1254     return;
1255   }
1256
1257
1258   //
1259   // Execute analysis for current event
1260   //
1261   static AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1262
1263
1264   // take all other information from the aod we take the tracks from
1265   static AliAODEvent *aod = 0;
1266   if(!aod){
1267     if(fUseAODTrackInput)aod = fAODIn;
1268     aod = fAODOut;
1269   }
1270   static AliAODEvent *fAOD = 0;
1271   if(!fAOD){
1272     if(fUseAODJetInput)fAOD = fAODIn;
1273     fAOD = fAODOut;
1274   }
1275
1276   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
1277
1278
1279   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1280
1281   if(!aodH){
1282     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1283     return;
1284   }
1285
1286   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1287   TClonesArray *aodRecJets = 0;
1288
1289
1290
1291   if(fAOD&&!aodRecJets){
1292     aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
1293   }
1294   if(fAODExtension&&!aodRecJets){
1295     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
1296   }
1297
1298   if(!aodRecJets){
1299     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
1300     return;
1301   }
1302
1303   Int_t nJets = aodRecJets->GetEntriesFast();
1304
1305   TClonesArray *aodGenJets = 0;
1306   if(fBranchGen.Length()>0){
1307     if(fAOD&&!aodGenJets){
1308       aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
1309     }
1310     if(fAODExtension&&!aodGenJets){
1311       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
1312     }
1313
1314     if(!aodGenJets){
1315       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
1316       return;
1317     }
1318   }
1319
1320
1321   // ==== General variables needed
1322   // We use statice array, not to fragment the memory
1323   AliAODJet genJets[kMaxJets];
1324   Int_t nGenJets = 0;
1325   AliAODJet recJets[kMaxJets];
1326   Int_t nRecJets = 0;
1327
1328   // These will replace the arrays...
1329   TList genJetsList;         // full acceptance
1330   TList genJetsListCut;   // acceptance cut
1331   TList recJetsList;         // full acceptance
1332   TList recJetsListCut;   // acceptance cut
1333
1334
1335   GetListOfJets(&genJetsList,aodGenJets,0);
1336   GetListOfJets(&genJetsListCut,aodGenJets,1);
1337
1338   if(fBranchGen.Length()>0){
1339     GetListOfJets(&recJetsList,aodRecJets,0);
1340     GetListOfJets(&recJetsListCut,aodRecJets,1);
1341   }
1342   ///////////////////////////
1343
1344
1345
1346
1347   if(fBkgSubtraction){
1348      AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
1349      if(!evBkg){
1350        Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
1351        return;
1352      }
1353
1354
1355      ///just to start: some very simple plots containing rho, sigma and area of the background.
1356
1357
1358      Float_t pthardest=0.;
1359
1360      if(nJets!=0){
1361
1362        Float_t bkg1=evBkg->GetBackground(0);
1363        Float_t bkg2=evBkg->GetBackground(1);
1364        Float_t bkg3=evBkg->GetBackground(2);
1365        Float_t sigma1=evBkg->GetSigma(0);
1366        Float_t sigma2=evBkg->GetSigma(1);
1367        Float_t sigma3=evBkg->GetSigma(2);
1368        Float_t area1=evBkg->GetMeanarea(0);
1369        Float_t area2=evBkg->GetMeanarea(1);
1370        Float_t area3=evBkg->GetMeanarea(2);
1371        fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
1372        fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
1373        fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
1374        fh1Sigma1->Fill(sigma1);
1375        fh1Sigma2->Fill(sigma2);
1376        fh1Sigma3->Fill(sigma3);
1377        fh1Area1->Fill(area1);
1378        fh1Area2->Fill(area2);
1379        fh1Area3->Fill(area3);
1380
1381        Int_t iSubJetCounter = 0;
1382        for(Int_t k=0;k<nJets;k++){
1383          AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
1384          fh1Ptjet->Fill(jet->Pt());
1385          Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
1386          Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
1387          Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
1388          if(ptsub2<0.) ptsub2=0.;
1389          if(ptsub1<0.) ptsub1=0.;
1390          if(ptsub3<0.) ptsub3=0.;
1391          Float_t err1=sigma1*sqrt(area1);
1392          Float_t err2=sigma2*sqrt(area2);
1393          Float_t err3=sigma3*sqrt(area3);
1394          fh1Ptjetsub1->Fill(ptsub1);
1395          fh1Ptjetsub2->Fill(ptsub2);
1396          fh1Ptjetsub3->Fill(ptsub3);
1397          if(k==0) {
1398            pthardest=jet->Pt();
1399            fh1Ptjethardest->Fill(pthardest);
1400            fh1Ptjetsubhardest1->Fill(ptsub1);
1401            fh1Ptjetsubhardest2->Fill(ptsub2);
1402            fh1Ptjetsubhardest3->Fill(ptsub3);
1403            if (ptsub1 > 0)
1404              fh2Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
1405            if (ptsub2 > 0)
1406              fh2Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
1407            if (ptsub3 > 0)
1408              fh2Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
1409          }
1410
1411          Float_t ptsub=0.;
1412          if(fFillCorrBkg==1) ptsub=ptsub1;
1413          if(fFillCorrBkg==2) ptsub=ptsub2;
1414          if(fFillCorrBkg==3) ptsub=ptsub3;
1415          if(ptsub>0){// avoid unphysical jets pT
1416            Float_t subphi=jet->Phi();
1417            Float_t subtheta=jet->Theta();
1418            Float_t subpz = ptsub/TMath::Tan(subtheta);
1419            Float_t subpx=ptsub*TMath::Cos(subphi);
1420            Float_t subpy=ptsub * TMath::Sin(subphi);
1421            Float_t subp  = TMath::Sqrt(ptsub*ptsub+subpz*subpz);
1422            if(k<kMaxJets){// only store the jets which are assoicated to anohter one
1423              genJets[iSubJetCounter].SetPxPyPzE(subpx,subpy,subpz,subp);
1424              iSubJetCounter++;
1425            }
1426          }
1427        }
1428        nGenJets = iSubJetCounter;
1429        fh2Rhovspthardest1->Fill(pthardest,bkg1);
1430        fh2Rhovspthardest2->Fill(pthardest,bkg2);
1431        fh2Rhovspthardest3->Fill(pthardest,bkg3);
1432      }
1433   }// background subtraction
1434
1435
1436   Double_t eventW = 1;
1437   Double_t ptHard = 0;
1438   Double_t nTrials = 1; // Trials for MC trigger
1439
1440   if(fUseExternalWeightOnly){
1441     eventW = fExternalWeight;
1442   }
1443
1444   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
1445   //  if(fDebug>0)aodH->SetFillAOD(kFALSE);
1446   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1447   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
1448     // this is the part we only use when we have MC information
1449     AliMCEvent* mcEvent = MCEvent();
1450     //    AliStack *pStack = 0;
1451     if(!mcEvent){
1452       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
1453       return;
1454     }
1455     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1456     Int_t iCount = 0;
1457     if(pythiaGenHeader){
1458       nTrials = pythiaGenHeader->Trials();
1459       ptHard  = pythiaGenHeader->GetPtHard();
1460       int iProcessType = pythiaGenHeader->ProcessType();
1461       // 11 f+f -> f+f
1462       // 12 f+barf -> f+barf
1463       // 13 f+barf -> g+g
1464       // 28 f+g -> f+g
1465       // 53 g+g -> f+barf
1466       // 68 g+g -> g+g
1467       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
1468       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
1469
1470       // fetch the pythia generated jets only to be used here
1471       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
1472       AliAODJet pythiaGenJets[kMaxJets];
1473       for(int ip = 0;ip < nPythiaGenJets;++ip){
1474         if(iCount>=kMaxJets)continue;
1475         Float_t p[4];
1476         pythiaGenHeader->TriggerJet(ip,p);
1477         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
1478
1479         if(fBranchGen.Length()==0){
1480           /*
1481           if(fLimitGenJetEta){
1482             if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
1483                pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
1484           }
1485           */
1486           if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
1487           // if we have MC particles and we do not read from the aod branch
1488           // use the pythia jets
1489           if(!fBkgSubtraction){
1490             genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
1491             iCount++;
1492           }
1493         }
1494       }
1495     }
1496     if(fBranchGen.Length()==0&&!fBkgSubtraction)nGenJets = iCount;
1497   }// (fAnalysisType&kMCESD)==kMCESD)
1498
1499
1500   // we simply fetch the tracks/mc particles as a list of AliVParticles
1501
1502   TList recParticles;
1503   TList genParticles;
1504
1505
1506
1507
1508   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
1509   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1510   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1511   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
1512
1513
1514   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1515   fh1PtHard->Fill(ptHard,eventW);
1516   fh1PtHardNoW->Fill(ptHard,1);
1517   fh1PtHardTrials->Fill(ptHard,nTrials);
1518   fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
1519
1520
1521   // If we set a second branch for the input jets fetch this
1522   if(fBranchGen.Length()>0&&!fBkgSubtraction){
1523     if(aodGenJets){
1524       Int_t iCount = 0;
1525       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
1526         if(iCount>=kMaxJets)continue;
1527         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
1528         if(!tmp)continue;
1529         /*
1530         if(fLimitGenJetEta){
1531           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
1532              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
1533         }
1534         */
1535         if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
1536         genJets[iCount] = *tmp;
1537         iCount++;
1538       }
1539       nGenJets = iCount;
1540     }
1541     else{
1542       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
1543       if(fDebug>2)fAOD->Print();
1544     }
1545   }
1546
1547   fh1NGenJets->Fill(nGenJets);
1548   // We do not want to exceed the maximum jet number
1549   nGenJets = TMath::Min(nGenJets,kMaxJets);
1550
1551   // Fetch the reconstructed jets...
1552
1553   nRecJets = aodRecJets->GetEntries();
1554
1555   nRecJets = aodRecJets->GetEntries();
1556   fh1NRecJets->Fill(nRecJets);
1557
1558   // Do something with the all rec jets
1559   Int_t nRecOver = nRecJets;
1560
1561   // check that the jets are sorted
1562   Float_t ptOld = 999999;
1563   for(int ir = 0;ir < nRecJets;ir++){
1564     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
1565     Float_t tmpPt = tmp->Pt();
1566     if(tmpPt>ptOld){
1567       Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
1568     }
1569     ptOld = tmpPt;
1570   }
1571
1572
1573   if(nRecOver>0){
1574     TIterator *recIter = aodRecJets->MakeIterator();
1575     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
1576     Float_t pt = tmpRec->Pt();
1577     if(tmpRec){
1578       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1579         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1580         while(pt<ptCut&&tmpRec){
1581           nRecOver--;
1582           tmpRec = (AliAODJet*)(recIter->Next());
1583           if(tmpRec){
1584             pt = tmpRec->Pt();
1585           }
1586         }
1587         if(nRecOver<=0)break;
1588         fh2NRecJetsPt->Fill(ptCut,nRecOver);
1589       }
1590     }
1591     recIter->Reset();
1592     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
1593     Float_t phi = leading->Phi();
1594     if(phi<0)phi+=TMath::Pi()*2.;
1595     pt = leading->Pt();
1596     while((tmpRec = (AliAODJet*)(recIter->Next()))){
1597       Float_t tmpPt = tmpRec->Pt();
1598       fh1PtJetsRecIn->Fill(tmpPt);
1599       if(tmpRec==leading){
1600         fh1PtJetsLeadingRecIn->Fill(tmpPt);
1601         continue;
1602       }
1603       // correlation
1604       Float_t tmpPhi =  tmpRec->Phi();
1605
1606       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1607       Float_t dPhi = phi - tmpRec->Phi();
1608       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1609       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1610       //  Float_t dEta = eta - tmpRec->Eta();
1611       //      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1612       //      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1613     }
1614     delete recIter;
1615   }
1616
1617   Int_t nTrackOver = recParticles.GetSize();
1618   // do the same for tracks and jets
1619   if(nTrackOver>0){
1620     TIterator *recIter = recParticles.MakeIterator();
1621     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1622     Float_t pt = tmpRec->Pt();
1623     //    Printf("Leading track p_t %3.3E",pt);
1624     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1625       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1626       while(pt<ptCut&&tmpRec){
1627         nTrackOver--;
1628         tmpRec = (AliVParticle*)(recIter->Next());
1629         if(tmpRec){
1630           pt = tmpRec->Pt();
1631         }
1632       }
1633       if(nTrackOver<=0)break;
1634       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1635     }
1636
1637     recIter->Reset();
1638     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1639     Float_t phi = leading->Phi();
1640     if(phi<0)phi+=TMath::Pi()*2.;
1641     pt  = leading->Pt();
1642     while((tmpRec = (AliVParticle*)(recIter->Next()))){
1643       Float_t tmpPt = tmpRec->Pt();
1644       fh1PtTracksRecIn->Fill(tmpPt);
1645       if(tmpRec==leading){
1646         fh1PtTracksLeadingRecIn->Fill(tmpPt);
1647         continue;
1648       }
1649       // correlation
1650       Float_t tmpPhi =  tmpRec->Phi();
1651
1652       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1653       Float_t dPhi = phi - tmpRec->Phi();
1654       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1655       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1656       //      Float_t dEta = eta - tmpRec->Eta();
1657       //      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1658       //      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1659     }
1660     delete recIter;
1661   }
1662
1663   if(genParticles.GetSize()){
1664     TIterator *genIter = genParticles.MakeIterator();
1665     AliVParticle *tmpGen = 0;
1666     while((tmpGen = (AliVParticle*)(genIter->Next()))){
1667       if(TMath::Abs(tmpGen->Eta())<0.9){
1668         Float_t tmpPt = tmpGen->Pt();
1669         fh1PtTracksGenIn->Fill(tmpPt);
1670       }
1671     }
1672     delete genIter;
1673   }
1674
1675   nRecJets = TMath::Min(nRecJets,kMaxJets);
1676   nJets = TMath::Min(nJets,kMaxJets);
1677   Int_t iCountRec = 0;
1678   for(int ir = 0;ir < nRecJets;++ir){
1679     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1680     if(!tmp)continue;
1681     if(tmp->Pt()<fMinJetPt)continue;
1682     recJets[iCountRec] = *tmp;
1683     iCountRec++;
1684   }
1685   nRecJets = iCountRec;
1686
1687
1688   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1689   // Relate the jets
1690   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
1691   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
1692
1693
1694   for(int i = 0;i<kMaxJets;++i){
1695     iGenIndex[i] = iRecIndex[i] = -1;
1696   }
1697
1698
1699   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1700                                             iGenIndex,iRecIndex,fDebug);
1701
1702   static TArrayI aGenIndex(recJetsListCut.GetEntries());
1703   if(aGenIndex.GetSize()<recJetsListCut.GetEntries())aGenIndex.Set(recJetsListCut.GetEntries());
1704
1705   static TArrayI aRecIndex(genJetsList.GetEntries());
1706   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1707
1708   if(fDebug){
1709     Printf("New Gens List %d rec index Array %d count %d",genJetsList.GetEntries(),aRecIndex.GetSize(),nGenJets);
1710     Printf("New Rec List %d gen indey Array %d count %d",recJetsListCut.GetEntries(),aGenIndex.GetSize(),nRecJets);
1711   }
1712   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,nGenJets,&recJetsListCut,nRecJets,
1713                                             aGenIndex,aRecIndex,fDebug);
1714
1715
1716   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1717
1718   if(fDebug){
1719     for(int i = 0;i<kMaxJets;++i){
1720       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1721       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1722     }
1723   }
1724
1725
1726
1727
1728   Double_t container[6];
1729   Double_t containerPhiZ[6];
1730
1731   // loop over generated jets
1732
1733   // radius; tmp, get from the jet header itself
1734   // at some point, how todeal woht FastJet on MC?
1735   Float_t radiusGen =0.4;
1736   Float_t radiusRec =0.4;
1737
1738   for(int ig = 0;ig < nGenJets;++ig){
1739     Double_t ptGen = genJets[ig].Pt();
1740     Double_t phiGen = genJets[ig].Phi();
1741     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1742     Double_t etaGen = genJets[ig].Eta();
1743     fh1PtJetsGenIn->Fill(ptGen);
1744
1745     container[3] = ptGen;
1746     container[4] = etaGen;
1747     container[5] = phiGen;
1748     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1749     Int_t ir = iRecIndex[ig];
1750
1751     if(TMath::Abs(etaGen)<fJetRecEtaWindow){
1752       fh1TmpRho->Reset();
1753
1754       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1755       fh1PtGenIn[ig]->Fill(ptGen,eventW);
1756       // fill the fragmentation function
1757       for(int it = 0;it<genParticles.GetEntries();++it){
1758         AliVParticle *part = (AliVParticle*)genParticles.At(it);
1759         Float_t deltaR = genJets[ig].DeltaR(part);
1760         if(ptGen>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1761         if(deltaR<radiusGen&&ptGen>0){
1762           Float_t z = part->Pt()/ptGen;
1763           Float_t lnz =  -1.*TMath::Log(z);
1764           fh2FragGen[ig]->Fill(z,ptGen,eventW);
1765           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1766         }
1767
1768       }
1769       Float_t rhoSum = 0;
1770       for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
1771         Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1772         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1773         rhoSum += rho;
1774         fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1775         fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
1776       }
1777     }
1778     if(ir>=0&&ir<nRecJets){
1779       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1780       Double_t etaRec = recJets[ir].Eta();
1781       if(TMath::Abs(etaRec)<fJetRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1782       if(TMath::Abs(etaRec)<fJetRecEtaWindow&&TMath::Abs(etaGen)<fJetRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
1783     }
1784   }// loop over generated jets
1785
1786   Float_t sumPt = 0;
1787   for(int it = 0;it<recParticles.GetEntries();++it){
1788     AliVParticle *part = (AliVParticle*)recParticles.At(it);
1789     // fill sum pt and P_t of all paritles
1790     if(TMath::Abs(part->Eta())<0.9){
1791       Float_t pt = part->Pt();
1792       fh1PtTrackRec->Fill(pt,eventW);
1793       //      fh2TrackPtTrackPhi->Fill(pt,part->Phi());
1794       sumPt += pt;
1795     }
1796   }
1797   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1798   fh1SumPtTrackRec->Fill(sumPt,eventW);
1799
1800
1801   // loop over reconstructed jets
1802   for(int ir = 0;ir < nRecJets;++ir){
1803     Double_t etaRec = recJets[ir].Eta();
1804     Double_t ptRec = recJets[ir].Pt();
1805     Double_t phiRec = recJets[ir].Phi();
1806     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1807
1808
1809   // do something with dijets...
1810
1811     container[0] = ptRec;
1812     container[1] = etaRec;
1813     container[2] = phiRec;
1814     containerPhiZ[0] = ptRec;
1815     containerPhiZ[1] = phiRec;
1816     if(ptRec>30.&&fDebug>2){
1817       // need to cast to int, otherwise the printf overwrites
1818       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1819       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1820       if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
1821       //  aodH->SetFillAOD(kTRUE);
1822       fAOD->GetHeader()->Print();
1823       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
1824       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1825         AliAODTrack *tr = fAOD->GetTrack(it);
1826         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1827         tr->Print();
1828         //      tr->Dump();
1829         if(fESD){
1830           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1831           esdTr->Print("");
1832           // esdTr->Dump();
1833         }
1834       }
1835     }
1836
1837
1838     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1839     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1840
1841     Float_t zLeading = -1;
1842     if(TMath::Abs(etaRec)<fJetRecEtaWindow){
1843       //      fh2JetPtJetPhi->Fill(ptRec,phiRec);
1844       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1845       fh1PtRecIn[ir]->Fill(ptRec,eventW);
1846       // fill the fragmentation function
1847
1848       fh1TmpRho->Reset();
1849
1850       for(int it = 0;it<recParticles.GetEntries();++it){
1851         AliVParticle *part = (AliVParticle*)recParticles.At(it);
1852         Float_t eta = part->Eta();
1853         if(TMath::Abs(eta)<0.9){
1854           Float_t phi = part->Phi();
1855           if(phi<0)phi+=TMath::Pi()*2.;
1856           Float_t dPhi = phi - phiRec;
1857           Float_t dEta = eta - etaRec;
1858           if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1859           if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1860           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1861           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1862         }
1863
1864         Float_t deltaR = recJets[ir].DeltaR(part);
1865         if(ptRec>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1866
1867
1868         if(deltaR<radiusRec&&ptRec>0){
1869           Float_t z = part->Pt()/ptRec;
1870           if(z>zLeading)zLeading=z;
1871           Float_t lnz =  -1.*TMath::Log(z);
1872           fh2FragRec[ir]->Fill(z,ptRec,eventW);
1873           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1874         }
1875       }
1876       // fill the jet shapes
1877       Float_t rhoSum = 0;
1878       for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1879         Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1880         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1881         rhoSum += rho;
1882         fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1883         fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1884       }
1885     }
1886
1887     // Fill Correlation
1888     Int_t ig = iGenIndex[ir];
1889     if(ig>=0 && ig<nGenJets){
1890       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1891       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1892       Double_t ptGen  = genJets[ig].Pt();
1893       Double_t phiGen = genJets[ig].Phi();
1894       if(phiGen<0)phiGen+=TMath::Pi()*2.;
1895       Double_t etaGen = genJets[ig].Eta();
1896
1897       container[3] = ptGen;
1898       container[4] = etaGen;
1899       container[5] = phiGen;
1900       containerPhiZ[3] = ptGen;
1901       //
1902       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1903       //
1904
1905       if(TMath::Abs(etaGen)<fJetRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1906       if(TMath::Abs(etaRec)<fJetRecEtaWindow){
1907         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1908         fhnCorrelation->Fill(container);
1909         if(ptGen>0){
1910           Float_t delta = (ptRec-ptGen)/ptGen;
1911           fh2RelPtFGen->Fill(ptGen,delta,eventW);
1912         }
1913         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1914       }// if etarec in window
1915     }
1916     else{
1917       containerPhiZ[3] = 0;
1918       if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1919     }
1920   }// loop over reconstructed jets
1921
1922
1923   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1924   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1925   PostData(1, fHistList);
1926 }
1927
1928
1929 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1930   //
1931   // Create the particle container for the correction framework manager and
1932   // link it
1933   //
1934   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1935   const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1936   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1937   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1938   const Double_t kZmin = 0., kZmax = 1;
1939
1940   // can we neglect migration in eta and phi?
1941   // phi should be no problem since we cover full phi and are phi symmetric
1942   // eta migration is more difficult  due to needed acceptance correction
1943   // in limited eta range
1944
1945   //arrays for the number of bins in each dimension
1946   Int_t iBin[kNvar];
1947   iBin[0] = 320; //bins in pt
1948   iBin[1] =  1; //bins in eta
1949   iBin[2] = 1; // bins in phi
1950
1951   //arrays for lower bounds :
1952   Double_t* binEdges[kNvar];
1953   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1954     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1955
1956   //values for bin lower bounds
1957   //  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);
1958   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1959   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1960   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1961
1962
1963   for(int i = 0;i < kMaxStep*2;++i){
1964     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1965     for (int k=0; k<kNvar; k++) {
1966       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1967     }
1968   }
1969   //create correlation matrix for unfolding
1970   Int_t thnDim[2*kNvar];
1971   for (int k=0; k<kNvar; k++) {
1972     //first half  : reconstructed
1973     //second half : MC
1974     thnDim[k]      = iBin[k];
1975     thnDim[k+kNvar] = iBin[k];
1976   }
1977
1978   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1979   for (int k=0; k<kNvar; k++) {
1980     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1981     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1982   }
1983   fhnCorrelation->Sumw2();
1984
1985   // for second correlation histogram
1986
1987
1988   const Int_t kNvarPhiZ   = 4;
1989   //arrays for the number of bins in each dimension
1990   Int_t iBinPhiZ[kNvarPhiZ];
1991   iBinPhiZ[0] = 80; //bins in pt
1992   iBinPhiZ[1] = 72; //bins in phi
1993   iBinPhiZ[2] = 20; // bins in Z
1994   iBinPhiZ[3] = 80; //bins in ptgen
1995
1996   //arrays for lower bounds :
1997   Double_t* binEdgesPhiZ[kNvarPhiZ];
1998   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1999     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
2000
2001   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
2002   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
2003   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
2004   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
2005
2006   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
2007   for (int k=0; k<kNvarPhiZ; k++) {
2008     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
2009   }
2010   fhnCorrelationPhiZRec->Sumw2();
2011
2012
2013   // Add a histogram for Fake jets
2014
2015   for(Int_t ivar = 0; ivar < kNvar; ivar++)
2016     delete [] binEdges[ivar];
2017
2018   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
2019     delete [] binEdgesPhiZ[ivar];
2020
2021 }
2022
2023 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
2024 {
2025 // Terminate analysis
2026 //
2027     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
2028 }
2029
2030
2031 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
2032
2033   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
2034
2035   Int_t iCount = 0;
2036   if(type==kTrackAOD){
2037     AliAODEvent *aod = 0;
2038     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2039     else aod = AODEvent();
2040     if(!aod){
2041       return iCount;
2042     }
2043     for(int it = 0;it < aod->GetNumberOfTracks();++it){
2044       AliAODTrack *tr = aod->GetTrack(it);
2045       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
2046       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
2047       if(tr->Pt()<fMinTrackPt)continue;
2048       if(fDebug>0){
2049         if(tr->Pt()>20){
2050           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
2051           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
2052           tr->Print();
2053           //    tr->Dump();
2054           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
2055           if(fESD){
2056             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
2057             esdTr->Print("");
2058             // esdTr->Dump();
2059           }
2060         }
2061       }
2062       list->Add(tr);
2063       iCount++;
2064     }
2065   }
2066   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
2067     AliMCEvent* mcEvent = MCEvent();
2068     if(!mcEvent)return iCount;
2069     // we want to have alivpartilces so use get track
2070     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
2071       if(!mcEvent->IsPhysicalPrimary(it))continue;
2072       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
2073       if(part->Pt()<fMinTrackPt)continue;
2074       if(type == kTrackKineAll){
2075         list->Add(part);
2076         iCount++;
2077       }
2078       else if(type == kTrackKineCharged){
2079         if(part->Particle()->GetPDG()->Charge()==0)continue;
2080         list->Add(part);
2081         iCount++;
2082       }
2083     }
2084   }
2085   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
2086     AliAODEvent *aod = 0;
2087     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2088     else aod = AODEvent();
2089     if(!aod)return iCount;
2090     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2091     if(!tca)return iCount;
2092     for(int it = 0;it < tca->GetEntriesFast();++it){
2093       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
2094       if(part->Pt()<fMinTrackPt)continue;
2095       if(!part->IsPhysicalPrimary())continue;
2096       if(type == kTrackAODMCAll){
2097         list->Add(part);
2098         iCount++;
2099       }
2100       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
2101         if(part->Charge()==0)continue;
2102         if(kTrackAODMCCharged){
2103           list->Add(part);
2104         }
2105         else {
2106           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
2107           list->Add(part);
2108         }
2109         iCount++;
2110       }
2111     }
2112   }// AODMCparticle
2113   list->Sort();
2114   return iCount;
2115
2116 }
2117
2118
2119
2120 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
2121   Bool_t selected = false;
2122
2123   if(!jet)return selected;
2124
2125   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
2126     selected = kTRUE;
2127   }
2128   return selected;
2129
2130 }
2131
2132 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
2133
2134   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
2135   Int_t iCount = 0;
2136
2137   if(!jarray){
2138     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
2139     return iCount;
2140   }
2141
2142
2143   for(int ij=0;ij<jarray->GetEntries();ij++){
2144     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
2145     if(!jet)continue;
2146     if(type==0){
2147       // No cut at all, main purpose here is sorting
2148       list->Add(jet);
2149       iCount++;
2150     }
2151     else if(type == 1){
2152       // eta cut
2153       if(JetSelected(jet)){
2154         list->Add(jet);
2155         iCount++;
2156       }
2157     }
2158   }
2159
2160   list->Sort();
2161   return iCount;
2162
2163 }
2164
2165