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   fNMatchJets(5),
  88   fFillCorrBkg(0),
  89   fFilterMask(0),
  90   fEventSelectionMask(0),
  91   fAnalysisType(0),
  92   fTrackTypeRec(kTrackUndef),
  93   fTrackTypeGen(kTrackUndef),
  94   fEventClass(0),
  95   fAvgTrials(1),
  96   fExternalWeight(1),
  97   fJetRecEtaWindow(0.5),
  98   fTrackRecEtaWindow(0.5),
  99   fMinJetPt(0),
 100   fMinTrackPt(0.15),
 101   fDeltaPhiWindow(90./180.*TMath::Pi()),
 102   fh1Xsec(0x0),
 103   fh1Trials(0x0),
 104   fh1PtHard(0x0),
 105   fh1PtHardNoW(0x0),
 106   fh1PtHardTrials(0x0),
 107   fh1ZVtx(0x0),
 108   fh1NGenJets(0x0),
 109   fh1NRecJets(0x0),
 110   fh1PtTrackRec(0x0),
 111   fh1SumPtTrackRec(0x0),
 112   fh1SumPtTrackAreaRec(0x0),
 113   fh1TmpRho(0x0),
 114   fh1PtJetsRecIn(0x0),
 115   fh1PtJetsGenIn(0x0),
 116   fh1PtJetsLeadingRecIn(0x0),
 117   fh1PtTracksRecIn(0x0),
 118   fh1PtTracksLeadingRecIn(0x0),
 119   fh1PtTracksGenIn(0x0),
 120   fh2NRecJetsPt(0x0),
 121   fh2NRecTracksPt(0x0),
 122   fh2NGenJetsPt(0x0),
 123   fh2NGenTracksPt(0x0),
 124   fh2PtFGen(0x0),
 125   fh2RelPtFGen(0x0),
 126   fh1Bkg1(0x0),
 127   fh1Bkg2(0x0),
 128   fh1Bkg3(0x0),
 129   fh1Sigma1(0x0),
 130   fh1Sigma2(0x0),
 131   fh1Sigma3(0x0),
 132   fh1Area1(0x0),
 133   fh1Area2(0x0),
 134   fh1Area3(0x0),
 135   fh1Ptjet(0x0),
 136   fh1Ptjetsub1(0x0),
 137   fh1Ptjetsub2(0x0),
 138   fh1Ptjetsub3(0x0),
 139   fh1Ptjethardest(0x0),
 140   fh1Ptjetsubhardest1(0x0),
 141   fh1Ptjetsubhardest2(0x0),
 142   fh1Ptjetsubhardest3(0x0),
 143   fh2Rhovspthardest1(0x0),
 144   fh2Rhovspthardest2(0x0),
 145   fh2Rhovspthardest3(0x0),
 146   fh2Errorvspthardest1(0x0),
 147   fh2Errorvspthardest2(0x0),
 148   fh2Errorvspthardest3(0x0),
 149   fHistList(0x0)
 150 {
 151   for(int i = 0;i < kMaxStep*2;++i){
 152     fhnJetContainer[i] = 0;
 153   }
 154   for(int i = 0;i < kMaxJets;++i){
 155     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
 156
 157     fh2RhoPtRec[i] = 0;
 158     fh2RhoPtGen[i] = 0;
 159     fh2PsiPtGen[i] = 0;
 160     fh2PsiPtRec[i] = 0;
 161     fh2FragRec[i] = 0;
 162     fh2FragLnRec[i] = 0;
 163     fh2FragGen[i] = 0;
 164     fh2FragLnGen[i] = 0;
 165   }
 166
 167
 168   for(int ij = 0;ij <kJetTypes;++ij){
 169     fh1NJets[ij] = 0;
 170     fh1SumPtTrack[ij] = 0;
 171     fh1PtJetsIn[ij] = 0;
 172     fh1PtTracksIn[ij] = 0;
 173     fh1PtTracksLeadingIn[ij] = 0;
 174     fh2NJetsPt[ij]  = 0;
 175     fh2NTracksPt[ij]  = 0;
 176     fh2LeadingJetPtJetPhi[ij] = 0;
 177     fh2LeadingTrackPtTrackPhi[ij] = 0;
 178     for(int i = 0;i < kMaxJets;++i){
 179       fh2PhiPt[ij][i] = 0;
 180       fh2EtaPt[ij][i] = 0;
 181       fh2PhiEta[ij][i] = 0;
 182
 183       fh1PtIn[ij][i] = 0;
 184       fh2RhoPt[ij][i] = 0;
 185       fh2PsiPt[ij][i] = 0;
 186     }
 187
 188     fh1DijetMinv[ij] = 0;
 189     fh2DijetDeltaPhiPt[ij] = 0;
 190     fh2DijetAsymPt[ij] = 0;
 191     fh2DijetPt2vsPt1[ij] = 0;
 192     fh2DijetDifvsSum[ij] = 0;
 193
 194   }
 195 }
 196
 197 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 198   AliAnalysisTaskSE(name),
 199   fJetHeaderRec(0x0),
 200   fJetHeaderGen(0x0),
 201   fAODIn(0x0),
 202   fAODOut(0x0),
 203   fAODExtension(0x0),
 204   fhnCorrelation(0x0),
 205   fhnCorrelationPhiZRec(0x0),
 206   f1PtScale(0x0),
 207   fBranchRec("jets"),
 208   fBranchGen(""),
 209   fBranchBkg(""),
 210   fNonStdFile(""),
 211   fUseAODJetInput(kFALSE),
 212   fUseAODTrackInput(kFALSE),
 213   fUseAODMCInput(kFALSE),
 214   fUseGlobalSelection(kFALSE),
 215   fUseExternalWeightOnly(kFALSE),
 216   fLimitGenJetEta(kFALSE),
 217   fBkgSubtraction(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     fh2RhoPtRec[i] = 0;
 290     fh2RhoPtGen[i] = 0;
 291     fh2PsiPtGen[i] = 0;
 292     fh2PsiPtRec[i] = 0;
 293     fh2FragRec[i] = 0;
 294     fh2FragLnRec[i] = 0;
 295     fh2FragGen[i] = 0;
 296     fh2FragLnGen[i] = 0;
 297   }
 298
 299   for(int ij = 0;ij <kJetTypes;++ij){
 300     fh1NJets[ij] = 0;
 301     fh1SumPtTrack[ij] = 0;
 302     fh1PtJetsIn[ij] = 0;
 303     fh1PtTracksIn[ij] = 0;
 304     fh1PtTracksLeadingIn[ij] = 0;
 305     fh2NJetsPt[ij]  = 0;
 306     fh2NTracksPt[ij]  = 0;
 307     fh2LeadingJetPtJetPhi[ij] = 0;
 308     fh2LeadingTrackPtTrackPhi[ij] = 0;
 309     for(int i = 0;i < kMaxJets;++i){
 310       fh2PhiPt[ij][i] = 0;
 311       fh2EtaPt[ij][i] = 0;
 312       fh2PhiEta[ij][i] = 0;
 313
 314       fh1PtIn[ij][i] = 0;
 315       fh2RhoPt[ij][i] = 0;
 316       fh2PsiPt[ij][i] = 0;
 317     }
 318
 319     fh1DijetMinv[ij] = 0;
 320     fh2DijetDeltaPhiPt[ij] = 0;
 321     fh2DijetAsymPt[ij] = 0;
 322     fh2DijetPt2vsPt1[ij] = 0;
 323     fh2DijetDifvsSum[ij] = 0;
 324   }
 325
 326  DefineOutput(1, TList::Class());
 327 }
 328
 329
 330
 331 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 332 {
 333
 334
 335
 336   //
 337   // Implemented Notify() to read the cross sections
 338   // and number of trials from pyxsec.root
 339   //
 340
 341   // Fetch the aod also from the input in,
 342   // have todo it in notify
 343
 344
 345   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 346   //  assume that the AOD is in the general output...
 347   fAODOut  = AODEvent();
 348
 349   if(fNonStdFile.Length()!=0){
 350     // case that we have an AOD extension we need can fetch the jets from the extended output
 351     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 352     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 353     if(!fAODExtension){
 354       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 355     }
 356   }
 357
 358
 359   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 360   Float_t xsection = 0;
 361   Float_t ftrials  = 1;
 362
 363   fAvgTrials = 1;
 364   if(tree){
 365     TFile *curfile = tree->GetCurrentFile();
 366     if (!curfile) {
 367       Error("Notify","No current file");
 368       return kFALSE;
 369     }
 370     if(!fh1Xsec||!fh1Trials){
 371       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 372       return kFALSE;
 373     }
 374     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 375     fh1Xsec->Fill("<#sigma>",xsection);
 376     // construct a poor man average trials
 377     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 378     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 379   }
 380
 381
 382   return kTRUE;
 383 }
 384
 385 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 386 {
 387
 388
 389   // Connect the AOD
 390
 391
 392
 393
 394
 395   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 396
 397   OpenFile(1);
 398   if(!fHistList)fHistList = new TList();
 399   fHistList->SetOwner(kTRUE);
 400   Bool_t oldStatus = TH1::AddDirectoryStatus();
 401   TH1::AddDirectory(kFALSE);
 402
 403
 404   MakeJetContainer();
 405   fHistList->Add(fhnCorrelation);
 406   fHistList->Add(fhnCorrelationPhiZRec);
 407   for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
 408
 409   //
 410   //  Histogram
 411
 412   const Int_t nBinPt = 320;
 413   Double_t binLimitsPt[nBinPt+1];
 414   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 415     if(iPt == 0){
 416       binLimitsPt[iPt] = 0.0;
 417     }
 418     else {// 1.0
 419       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 420     }
 421   }
 422
 423   const Int_t nBinPhi = 90;
 424   Double_t binLimitsPhi[nBinPhi+1];
 425   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 426     if(iPhi==0){
 427       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 428     }
 429     else{
 430       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 431     }
 432   }
 433
 434
 435   const Int_t nBinPhi2 = 360;
 436   Double_t binLimitsPhi2[nBinPhi2+1];
 437   for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
 438     if(iPhi2==0){
 439       binLimitsPhi2[iPhi2] = 0.;
 440     }
 441     else{
 442       binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
 443     }
 444   }
 445
 446   const Int_t nBinEta = 40;
 447   Double_t binLimitsEta[nBinEta+1];
 448   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 449     if(iEta==0){
 450       binLimitsEta[iEta] = -2.0;
 451     }
 452     else{
 453       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 454     }
 455   }
 456
 457
 458   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 459   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 460    fHistList->Add(fh1Xsec);
 461
 462   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 463   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 464   fHistList->Add(fh1Trials);
 465
 466   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 467   fHistList->Add(fh1PtHard);
 468
 469   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 470   fHistList->Add(fh1PtHardNoW);
 471
 472   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 473   fHistList->Add(fh1PtHardTrials);
 474
 475
 476   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 477   fHistList->Add(fh1ZVtx);
 478
 479   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 480   fHistList->Add(fh2PtFGen);
 481
 482   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
 483   fHistList->Add(fh2RelPtFGen);
 484
 485
 486     for(int ij = 0;ij <kJetTypes;++ij){
 487       TString cAdd = "";
 488       TString cJetBranch = "";
 489       if(ij==kJetRec){
 490         cAdd = "Rec";
 491         cJetBranch = fBranchRec.Data();
 492       }
 493       else if (ij==kJetGen){
 494         cAdd = "Gen";
 495         cJetBranch = fBranchGen.Data();
 496       }
 497
 498       fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
 499       fHistList->Add(fh1NJets[ij]);
 500
 501       fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 502       fHistList->Add(fh1PtJetsIn[ij]);
 503
 504       fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 505       fHistList->Add(fh1PtTracksIn[ij]);
 506
 507       fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 508       fHistList->Add(fh1PtTracksLeadingIn[ij]);
 509
 510       fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 511       fHistList->Add(fh1SumPtTrack[ij]);
 512
 513       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);
 514       fHistList->Add(fh2NJetsPt[ij]);
 515
 516       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);
 517       fHistList->Add(fh2NTracksPt[ij]);
 518
 519       fh2LeadingJetPtJetPhi[ij] = new TH2F(Form("fh2Leading%sJetPtJetPhi",cAdd.Data()),Form("phi of leading %s jet;p_{T};#phi;",cAdd.Data()),
 520                                            nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
 521       fHistList->Add(fh2LeadingJetPtJetPhi[ij]);
 522
 523       fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
 524                                                nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
 525       fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
 526
 527       for(int i = 0;i < kMaxJets;++i){
 528
 529         fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
 530         fHistList->Add(fh1PtIn[ij][i]);
 531
 532         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()),
 533                                    40,0.,2.,nBinPt,binLimitsPt);
 534         fHistList->Add(fh2RhoPt[ij][i]);
 535         if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 536
 537
 538         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()),
 539                                      40,0.,2.,nBinPt,binLimitsPt);
 540         fHistList->Add(fh2PsiPt[ij][i]);
 541
 542         fh2PhiPt[ij][i] = new TH2F(Form("fh2PhiPt%s_j%d",cAdd.Data(),i),Form("pt vs phi %s;phi;p_{T};",cAdd.Data()),
 543                                    nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 544         fHistList->Add(fh2PhiPt[ij][i]);
 545         fh2EtaPt[ij][i] = new TH2F(Form("fh2EtaPt%s_j%d",cAdd.Data(),i),Form("pt vs phi %s;phi;p_{T};",cAdd.Data()),
 546                                    20,-1.,1.,nBinPt,binLimitsPt);
 547         fHistList->Add(fh2EtaPt[ij][i]);
 548         fh2PhiEta[ij][i] =  new TH2F(Form("fh2PhiEta%s_j%d",cAdd.Data(),i),Form("phi vs eta %s ;phi;p_{T};",cAdd.Data()),
 549                                      nBinPhi,binLimitsPhi,20,-1.,1.);
 550         fHistList->Add(fh2PhiEta[ij][i]);
 551
 552       }
 553
 554
 555       fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 556       fHistList->Add(fh1DijetMinv[ij]);
 557
 558       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);
 559       fHistList->Add(fh2DijetDeltaPhiPt[ij]);
 560
 561       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);
 562       fHistList->Add(fh2DijetAsymPt[ij]);
 563
 564       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.);
 565       fHistList->Add(fh2DijetPt2vsPt1[ij]);
 566
 567       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.);
 568       fHistList->Add( fh2DijetDifvsSum[ij]);
 569     }
 570
 571   // =========== Switch on Sumw2 for all histos ===========
 572   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 573     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 574     if (h1){
 575       h1->Sumw2();
 576       continue;
 577     }
 578     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 579     if(hn)hn->Sumw2();
 580   }
 581   TH1::AddDirectory(oldStatus);
 582 }
 583
 584 void AliAnalysisTaskJetSpectrum2::Init()
 585 {
 586   //
 587   // Initialization
 588   //
 589
 590   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 591
 592 }
 593
 594 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
 595
 596
 597   Bool_t selected = kTRUE;
 598
 599   if(fUseGlobalSelection&&fEventSelectionMask==0){
 600     selected = AliAnalysisHelperJetTasks::Selected();
 601   }
 602   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 603     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 604   }
 605
 606   if(fEventClass>0){
 607     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 608   }
 609
 610   if(!selected){
 611     // no selection by the service task, we continue
 612     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 613     PostData(1, fHistList);
 614     return;
 615   }
 616
 617
 618   static AliAODEvent* aod = 0;
 619
 620   // take all other information from the aod we take the tracks from
 621   if(!aod){
 622    if(fUseAODTrackInput)aod = fAODIn;
 623    else aod = fAODOut;
 624   }
 625
 626
 627   //
 628   // Execute analysis for current event
 629   //
 630   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 631   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 632   if(!aodH){
 633     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 634     return;
 635   }
 636
 637   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 638   TClonesArray *aodRecJets = 0;
 639
 640   if(fAODOut&&!aodRecJets){
 641     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
 642   }
 643   if(fAODExtension&&!aodRecJets){
 644     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
 645   }
 646   if(fAODIn&&!aodRecJets){
 647     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
 648   }
 649
 650
 651
 652   if(!aodRecJets){
 653     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 654     return;
 655   }
 656
 657   TClonesArray *aodGenJets = 0;
 658   if(fBranchGen.Length()>0){
 659     if(fAODOut&&!aodGenJets){
 660       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
 661     }
 662     if(fAODExtension&&!aodGenJets){
 663       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
 664     }
 665     if(fAODIn&&!aodGenJets){
 666       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
 667     }
 668
 669     if(!aodGenJets){
 670       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
 671       return;
 672     }
 673   }
 674
 675
 676   // new Scheme
 677   // first fill all the pure  histograms, i.e. full jets
 678   // in case of the correaltion limit to the n-leading jets
 679
 680   // reconstructed
 681
 682
 683   // generated
 684
 685
 686   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
 687
 688
 689   TList genJetsList;         // full acceptance
 690   TList genJetsListCut;      // acceptance cut
 691   TList recJetsList;         // full acceptance
 692   TList recJetsListCut;      // acceptance cut
 693
 694
 695   GetListOfJets(&recJetsList,aodRecJets,0);
 696   GetListOfJets(&recJetsListCut,aodRecJets,1);
 697
 698   if(fBranchGen.Length()>0){
 699     GetListOfJets(&genJetsList,aodGenJets,0);
 700     GetListOfJets(&genJetsListCut,aodGenJets,1);
 701   }
 702
 703   Double_t eventW = 1;
 704   Double_t ptHard = 0;
 705   Double_t nTrials = 1; // Trials for MC trigger
 706   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 707
 708   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 709     // this is the part we only use when we have MC information
 710     AliMCEvent* mcEvent = MCEvent();
 711     //    AliStack *pStack = 0;
 712     if(!mcEvent){
 713       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 714       return;
 715     }
 716     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 717     if(pythiaGenHeader){
 718       nTrials = pythiaGenHeader->Trials();
 719       ptHard  = pythiaGenHeader->GetPtHard();
 720       int iProcessType = pythiaGenHeader->ProcessType();
 721       // 11 f+f -> f+f
 722       // 12 f+barf -> f+barf
 723       // 13 f+barf -> g+g
 724       // 28 f+g -> f+g
 725       // 53 g+g -> f+barf
 726       // 68 g+g -> g+g
 727       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 728       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 729     }
 730   }// (fAnalysisType&kMCESD)==kMCESD)
 731
 732
 733   // we simply fetch the tracks/mc particles as a list of AliVParticles
 734
 735   TList recParticles;
 736   TList genParticles;
 737
 738   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 739   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 740   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 741   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 742
 743   // Event control and counting ...
 744   // MC
 745   fh1PtHard->Fill(ptHard,eventW);
 746   fh1PtHardNoW->Fill(ptHard,1);
 747   fh1PtHardTrials->Fill(ptHard,nTrials);
 748
 749   // Real
 750   if(aod->GetPrimaryVertex()){// No vtx for pure MC
 751     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
 752   }
 753
 754   // the loops for rec and gen should be indentical... pass it to a separate
 755   // function ...
 756   // Jet Loop
 757   // Track Loop
 758   // Jet Jet Loop
 759   // Jet Track loop
 760
 761   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 762   FillTrackHistos(recParticles,kJetRec);
 763
 764   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 765   FillTrackHistos(genParticles,kJetGen);
 766
 767   // Here follows the jet matching part
 768   // delegate to a separated method?
 769
 770   FillMatchHistos(recJetsList,genJetsList);
 771
 772   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 773   PostData(1, fHistList);
 774 }
 775
 776 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 777
 778   if(iType>=kJetTypes){
 779     return;
 780   }
 781
 782   Int_t nJets = jetsList.GetEntries();
 783   fh1NJets[iType]->Fill(nJets);
 784
 785   if(nJets<=0)return;
 786
 787   Float_t ptOld = 1.E+32;
 788   Float_t pT0 = 0;
 789   Float_t pT1 = 0;
 790   Float_t phi0 = 0;
 791   Float_t phi1 = 0;
 792   Int_t ij0 = -1;
 793   Int_t ij1 = -1;
 794
 795
 796   for(int ij = 0;ij < nJets;ij++){
 797     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
 798     Float_t ptJet = jet->Pt();
 799     fh1PtJetsIn[iType]->Fill(ptJet);
 800     if(ptJet>ptOld){
 801       Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
 802     }
 803     ptOld = ptJet;
 804
 805     // find the dijets assume sorting and acceptance cut...
 806     if(ij==0){
 807       ij0 = ij;
 808       pT0 = ptJet;
 809       phi0 = jet->Phi();
 810       if(phi0<0)phi0+=TMath::Pi()*2.;
 811     }
 812     else if(ptJet>pT1){
 813       // take only the backward hemisphere??
 814       phi1 = jet->Phi();
 815       if(phi1<0)phi1+=TMath::Pi()*2.;
 816       Float_t dPhi = phi1 - phi0;
 817       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 818       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 819       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
 820         ij1 = ij;
 821         pT1 = ptJet;
 822       }
 823     }
 824     // fill jet histos for kmax jets
 825     if(ij<kMaxJets){
 826       Float_t phiJet = jet->Phi();
 827       Float_t etaJet = jet->Eta();
 828       if(phiJet<0)phiJet+=TMath::Pi()*2.;
 829       fh1TmpRho->Reset();
 830       fh1PtIn[iType][ij]->Fill(ptJet);
 831       // fill leading jets...
 832       fh2PhiEta[iType][ij]->Fill(phiJet,etaJet);
 833       fh2PhiPt[iType][ij]->Fill(phiJet,ptJet);
 834       if(ptJet>10)fh2EtaPt[iType][ij]->Fill(etaJet,ptJet);
 835       if(ij==0){
 836         fh2LeadingJetPtJetPhi[iType]->Fill(ptJet,phiJet);
 837         if(ij==0&&iType==0&&fDebug>1){
 838           Printf("%d %3.3f %p %s",iType,ptJet,jet,fBranchRec.Data());
 839         }
 840       }
 841       if(particlesList.GetSize()){
 842
 843         // Particles... correlated with jets...
 844         for(int it = 0;it<particlesList.GetEntries();++it){
 845           AliVParticle *part = (AliVParticle*)particlesList.At(it);
 846           Float_t deltaR = jet->DeltaR(part);
 847           if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
 848         }
 849         // fill the jet shapes
 850         Float_t rhoSum = 0;
 851         for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
 852           Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
 853           Float_t rho = fh1TmpRho->GetBinContent(ibx);
 854           rhoSum += rho;
 855           fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
 856           fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
 857         }
 858       }// if we have particles
 859     }// ij < kMaxJets
 860   }// Jet Loop
 861
 862
 863   // do something with dijets...
 864   if(ij0>=0&&ij1>0){
 865     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
 866     Double_t ptJet0 = jet0->Pt();
 867     Double_t phiJet0 = jet0->Phi();
 868     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
 869
 870     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
 871     Double_t ptJet1 = jet1->Pt();
 872     Double_t phiJet1 = jet1->Phi();
 873     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
 874
 875     Float_t deltaPhi = phiJet0 - phiJet1;
 876     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
 877     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
 878     deltaPhi = TMath::Abs(deltaPhi);
 879     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
 880
 881     Float_t asym = 9999;
 882     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
 883       fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
 884       fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
 885       fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
 886       Float_t minv = 2.*(jet0->P()*jet1->P()-
 887                          jet0->Px()*jet1->Px()-
 888                          jet0->Py()*jet1->Py()-
 889                          jet0->Pz()*jet1->Pz());    // assume mass == 0;
 890       if(minv<0)minv=0; // prevent numerical instabilities
 891       minv = TMath::Sqrt(minv);
 892       fh1DijetMinv[iType]->Fill(minv);
 893   }
 894
 895
 896
 897   // count down the jets above thrueshold
 898   Int_t nOver = nJets;
 899   if(nOver>0){
 900     TIterator *jetIter = jetsList.MakeIterator();
 901     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
 902     if(tmpJet){
 903       Float_t pt = tmpJet->Pt();
 904       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
 905         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
 906         while(pt<ptCut&&tmpJet){
 907           nOver--;
 908           tmpJet = (AliAODJet*)(jetIter->Next());
 909           if(tmpJet){
 910             pt = tmpJet->Pt();
 911           }
 912         }
 913         if(nOver<=0)break;
 914         fh2NJetsPt[iType]->Fill(ptCut,nOver);
 915       }
 916     }
 917     delete jetIter;
 918   }
 919 }
 920
 921 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
 922
 923   Int_t nTrackOver = particlesList.GetSize();
 924   // do the same for tracks and jets
 925   if(nTrackOver>0){
 926     TIterator *trackIter = particlesList.MakeIterator();
 927     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
 928     Float_t pt = tmpTrack->Pt();
 929     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
 930       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
 931       while(pt<ptCut&&tmpTrack){
 932         nTrackOver--;
 933         tmpTrack = (AliVParticle*)(trackIter->Next());
 934         if(tmpTrack){
 935           pt = tmpTrack->Pt();
 936         }
 937       }
 938       if(nTrackOver<=0)break;
 939       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
 940     }
 941
 942
 943     trackIter->Reset();
 944     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
 945     Float_t sumPt = 0;
 946
 947     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
 948       Float_t tmpPt = tmpTrack->Pt();
 949       fh1PtTracksIn[iType]->Fill(tmpPt);
 950       sumPt += tmpPt;
 951       Float_t tmpPhi = tmpTrack->Phi();
 952       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 953       if(tmpTrack==leading){
 954         fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
 955         fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
 956         continue;
 957       }
 958     }
 959     fh1SumPtTrack[iType]->Fill(sumPt);
 960     delete trackIter;
 961   }
 962
 963 }
 964
 965
 966 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
 967
 968
 969   // Fill al the matching histos
 970   // It is important that the acceptances for the mathing are as large as possible
 971   // to avoid false matches on the edge of acceptance
 972   // therefore we add some extra matching jets as overhead
 973
 974   static TArrayI aGenIndex(recJetsList.GetEntries());
 975   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
 976
 977   static TArrayI aRecIndex(genJetsList.GetEntries());
 978   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
 979
 980   if(fDebug){
 981     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
 982     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
 983   }
 984   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
 985                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
 986                                             aGenIndex,aRecIndex,fDebug);
 987
 988   if(fDebug){
 989     for(int i = 0;i< aGenIndex.GetSize();++i){
 990       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
 991     }
 992     for(int i = 0;i< aRecIndex.GetSize();++i){
 993       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
 994     }
 995   }
 996
 997   Double_t container[6];
 998   Double_t containerPhiZ[6];
 999
