PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx

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