PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx

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