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