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                                 40,0.,1.,nBinPt,binLimitsPt);
 427     fh2PsiPtRec[ij] =  new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
 428                                 40,0.,2.,nBinPt,binLimitsPt);
 429
 430     fh2RhoPtGen[ij] =  new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
 431                                40,0.,2.,nBinPt,binLimitsPt);
 432     fh2PsiPtGen[ij] =  new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
 433                                 40,0.,2.,nBinPt,binLimitsPt);
 434     if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 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("fh2DijetDeltaPhiPt","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     fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 624   }
 625
 626
 627
 628
 629   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 630
 631
 632   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 633
 634   if(!aodH){
 635     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 636     return;
 637   }
 638
 639   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 640   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 641   if(!aodRecJets){
 642     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 643     return;
 644   }
 645
 646   if(fBkgSubtraction){
 647      AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
 648
 649
 650      if(!evBkg){
 651        Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
 652        return;
 653      }
 654
 655
 656      ///just to start: some very simple plots containing rho, sigma and area of the background.
 657
 658      Int_t nJets = aodRecJets->GetEntriesFast();
 659      Float_t pthardest=0;
 660      if(nJets!=0){
 661
 662
 663        Float_t bkg1=evBkg->GetBackground(0);
 664        Float_t bkg2=evBkg->GetBackground(1);
 665        Float_t sigma1=evBkg->GetSigma(0);
 666        Float_t sigma2=evBkg->GetSigma(1);
 667        Float_t area1=evBkg->GetMeanarea(0);
 668        Float_t area2=evBkg->GetMeanarea(1);
 669        fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
 670        fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the hardest.
 671        fh1Sigma1->Fill(sigma1);
 672        fh1Sigma2->Fill(sigma2);
 673        fh1Area1->Fill(area1);
 674        fh1Area2->Fill(area2);
 675
 676
 677
 678        for(Int_t k=0;k<nJets;k++){
 679          AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
 680          fh1Ptjet->Fill(jet->Pt());
 681          Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
 682          Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
 683          Float_t err1=sigma1*sqrt(area1);
 684          Float_t err2=sigma2*sqrt(area2);
 685
 686          fh1Ptjetsub1->Fill(ptsub1);
 687          fh1Ptjetsub2->Fill(ptsub2);
 688
 689          if(k==0) {
 690            pthardest=jet->Pt();
 691            fh1Ptjethardest->Fill(pthardest);
 692            fh1Ptjetsubhardest1->Fill(ptsub1);
 693            fh1Ptjetsubhardest2->Fill(ptsub2);
 694            fh1Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
 695            fh1Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
 696          }
 697        }
 698        fh1Rhovspthardest1->Fill(pthardest,bkg1);
 699        fh1Rhovspthardest2->Fill(pthardest,bkg2);
 700      }
 701   }// background subtraction
 702
 703
 704   // ==== General variables needed
 705
 706
 707   // We use statice array, not to fragment the memory
 708   AliAODJet genJets[kMaxJets];
 709   Int_t nGenJets = 0;
 710   AliAODJet recJets[kMaxJets];
 711   Int_t nRecJets = 0;
 712   ///////////////////////////
 713
 714
 715   Double_t eventW = 1;
 716   Double_t ptHard = 0;
 717   Double_t nTrials = 1; // Trials for MC trigger
 718
 719   if(fUseExternalWeightOnly){
 720     eventW = fExternalWeight;
 721   }
 722
 723   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 724   //  if(fDebug>0)aodH->SetFillAOD(kFALSE);
 725   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 726   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 727     // this is the part we only use when we have MC information
 728     AliMCEvent* mcEvent = MCEvent();
 729     //    AliStack *pStack = 0;
 730     if(!mcEvent){
 731       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 732       return;
 733     }
 734     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 735     Int_t iCount = 0;
 736     if(pythiaGenHeader){
 737       nTrials = pythiaGenHeader->Trials();
 738       ptHard  = pythiaGenHeader->GetPtHard();
 739       int iProcessType = pythiaGenHeader->ProcessType();
 740       // 11 f+f -> f+f
 741       // 12 f+barf -> f+barf
 742       // 13 f+barf -> g+g
 743       // 28 f+g -> f+g
 744       // 53 g+g -> f+barf
 745       // 68 g+g -> g+g
 746       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 747       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 748
 749       // fetch the pythia generated jets only to be used here
 750       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 751       AliAODJet pythiaGenJets[kMaxJets];
 752       for(int ip = 0;ip < nPythiaGenJets;++ip){
 753         if(iCount>=kMaxJets)continue;
 754         Float_t p[4];
 755         pythiaGenHeader->TriggerJet(ip,p);
 756         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 757
 758         if(fBranchGen.Length()==0){
 759           /*
 760           if(fLimitGenJetEta){
 761             if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
 762                pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 763           }
 764           */
 765           if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
 766           // if we have MC particles and we do not read from the aod branch
 767           // use the pythia jets
 768           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 769           iCount++;
 770         }
 771       }
 772     }
 773     if(fBranchGen.Length()==0)nGenJets = iCount;
 774   }// (fAnalysisType&kMCESD)==kMCESD)
 775
 776
 777   // we simply fetch the tracks/mc particles as a list of AliVParticles
 778
 779   TList recParticles;
 780   TList genParticles;
 781
 782
 783
 784
 785   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 786   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 787   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 788   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 789
 790
 791   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 792   fh1PtHard->Fill(ptHard,eventW);
 793   fh1PtHardNoW->Fill(ptHard,1);
 794   fh1PtHardTrials->Fill(ptHard,nTrials);
 795
 796   // If we set a second branch for the input jets fetch this
 797   if(fBranchGen.Length()>0){
 798     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 799     if(aodGenJets){
 800       Int_t iCount = 0;
 801       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
 802         if(iCount>=kMaxJets)continue;
 803         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 804         if(!tmp)continue;
 805         /*
 806         if(fLimitGenJetEta){
 807           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
 808              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
 809         }
 810         */
 811         if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
 812         genJets[iCount] = *tmp;
 813         iCount++;
 814       }
 815       nGenJets = iCount;
 816     }
 817     else{
 818       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
 819       if(fDebug>2)fAOD->Print();
 820     }
 821   }
 822
 823   fh1NGenJets->Fill(nGenJets);
 824   // We do not want to exceed the maximum jet number
 825   nGenJets = TMath::Min(nGenJets,kMaxJets);
 826
 827   // Fetch the reconstructed jets...
 828
 829   nRecJets = aodRecJets->GetEntries();
 830
 831   nRecJets = aodRecJets->GetEntries();
 832   fh1NRecJets->Fill(nRecJets);
 833
 834   // Do something with the all rec jets
 835   Int_t nRecOver = nRecJets;
 836
 837   // check that the jets are sorted
 838   Float_t ptOld = 999999;
 839   for(int ir = 0;ir < nRecJets;ir++){
 840     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
 841     Float_t tmpPt = tmp->Pt();
 842     if(tmpPt>ptOld){
 843       Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
 844     }
 845     ptOld = tmpPt;
 846   }
 847
 848
 849   if(nRecOver>0){
 850     TIterator *recIter = aodRecJets->MakeIterator();
 851     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
 852     Float_t pt = tmpRec->Pt();
 853     if(tmpRec){
 854       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 855         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 856         while(pt<ptCut&&tmpRec){
 857           nRecOver--;
 858           tmpRec = (AliAODJet*)(recIter->Next());
 859           if(tmpRec){
 860             pt = tmpRec->Pt();
 861           }
 862         }
 863         if(nRecOver<=0)break;
 864         fh2NRecJetsPt->Fill(ptCut,nRecOver);
 865       }
 866     }
 867     recIter->Reset();
 868     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
 869     Float_t phi = leading->Phi();
 870     if(phi<0)phi+=TMath::Pi()*2.;
 871     Float_t eta = leading->Eta();
 872     pt = leading->Pt();
 873     while((tmpRec = (AliAODJet*)(recIter->Next()))){
 874       Float_t tmpPt = tmpRec->Pt();
 875       fh1PtJetsRecIn->Fill(tmpPt);
 876       if(tmpRec==leading){
 877         fh1PtJetsLeadingRecIn->Fill(tmpPt);
 878         continue;
 879       }
 880       // correlation
 881       Float_t tmpPhi =  tmpRec->Phi();
 882
 883       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 884       Float_t dPhi = phi - tmpRec->Phi();
 885       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 886       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 887       Float_t dEta = eta - tmpRec->Eta();
 888       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 889       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 890     }
 891     delete recIter;
 892   }
 893
 894   Int_t nTrackOver = recParticles.GetSize();
 895   // do the same for tracks and jets
 896   if(nTrackOver>0){
 897     TIterator *recIter = recParticles.MakeIterator();
 898     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
 899     Float_t pt = tmpRec->Pt();
 900     //    Printf("Leading track p_t %3.3E",pt);
 901     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
 902       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
 903       while(pt<ptCut&&tmpRec){
 904         nTrackOver--;
 905         tmpRec = (AliVParticle*)(recIter->Next());
 906         if(tmpRec){
 907           pt = tmpRec->Pt();
 908         }
 909       }
 910       if(nTrackOver<=0)break;
 911       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
 912     }
 913
 914     recIter->Reset();
 915     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
 916     Float_t phi = leading->Phi();
 917     if(phi<0)phi+=TMath::Pi()*2.;
 918     Float_t eta = leading->Eta();
 919     pt  = leading->Pt();
 920     while((tmpRec = (AliVParticle*)(recIter->Next()))){
 921       Float_t tmpPt = tmpRec->Pt();
 922       fh1PtTracksRecIn->Fill(tmpPt);
 923       if(tmpRec==leading){
 924         fh1PtTracksLeadingRecIn->Fill(tmpPt);
 925         continue;
 926       }
 927       // correlation
 928       Float_t tmpPhi =  tmpRec->Phi();
 929
 930       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 931       Float_t dPhi = phi - tmpRec->Phi();
 932       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 933       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 934       Float_t dEta = eta - tmpRec->Eta();
 935       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
 936       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
 937     }
 938     delete recIter;
 939   }
 940
 941   if(genParticles.GetSize()){
 942     TIterator *genIter = genParticles.MakeIterator();
 943     AliVParticle *tmpGen = 0;
 944     while((tmpGen = (AliVParticle*)(genIter->Next()))){
 945       if(TMath::Abs(tmpGen->Eta())<0.9){
 946         Float_t tmpPt = tmpGen->Pt();
 947         fh1PtTracksGenIn->Fill(tmpPt);
 948       }
 949     }
 950     delete genIter;
 951   }
 952
 953   nRecJets = TMath::Min(nRecJets,kMaxJets);
 954
 955   Int_t iCountRec = 0;
 956   for(int ir = 0;ir < nRecJets;++ir){
 957     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 958     if(!tmp)continue;
 959     if(tmp->Pt()<fMinJetPt)continue;
 960     recJets[ir] = *tmp;
 961     iCountRec++;
 962   }
 963   nRecJets = iCountRec;
 964
 965
 966   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 967   // Relate the jets
 968   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
 969   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
 970
 971
 972   for(int i = 0;i<kMaxJets;++i){
 973     iGenIndex[i] = iRecIndex[i] = -1;
 974   }
 975
 976   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
 977                                             iGenIndex,iRecIndex,fDebug);
 978   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 979
 980   if(fDebug){
 981     for(int i = 0;i<kMaxJets;++i){
 982       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
 983       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
 984     }
 985   }
 986
 987
 988
 989
 990   Double_t container[6];
 991   Double_t containerPhiZ[6];
 992
 993   // loop over generated jets
 994
 995   // radius; tmp, get from the jet header itself
 996   // at some point, how todeal woht FastJet on MC?
 997   Float_t radiusGen =0.4;
 998   Float_t radiusRec =0.4;
 999
