PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx

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