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