1000   for(int ig = 0;ig < nGenJets;++ig){
1001     Double_t ptGen = genJets[ig].Pt();
1002     Double_t phiGen = genJets[ig].Phi();
1003     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1004     Double_t etaGen = genJets[ig].Eta();
1005
1006     container[3] = ptGen;
1007     container[4] = etaGen;
1008     container[5] = phiGen;
1009     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1010     Int_t ir = iRecIndex[ig];
1011
1012     if(TMath::Abs(etaGen)<fRecEtaWindow){
1013       fh1TmpRho->Reset();
1014
1015       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1016       fh1PtGenIn[ig]->Fill(ptGen,eventW);
1017       // fill the fragmentation function
1018       for(int it = 0;it<genParticles.GetEntries();++it){
1019         AliVParticle *part = (AliVParticle*)genParticles.At(it);
1020         Float_t deltaR = genJets[ig].DeltaR(part);
1021         fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1022         if(deltaR<radiusGen){
1023           Float_t z = part->Pt()/ptGen;
1024           Float_t lnz =  -1.*TMath::Log(z);
1025           fh2FragGen[ig]->Fill(z,ptGen,eventW);
1026           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1027         }
1028
1029       }
1030       Float_t rhoSum = 0;
1031       for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
1032         Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1033         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1034         rhoSum += rho;
1035         fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1036         fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
1037       }
1038     }
1039     if(ir>=0&&ir<nRecJets){
1040       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1041       Double_t etaRec = recJets[ir].Eta();
1042       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1043     }
1044   }// loop over generated jets
1045
1046   Float_t sumPt = 0;
1047   for(int it = 0;it<recParticles.GetEntries();++it){
1048     AliVParticle *part = (AliVParticle*)recParticles.At(it);
1049     // fill sum pt and P_t of all paritles
1050     if(TMath::Abs(part->Eta())<0.9){
1051       Float_t pt = part->Pt();
1052       fh1PtTrackRec->Fill(pt,eventW);
1053       fh2TrackPtTrackPhi->Fill(pt,part->Phi());
1054       sumPt += pt;
1055     }
1056   }
1057   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1058   fh1SumPtTrackRec->Fill(sumPt,eventW);
1059
1060
1061   // loop over reconstructed jets
1062   for(int ir = 0;ir < nRecJets;++ir){
1063     Double_t etaRec = recJets[ir].Eta();
1064     Double_t ptRec = recJets[ir].Pt();
1065     Double_t phiRec = recJets[ir].Phi();
1066     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1067
1068
1069   // do something with dijets...
1070   if(ir==1){
1071     Double_t etaRec1 = recJets[0].Eta();
1072     Double_t ptRec1 = recJets[0].Pt();
1073     Double_t phiRec1 = recJets[0].Phi();
1074     if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
1075
1076
1077     if(TMath::Abs(etaRec1)<fRecEtaWindow
1078       &&TMath::Abs(etaRec)<fRecEtaWindow){
1079
1080       Float_t deltaPhi = phiRec1 - phiRec;
1081
1082       if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1083       if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1084       deltaPhi = TMath::Abs(deltaPhi);
1085       fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
1086       Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
1087       fh2DijetAsymPt->Fill(asym,ptRec1);
1088       fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
1089       fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
1090       fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
1091       Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
1092                          recJets[0].Px()*recJets[1].Px()-
1093                          recJets[0].Py()*recJets[1].Py()-
1094                          recJets[0].Pz()*recJets[1].Py());
1095       minv = TMath::Sqrt(minv);
1096       // with mass == 0;
1097
1098       fh1DijetMinv->Fill(minv);
1099       if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
1100         fh1DijetMinvCut->Fill(minv);
1101         fh2DijetAsymPtCut->Fill(asym,ptRec1);
1102       }
1103     }
1104   }
1105
1106
1107     container[0] = ptRec;
1108     container[1] = etaRec;
1109     container[2] = phiRec;
1110     containerPhiZ[0] = ptRec;
1111     containerPhiZ[1] = phiRec;
1112     if(ptRec>30.&&fDebug>0){
1113       // need to cast to int, otherwise the printf overwrites
1114       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1115       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1116       if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
1117       //  aodH->SetFillAOD(kTRUE);
1118       fAOD->GetHeader()->Print();
1119       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
1120       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1121         AliAODTrack *tr = fAOD->GetTrack(it);
1122         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1123         tr->Print();
1124         //      tr->Dump();
1125         if(fESD){
1126           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1127           esdTr->Print("");
1128           // esdTr->Dump();
1129         }
1130       }
1131     }
1132
1133
1134     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1135     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1136
1137     Float_t zLeading = -1;
1138     if(TMath::Abs(etaRec)<fRecEtaWindow){
1139       fh2JetPtJetPhi->Fill(ptRec,phiRec);
1140       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1141       fh1PtRecIn[ir]->Fill(ptRec,eventW);
1142       // fill the fragmentation function
1143
1144       fh1TmpRho->Reset();
1145
1146       for(int it = 0;it<recParticles.GetEntries();++it){
1147         AliVParticle *part = (AliVParticle*)recParticles.At(it);
1148         Float_t eta = part->Eta();
1149         if(TMath::Abs(eta)<0.9){
1150           Float_t phi = part->Phi();
1151           if(phi<0)phi+=TMath::Pi()*2.;
1152           Float_t dPhi = phi - phiRec;
1153           Float_t dEta = eta - etaRec;
1154           if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1155           if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1156           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1157           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1158         }
1159
1160         Float_t deltaR = recJets[ir].DeltaR(part);
1161         fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1162
1163
1164         if(deltaR<radiusRec){
1165           Float_t z = part->Pt()/ptRec;
1166           if(z>zLeading)zLeading=z;
1167           Float_t lnz =  -1.*TMath::Log(z);
1168           fh2FragRec[ir]->Fill(z,ptRec,eventW);
1169           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1170         }
1171       }
1172       // fill the jet shapes
1173       Float_t rhoSum = 0;
1174       for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1175         Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1176         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1177         rhoSum += rho;
1178         fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1179         fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1180       }
1181     }
1182
1183
1184     containerPhiZ[2] = zLeading;
1185
1186     // Fill Correlation
1187     Int_t ig = iGenIndex[ir];
1188     if(ig>=0 && ig<nGenJets){
1189       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1190       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1191       Double_t ptGen  = genJets[ig].Pt();
1192       Double_t phiGen = genJets[ig].Phi();
1193       if(phiGen<0)phiGen+=TMath::Pi()*2.;
1194       Double_t etaGen = genJets[ig].Eta();
1195
1196       container[3] = ptGen;
1197       container[4] = etaGen;
1198       container[5] = phiGen;
1199       containerPhiZ[3] = ptGen;
1200       //
1201       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1202       //
1203
1204       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1205       if(TMath::Abs(etaRec)<fRecEtaWindow){
1206         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1207         fhnCorrelation->Fill(container);
1208         if(ptGen>0){
1209           Float_t delta = (ptRec-ptGen)/ptGen;
1210           fh2RelPtFGen->Fill(ptGen,delta,eventW);
1211         }
1212         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1213
1214       }// if etarec in window
1215
1216     }
1217     else{
1218       containerPhiZ[3] = 0;
1219       if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1220     }
1221   }// loop over reconstructed jets
1222
1223
1224   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1225   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1226   PostData(1, fHistList);
1227 }
1228
1229
1230 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1231   //
1232   // Create the particle container for the correction framework manager and
1233   // link it
1234   //
1235   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1236   const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1237   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1238   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1239   const Double_t kZmin = 0., kZmax = 1;
1240
1241   // can we neglect migration in eta and phi?
1242   // phi should be no problem since we cover full phi and are phi symmetric
1243   // eta migration is more difficult  due to needed acceptance correction
1244   // in limited eta range
1245
1246   //arrays for the number of bins in each dimension
1247   Int_t iBin[kNvar];
1248   iBin[0] = 320; //bins in pt
1249   iBin[1] =  1; //bins in eta
1250   iBin[2] = 1; // bins in phi
1251
1252   //arrays for lower bounds :
1253   Double_t* binEdges[kNvar];
1254   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1255     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1256
1257   //values for bin lower bounds
1258   //  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);
1259   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1260   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1261   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1262
1263
1264   for(int i = 0;i < kMaxStep*2;++i){
1265     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1266     for (int k=0; k<kNvar; k++) {
1267       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1268     }
1269   }
1270   //create correlation matrix for unfolding
1271   Int_t thnDim[2*kNvar];
1272   for (int k=0; k<kNvar; k++) {
1273     //first half  : reconstructed
1274     //second half : MC
1275     thnDim[k]      = iBin[k];
1276     thnDim[k+kNvar] = iBin[k];
1277   }
1278
1279   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1280   for (int k=0; k<kNvar; k++) {
1281     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1282     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1283   }
1284   fhnCorrelation->Sumw2();
1285
1286   // for second correlation histogram
1287
1288
1289   const Int_t kNvarPhiZ   = 4;
1290   //arrays for the number of bins in each dimension
1291   Int_t iBinPhiZ[kNvarPhiZ];
1292   iBinPhiZ[0] = 80; //bins in pt
1293   iBinPhiZ[1] = 72; //bins in phi
1294   iBinPhiZ[2] = 20; // bins in Z
1295   iBinPhiZ[3] = 80; //bins in ptgen
1296
1297   //arrays for lower bounds :
1298   Double_t* binEdgesPhiZ[kNvarPhiZ];
1299   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1300     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1301
1302   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1303   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1304   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1305   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1306
1307   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1308   for (int k=0; k<kNvarPhiZ; k++) {
1309     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1310   }
1311   fhnCorrelationPhiZRec->Sumw2();
1312
1313
1314   // Add a histogram for Fake jets
1315
1316   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1317     delete [] binEdges[ivar];
1318
1319   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1320     delete [] binEdgesPhiZ[ivar];
1321
1322 }
1323
1324 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1325 {
1326 // Terminate analysis
1327 //
1328     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1329 }
1330
1331
1332 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1333
1334   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1335
1336   Int_t iCount = 0;
1337   if(type==kTrackAOD){
1338     AliAODEvent *aod = 0;
1339     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1340     else aod = AODEvent();
1341     if(!aod){
1342       return iCount;
1343     }
1344     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1345       AliAODTrack *tr = aod->GetTrack(it);
1346       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1347       if(TMath::Abs(tr->Eta())>0.9)continue;
1348       if(fDebug>0){
1349         if(tr->Pt()>20){
1350           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1351           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1352           tr->Print();
1353           //    tr->Dump();
1354           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1355           if(fESD){
1356             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1357             esdTr->Print("");
1358             // esdTr->Dump();
1359           }
1360         }
1361       }
1362       list->Add(tr);
1363       iCount++;
1364     }
1365   }
1366   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1367     AliMCEvent* mcEvent = MCEvent();
1368     if(!mcEvent)return iCount;
1369     // we want to have alivpartilces so use get track
1370     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1371       if(!mcEvent->IsPhysicalPrimary(it))continue;
1372       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1373       if(type == kTrackKineAll){
1374         list->Add(part);
1375         iCount++;
1376       }
1377       else if(type == kTrackKineCharged){
1378         if(part->Particle()->GetPDG()->Charge()==0)continue;
1379         list->Add(part);
1380         iCount++;
1381       }
1382     }
1383   }
1384   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1385     AliAODEvent *aod = 0;
1386     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1387     else aod = AODEvent();
1388     if(!aod)return iCount;
1389     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1390     if(!tca)return iCount;
1391     for(int it = 0;it < tca->GetEntriesFast();++it){
1392       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1393       if(!part->IsPhysicalPrimary())continue;
1394       if(type == kTrackAODMCAll){
1395         list->Add(part);
1396         iCount++;
1397       }
1398       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1399         if(part->Charge()==0)continue;
1400         if(kTrackAODMCCharged){
1401           list->Add(part);
1402         }
1403         else {
1404           if(TMath::Abs(part->Eta())>0.9)continue;
1405           list->Add(part);
1406         }
1407         iCount++;
1408       }
1409     }
1410   }// AODMCparticle
1411   list->Sort();
1412   return iCount;
1413
1414 }