1000   // loop over generated jets
1001   // consider the
1002   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1003     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1004     Double_t ptGen = genJet->Pt();
1005     Double_t phiGen = genJet->Phi();
1006     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1007     Double_t etaGen = genJet->Eta();
1008     container[3] = ptGen;
1009     container[4] = etaGen;
1010     container[5] = phiGen;
1011     fhnJetContainer[kStep0]->Fill(&container[3]);
1012     if(JetSelected(genJet)){
1013       fhnJetContainer[kStep1]->Fill(&container[3]);
1014       Int_t ir = aRecIndex[ig];
1015       if(ir>=0&&ir<recJetsList.GetEntries()){
1016         fhnJetContainer[kStep2]->Fill(&container[3]);
1017         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1018         if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
1019         if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
1020       }
1021     }
1022   }// loop over generated jets used for matching...
1023
1024
1025
1026   // loop over reconstructed jets
1027   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1028     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1029     Double_t etaRec = recJet->Eta();
1030     Double_t ptRec = recJet->Pt();
1031     Double_t phiRec = recJet->Phi();
1032     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1033     // do something with dijets...
1034
1035     container[0] = ptRec;
1036     container[1] = etaRec;
1037     container[2] = phiRec;
1038     containerPhiZ[0] = ptRec;
1039     containerPhiZ[1] = phiRec;
1040
1041     fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1042     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1043
1044     if(JetSelected(recJet)){
1045       fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1046       // Fill Correlation
1047       Int_t ig = aGenIndex[ir];
1048       if(ig>=0 && ig<genJetsList.GetEntries()){
1049         fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1050         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1051         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1052         Double_t ptGen  = genJet->Pt();
1053         Double_t phiGen = genJet->Phi();
1054         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1055         Double_t etaGen = genJet->Eta();
1056
1057         container[3] = ptGen;
1058         container[4] = etaGen;
1059         container[5] = phiGen;
1060         containerPhiZ[3] = ptGen;
1061         //
1062         // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1063         //
1064         if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1065         fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1066         fhnCorrelation->Fill(container);
1067         if(ptGen>0){
1068           Float_t delta = (ptRec-ptGen)/ptGen;
1069           fh2RelPtFGen->Fill(ptGen,delta);
1070           fh2PtFGen->Fill(ptGen,ptRec);
1071         }
1072         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1073       }
1074       else{
1075         containerPhiZ[3] = 0;
1076         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1077       }
1078     }// loop over reconstructed jets
1079   }
1080   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1081 }
1082
1083
1084 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1085   //
1086   // Create the particle container for the correction framework manager and
1087   // link it
1088   //
1089   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1090   const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1091   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1092   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1093   const Double_t kZmin = 0., kZmax = 1;
1094
1095   // can we neglect migration in eta and phi?
1096   // phi should be no problem since we cover full phi and are phi symmetric
1097   // eta migration is more difficult  due to needed acceptance correction
1098   // in limited eta range
1099
1100   //arrays for the number of bins in each dimension
1101   Int_t iBin[kNvar];
1102   iBin[0] = 320; //bins in pt
1103   iBin[1] =  1; //bins in eta
1104   iBin[2] = 1; // bins in phi
1105
1106   //arrays for lower bounds :
1107   Double_t* binEdges[kNvar];
1108   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1109     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1110
1111   //values for bin lower bounds
1112   //  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);
1113   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1114   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1115   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1116
1117
1118   for(int i = 0;i < kMaxStep*2;++i){
1119     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1120     for (int k=0; k<kNvar; k++) {
1121       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1122     }
1123   }
1124   //create correlation matrix for unfolding
1125   Int_t thnDim[2*kNvar];
1126   for (int k=0; k<kNvar; k++) {
1127     //first half  : reconstructed
1128     //second half : MC
1129     thnDim[k]      = iBin[k];
1130     thnDim[k+kNvar] = iBin[k];
1131   }
1132
1133   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1134   for (int k=0; k<kNvar; k++) {
1135     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1136     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1137   }
1138   fhnCorrelation->Sumw2();
1139
1140   // for second correlation histogram
1141
1142
1143   const Int_t kNvarPhiZ   = 4;
1144   //arrays for the number of bins in each dimension
1145   Int_t iBinPhiZ[kNvarPhiZ];
1146   iBinPhiZ[0] = 80; //bins in pt
1147   iBinPhiZ[1] = 72; //bins in phi
1148   iBinPhiZ[2] = 20; // bins in Z
1149   iBinPhiZ[3] = 80; //bins in ptgen
1150
1151   //arrays for lower bounds :
1152   Double_t* binEdgesPhiZ[kNvarPhiZ];
1153   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1154     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1155
1156   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1157   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1158   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1159   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1160
1161   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1162   for (int k=0; k<kNvarPhiZ; k++) {
1163     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1164   }
1165   fhnCorrelationPhiZRec->Sumw2();
1166
1167
1168   // Add a histogram for Fake jets
1169
1170   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1171     delete [] binEdges[ivar];
1172
1173   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1174     delete [] binEdgesPhiZ[ivar];
1175
1176 }
1177
1178 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1179 {
1180 // Terminate analysis
1181 //
1182     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1183 }
1184
1185
1186 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1187
1188   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1189
1190   Int_t iCount = 0;
1191   if(type==kTrackAOD){
1192     AliAODEvent *aod = 0;
1193     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1194     else aod = AODEvent();
1195     if(!aod){
1196       return iCount;
1197     }
1198     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1199       AliAODTrack *tr = aod->GetTrack(it);
1200       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1201       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1202       if(tr->Pt()<fMinTrackPt)continue;
1203       if(fDebug>0){
1204         if(tr->Pt()>20){
1205           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1206           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1207           tr->Print();
1208           //    tr->Dump();
1209           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1210           if(fESD){
1211             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1212             esdTr->Print("");
1213             // esdTr->Dump();
1214           }
1215         }
1216       }
1217       list->Add(tr);
1218       iCount++;
1219     }
1220   }
1221   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1222     AliMCEvent* mcEvent = MCEvent();
1223     if(!mcEvent)return iCount;
1224     // we want to have alivpartilces so use get track
1225     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1226       if(!mcEvent->IsPhysicalPrimary(it))continue;
1227       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1228       if(part->Pt()<fMinTrackPt)continue;
1229       if(type == kTrackKineAll){
1230         list->Add(part);
1231         iCount++;
1232       }
1233       else if(type == kTrackKineCharged){
1234         if(part->Particle()->GetPDG()->Charge()==0)continue;
1235         list->Add(part);
1236         iCount++;
1237       }
1238     }
1239   }
1240   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1241     AliAODEvent *aod = 0;
1242     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1243     else aod = AODEvent();
1244     if(!aod)return iCount;
1245     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1246     if(!tca)return iCount;
1247     for(int it = 0;it < tca->GetEntriesFast();++it){
1248       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1249       if(!part)continue;
1250       if(part->Pt()<fMinTrackPt)continue;
1251       if(!part->IsPhysicalPrimary())continue;
1252       if(type == kTrackAODMCAll){
1253         list->Add(part);
1254         iCount++;
1255       }
1256       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1257         if(part->Charge()==0)continue;
1258         if(kTrackAODMCCharged){
1259           list->Add(part);
1260         }
1261         else {
1262           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1263           list->Add(part);
1264         }
1265         iCount++;
1266       }
1267     }
1268   }// AODMCparticle
1269   list->Sort();
1270   return iCount;
1271
1272 }
1273
1274
1275
1276 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1277   Bool_t selected = false;
1278
1279   if(!jet)return selected;
1280
1281   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1282     selected = kTRUE;
1283   }
1284   return selected;
1285
1286 }
1287
1288 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1289
1290   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1291   Int_t iCount = 0;
1292
1293   if(!jarray){
1294     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1295     return iCount;
1296   }
1297
1298
1299   for(int ij=0;ij<jarray->GetEntries();ij++){
1300     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1301     if(!jet)continue;
1302     if(type==0){
1303       // No cut at all, main purpose here is sorting
1304       list->Add(jet);
1305       iCount++;
1306     }
1307     else if(type == 1){
1308       // eta cut
1309       if(JetSelected(jet)){
1310         list->Add(jet);
1311         iCount++;
1312       }
1313     }
1314   }
1315
1316   list->Sort();
1317   return iCount;
1318
1319 }
1320
1321