PWGJE/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 #include "AliCFContainer.h"
  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   fTRP(0),
  97   fDebug(0),
  98   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
  99   fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
 100   fFilterMask(0),
 101   fEventSelectionMask(0),
 102   fNTrigger(0),
 103   fTriggerBit(0x0),
 104   fNAcceptance(0),
 105   fNBinsLeadingTrackPt(10),
 106   fNBinsMult(20),
 107   fAnalysisType(0),
 108   fTrackTypeRec(kTrackUndef),
 109   fTrackTypeGen(kTrackUndef),
 110   fEventClass(0),
 111   fRPMethod(0),
 112   fAvgTrials(1),
 113   fExternalWeight(1),
 114   fJetRecEtaWindow(0.5),
 115   fTrackRecEtaWindow(0.5),
 116   fMinJetPt(0),
 117   fMinTrackPt(0.15),
 118   fDeltaPhiWindow(90./180.*TMath::Pi()),
 119   fAcceptancePhiMin(0x0),
 120   fAcceptancePhiMax(0x0),
 121   fAcceptanceEtaMin(0x0),
 122   fAcceptanceEtaMax(0x0),
 123   fCentrality(100),
 124   fRPAngle(0),
 125   fMultRec(0),
 126   fMultGen(0),
 127   fTriggerName(0x0),
 128   fh1Xsec(0x0),
 129   fh1Trials(0x0),
 130   fh1AvgTrials(0x0),
 131   fh1PtHard(0x0),
 132   fh1PtHardNoW(0x0),
 133   fh1PtHardTrials(0x0),
 134   fh1ZVtx(0x0),
 135   fh1RP(0x0),
 136   fh1Centrality(0x0),
 137   fh1TmpRho(0x0),
 138   fh2MultRec(0x0),
 139   fh2MultGen(0x0),
 140   fh2RPCentrality(0x0),
 141   fh2PtFGen(0x0),
 142   fh2deltaPt1Pt2(0x0),
 143   fh2RelPtFGen(0x0),
 144   fh3RelPtFGenLeadTrkPt(0x0),
 145   fHistList(0x0)
 146 {
 147
 148   for(int ij = 0;ij <kJetTypes;++ij){
 149     fFlagJetType[ij] = 1; // default = on
 150     fh1NJets[ij] = 0;
 151     fh1SumPtTrack[ij] = 0;
 152     fh1PtJetsIn[ij] = 0;
 153     fh1PtJetsInRej[ij] = 0;
 154     fh1PtJetsInBest[ij] = 0;
 155     fh1PtTracksIn[ij] = 0;
 156     fh1PtTracksInLow[ij] = 0;
 157     fh2NJetsPt[ij]  = 0;
 158     fh2NTracksPt[ij]  = 0;
 159     fp2MultRPPhiTrackPt[ij] = 0;
 160     fp2CentRPPhiTrackPt[ij] = 0;
 161     fhnJetPt[ij] = 0;
 162     fhnJetPtBest[ij] = 0;
 163     fhnJetPtRej[ij] = 0;
 164     fhnJetPtQA[ij] = 0;
 165     fhnTrackPt[ij] = 0;
 166     fhnTrackPtQA[ij] = 0;
 167     for(int i = 0;i <= kMaxJets;++i){
 168       fh2LTrackPtJetPt[ij][i] = 0;
 169       fh1PtIn[ij][i] = 0;
 170     }
 171
 172     fh1DijetMinv[ij] = 0;
 173     fh2DijetDeltaPhiPt[ij] = 0;
 174     fh2DijetAsymPt[ij] = 0;
 175     fh2DijetPt2vsPt1[ij] = 0;
 176     fh2DijetDifvsSum[ij] = 0;
 177
 178   }
 179 }
 180
 181 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
 182   AliAnalysisTaskSE(name),
 183   fJetHeaderRec(0x0),
 184   fJetHeaderGen(0x0),
 185   fAODIn(0x0),
 186   fAODOut(0x0),
 187   fAODExtension(0x0),
 188   fhnJetContainer(0x0),
 189   fhnCorrelation(0x0),
 190   fhnEvent(0x0),
 191   f1PtScale(0x0),
 192   fBranchRec("jets"),
 193   fBranchGen(""),
 194   fBranchBkgRec(""),
 195   fBranchBkgGen(""),
 196   fNonStdFile(""),
 197   fRandomizer(0x0),
 198   fUseAODJetInput(kFALSE),
 199   fUseAODTrackInput(kFALSE),
 200   fUseAODMCInput(kFALSE),
 201   fUseGlobalSelection(kFALSE),
 202   fUseExternalWeightOnly(kFALSE),
 203   fLimitGenJetEta(kFALSE),
 204   fDoMatching(kFALSE),
 205   fNMatchJets(5),
 206   fNRPBins(3),
 207   fTRP(0),
 208   fDebug(0),
 209   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 210   fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
 211   fFilterMask(0),
 212   fEventSelectionMask(0),
 213   fNTrigger(0),
 214   fTriggerBit(0x0),
 215   fNAcceptance(0),
 216   fNBinsLeadingTrackPt(10),
 217   fNBinsMult(20),
 218   fAnalysisType(0),
 219   fTrackTypeRec(kTrackUndef),
 220   fTrackTypeGen(kTrackUndef),
 221   fEventClass(0),
 222   fRPMethod(0),
 223   fAvgTrials(1),
 224   fExternalWeight(1),
 225   fJetRecEtaWindow(0.5),
 226   fTrackRecEtaWindow(0.5),
 227   fMinJetPt(0),
 228   fMinTrackPt(0.15),
 229   fDeltaPhiWindow(90./180.*TMath::Pi()),
 230   fAcceptancePhiMin(0x0),
 231   fAcceptancePhiMax(0x0),
 232   fAcceptanceEtaMin(0x0),
 233   fAcceptanceEtaMax(0x0),
 234   fCentrality(100),
 235   fRPAngle(0),
 236   fMultRec(0),
 237   fMultGen(0),
 238   fTriggerName(0x0),
 239   fh1Xsec(0x0),
 240   fh1Trials(0x0),
 241   fh1AvgTrials(0x0),
 242   fh1PtHard(0x0),
 243   fh1PtHardNoW(0x0),
 244   fh1PtHardTrials(0x0),
 245   fh1ZVtx(0x0),
 246   fh1RP(0x0),
 247   fh1Centrality(0x0),
 248   fh1TmpRho(0x0),
 249   fh2MultRec(0x0),
 250   fh2MultGen(0x0),
 251   fh2RPCentrality(0x0),
 252   fh2PtFGen(0x0),
 253   fh2deltaPt1Pt2(0x0),
 254   fh2RelPtFGen(0x0),
 255   fh3RelPtFGenLeadTrkPt(0x0),
 256   fHistList(0x0)
 257 {
 258
 259   for(int ij = 0;ij <kJetTypes;++ij){
 260     fFlagJetType[ij] = 1; // default = on
 261     fh1NJets[ij] = 0;
 262     fh1SumPtTrack[ij] = 0;
 263     fh1PtJetsIn[ij] = 0;
 264     fh1PtJetsInRej[ij] = 0;
 265     fh1PtJetsInBest[ij] = 0;
 266     fh1PtTracksIn[ij] = 0;
 267     fh1PtTracksInLow[ij] = 0;
 268     fh2NJetsPt[ij]  = 0;
 269     fh2NTracksPt[ij]  = 0;
 270     fp2MultRPPhiTrackPt[ij] = 0;
 271     fp2CentRPPhiTrackPt[ij] = 0;
 272     fhnJetPt[ij] = 0;
 273     fhnJetPtBest[ij] = 0;
 274     fhnJetPtRej[ij] = 0;
 275     fhnJetPtQA[ij] = 0;
 276     fhnTrackPt[ij] = 0;
 277     fhnTrackPtQA[ij] = 0;
 278     for(int i = 0;i <= kMaxJets;++i){
 279       fh2LTrackPtJetPt[ij][i] = 0;
 280       fh1PtIn[ij][i] = 0;
 281     }
 282
 283     fh1DijetMinv[ij] = 0;
 284     fh2DijetDeltaPhiPt[ij] = 0;
 285     fh2DijetAsymPt[ij] = 0;
 286     fh2DijetPt2vsPt1[ij] = 0;
 287     fh2DijetDifvsSum[ij] = 0;
 288   }
 289
 290   DefineOutput(1, TList::Class());
 291 }
 292
 293
 294
 295 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
 296 {
 297
 298
 299
 300   //
 301   // Implemented Notify() to read the cross sections
 302   // and number of trials from pyxsec.root
 303   //
 304
 305   // Fetch the aod also from the input in,
 306   // have todo it in notify
 307
 308
 309   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 310   //  assume that the AOD is in the general output...
 311   fAODOut  = AODEvent();
 312
 313   if(fNonStdFile.Length()!=0){
 314     // case that we have an AOD extension we need can fetch the jets from the extended output
 315     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 316     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 317     if(!fAODExtension){
 318       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 319     }
 320   }
 321
 322
 323   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 324   Float_t xsection = 0;
 325   Float_t ftrials  = 1;
 326
 327   fAvgTrials = 1;
 328   if(tree){
 329     TFile *curfile = tree->GetCurrentFile();
 330     if (!curfile) {
 331       Error("Notify","No current file");
 332       return kFALSE;
 333     }
 334     if(!fh1Xsec||!fh1Trials){
 335       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 336       return kFALSE;
 337     }
 338     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 339     fh1Xsec->Fill("<#sigma>",xsection);
 340     // construct a poor man average trials
 341     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 342     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 343     fh1Trials->Fill("#sum{ntrials}",ftrials);
 344   }
 345
 346   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 347
 348   return kTRUE;
 349 }
 350
 351 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
 352 {
 353
 354
 355   // Connect the AOD
 356
 357   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
 358   OpenFile(1);
 359   if(!fHistList)fHistList = new TList();
 360   PostData(1, fHistList); // post data in any case once
 361
 362   if(!fRandomizer)fRandomizer = new TRandom3(0);
 363
 364   fHistList->SetOwner(kTRUE);
 365   Bool_t oldStatus = TH1::AddDirectoryStatus();
 366   TH1::AddDirectory(kFALSE);
 367
 368
 369
 370   // event npsparse cent, mult
 371   const Int_t nBinsSparse0 = 4;
 372   const Int_t nBins0[nBinsSparse0] = {     100, 500,fNTrigger,125};
 373   const Double_t xmin0[nBinsSparse0]  = {    0,   0, -0.5,-2};
 374   const Double_t xmax0[nBinsSparse0]  = {  100,5000,fNTrigger-0.5,248};
 375
 376
 377   fhnEvent = new THnSparseF("fhnEvent",";cent;mult:trigger;#rho",nBinsSparse0,
 378                             nBins0,xmin0,xmax0);
 379   fHistList->Add(fhnEvent);
 380
 381   if(fDoMatching){
 382     MakeJetContainer();
 383     fHistList->Add(fhnCorrelation);
 384     fHistList->Add(fhnJetContainer);
 385   }
 386   //
 387   //  Histogram
 388
 389
 390
 391   const Int_t nBinPt = 150;
 392   Double_t binLimitsPt[nBinPt+1];
 393   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 394     if(iPt == 0){
 395       binLimitsPt[iPt] = -50.0;
 396     }
 397     else {// 1.0
 398       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
 399     }
 400   }
 401   const Int_t nBinPhi = 90;
 402   Double_t binLimitsPhi[nBinPhi+1];
 403   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 404     if(iPhi==0){
 405       binLimitsPhi[iPhi] = 0;
 406     }
 407     else{
 408       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 409     }
 410   }
 411
 412   const Int_t nBinEta = 40;
 413   Double_t binLimitsEta[nBinEta+1];
 414   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 415     if(iEta==0){
 416       binLimitsEta[iEta] = -2.0;
 417     }
 418     else{
 419       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 420     }
 421   }
 422
 423
 424   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 425   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 426   fHistList->Add(fh1Xsec);
 427   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 428   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 429   fHistList->Add(fh1Trials);
 430   fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
 431   fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
 432   fHistList->Add(fh1AvgTrials);
 433   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 434   fHistList->Add(fh1PtHard);
 435   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 436   fHistList->Add(fh1PtHardNoW);
 437   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 438   fHistList->Add(fh1PtHardTrials);
 439
 440   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
 441   fHistList->Add(fh1ZVtx);
 442
 443
 444   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
 445   fHistList->Add(fh1RP);
 446
 447   fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
 448   fHistList->Add(fh1Centrality);
 449
 450   fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
 451   fHistList->Add(fh2MultRec);
 452   fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
 453   fHistList->Add(fh2MultGen);
 454
 455
 456   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
 457   fHistList->Add(fh2RPCentrality);
 458
 459   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 460   fHistList->Add(fh2PtFGen);
 461
 462   //  fh2deltaPt1Pt2(0x0),
 463   fh2deltaPt1Pt2 = new TH3F("fh2deltaPt1Pt2",Form("%s vs. %s;delta;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 464   fHistList->Add(fh2deltaPt1Pt2);
 465
 466   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
 467   fHistList->Add(fh2RelPtFGen);
 468
 469   const Int_t nBinsLeadingTrackPt = 10;
 470   const Double_t binArrayLeadingTrackPt[nBinsLeadingTrackPt+1] = {0.,1.,2.,3.,4.,5.,6.,8.,10.,12.,200.}; //store pT of leading track in jet
 471
 472   const Int_t nBinDeltaPt = 241;
 473   Double_t binLimitsDeltaPt[nBinDeltaPt+1];
 474   for(Int_t iDeltaPt = 0;iDeltaPt <= nBinDeltaPt;iDeltaPt++){
 475     if(iDeltaPt == 0){
 476       binLimitsDeltaPt[iDeltaPt] = -2.41;
 477     }
 478     else {
 479       binLimitsDeltaPt[iDeltaPt] =  binLimitsDeltaPt[iDeltaPt-1] + 0.02;
 480     }
 481   }
 482
 483   fh3RelPtFGenLeadTrkPt = new TH3F("fh3RelPtFGenLeadTrkPt",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen};p_{T}^{leading track}",nBinPt,binLimitsPt,nBinDeltaPt,binLimitsDeltaPt,nBinsLeadingTrackPt,binArrayLeadingTrackPt);
 484   fHistList->Add(fh3RelPtFGenLeadTrkPt);
 485
 486   for(int ij = 0;ij <kJetTypes;++ij){
 487     TString cAdd = "";
 488     TString cJetBranch = "";
 489     if(ij==kJetRec){
 490       cAdd = "Rec";
 491       cJetBranch = fBranchRec.Data();
 492     }
 493     else if (ij==kJetGen){
 494       cAdd = "Gen";
 495       cJetBranch = fBranchGen.Data();
 496     }
 497     else if (ij==kJetRecFull){
 498       cAdd = "RecFull";
 499       cJetBranch = fBranchRec.Data();
 500     }
 501     else if (ij==kJetGenFull){
 502       cAdd = "GenFull";
 503       cJetBranch = fBranchGen.Data();
 504     }
 505
 506     if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
 507     if(!fFlagJetType[ij])continue;
 508
 509     fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
 510     fHistList->Add(fh1NJets[ij]);
 511
 512     fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 513     fHistList->Add(fh1PtJetsIn[ij]);
 514
 515     fh1PtJetsInRej[ij]  = new TH1F(Form("fh1PtJets%sInRej",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 516     fHistList->Add(fh1PtJetsInRej[ij]);
 517
 518     fh1PtJetsInBest[ij]  = new TH1F(Form("fh1PtJets%sInBest",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 519     fHistList->Add(fh1PtJetsInBest[ij]);
 520
 521     fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
 522     fHistList->Add(fh1PtTracksIn[ij]);
 523
 524     fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
 525     fHistList->Add(fh1PtTracksInLow[ij]);
 526
 527     fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),1000,0.,3000.);
 528     fHistList->Add(fh1SumPtTrack[ij]);
 529
 530     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);
 531     fHistList->Add(fh2NJetsPt[ij]);
 532
 533     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);
 534     fHistList->Add(fh2NTracksPt[ij]);
 535
 536
 537     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");
 538     fHistList->Add(fp2MultRPPhiTrackPt[ij]);
 539     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");
 540     fHistList->Add(fp2CentRPPhiTrackPt[ij]);
 541
 542     // Bins:  Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M
 543     const Int_t nBinsSparse1 = 9;
 544     const Int_t nBinsArea = 10;
 545     Int_t nbinr5=120;
 546     Int_t nbinr5min=-50;
 547
 548     if(cJetBranch.Contains("KT05")){
 549      nbinr5=132;
 550      nbinr5min=-80;}
 551
 552     Int_t nBins1[nBinsSparse1] = {     kMaxJets+1,nbinr5, 10,  fNBinsMult,    fNRPBins, nBinsArea,fNTrigger,fNBinsLeadingTrackPt,fNAcceptance+1};
 553     if(cJetBranch.Contains("RandomCone")){
 554       nBins1[1] = 600;
 555       nBins1[5] = 1;
 556     }
 557     const Double_t xmin1[nBinsSparse1]  = {        -0.5,nbinr5min,  0,   0,        -0.5, 0.,         -0.5,  0.,           -0.5,};
 558     const Double_t xmax1[nBinsSparse1]  = {kMaxJets+0.5,250,100,4000,fNRPBins-0.5,1.0,fNTrigger-0.5,200.,fNAcceptance+0.5};
 559
 560     const Double_t binArrayArea[nBinsArea+1] = {xmin1[5],0.07,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,xmax1[5]};
 561
 562
 563     fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leading track p_{T};acceptance bin",nBinsSparse1,nBins1,xmin1,xmax1);
 564     fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
 565     if(fNBinsLeadingTrackPt>1) fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
 566     fHistList->Add(fhnJetPt[ij]);
 567
 568
 569     // Bins:  pTJet, cent, trigger,
 570     const Int_t nBinsSparse1b = 3;
 571     Int_t nBins1b[nBinsSparse1b] = {120, 10,fNTrigger};
 572     const Double_t xmin1b[nBinsSparse1b]  = {-50,  0,-0.5};
 573     const Double_t xmax1b[nBinsSparse1b]  = {250,100,fNTrigger-0.5};
 574
 575     fhnJetPtBest[ij] = new THnSparseF(Form("fhnJetPtBest%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
 576     fHistList->Add(fhnJetPtBest[ij]);
 577
 578     fhnJetPtRej[ij] = new THnSparseF(Form("fhnJetPtRej%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
 579     fHistList->Add(fhnJetPtRej[ij]);
 580
 581     // Bins:  Jet number: pTJet, cent, eta, phi, Area, trigger, acceptance, signed pT leading
 582     const Int_t nBinsSparse2 = 9;
 583     Int_t nBins2[nBinsSparse2] = {     kMaxJets+1, 60,   8,  18, 18, 100,fNTrigger,fNAcceptance+1,20};
 584     if(cJetBranch.Contains("RandomCone")){
 585       nBins2[5] = 1;
 586     }
 587     const Double_t xmin2[nBinsSparse2]  = {        -0.5, -50,  0,-0.9,              0,  0.,        -0.5,            -0.5, -100};
 588     const Double_t xmax2[nBinsSparse2]  = {kMaxJets+0.5, 250, 80, 0.9, 2.*TMath::Pi(),1.0,fNTrigger-0.5,fNAcceptance+0.5,  100};
 589     fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading",nBinsSparse2,nBins2,xmin2,xmax2);
 590     //    fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
 591     fHistList->Add(fhnJetPtQA[ij]);
 592
 593     // Bins:track number  pTtrack, cent, mult, RP.  charge total bins = 224 k
 594     const Int_t nBinsSparse3 = 7;
 595     const Int_t nRPBinsSparse3 = 3;
 596     const Int_t nBins3[nBinsSparse3] = {       2,    100,     10,   1,    nRPBinsSparse3,fNTrigger,2};
 597     const Double_t xmin3[nBinsSparse3]  = { -0.5,      0,   0,      0,        -0.5,-0.5,-1.5};
 598     const Double_t xmax3[nBinsSparse3]  = { 1.5,     200, 100,   4000,nRPBinsSparse3-0.5,fNTrigger-0.5,1.5};
 599
 600     // change the binning of the pT axis:
 601     Double_t *xPt3 = new Double_t[nBins3[1]+1];
 602     xPt3[0] = 0.;
 603     for(int i = 1; i<=nBins3[1];i++){
 604       if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
 605       else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
 606       else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
 607       else if(xPt3[i-1]<15)xPt3[i] = xPt3[i-1] +  1.; // 62 - 67
 608       else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 2.; // 67 - 72
 609       else if(xPt3[i-1]<60)xPt3[i] = xPt3[i-1] + 5; // 72 - 76
 610       else xPt3[i] = xPt3[i-1] + 10; // 76 - 100 = 140
 611     }
 612
 613     fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP;trigger",nBinsSparse3,nBins3,xmin3,xmax3);
 614     fhnTrackPt[ij]->SetBinEdges(1,xPt3);
 615     fHistList->Add(fhnTrackPt[ij]);
 616     delete [] xPt3;
 617
 618     // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M
 619     const Int_t nBinsSparse4 = 6;
 620     const Int_t nBins4[nBinsSparse4] =    {    2, 50,  10,   20, 180,2};
 621     const Double_t xmin4[nBinsSparse4]  = { -0.5,  0,   0, -1.0,  0.,-1.5};
 622     const Double_t xmax4[nBinsSparse4]  = {  1.5,150, 100,  1.0,2.*TMath::Pi(),1.5};
 623
 624     // change the binning ot the pT axis:
 625     Double_t *xPt4 = new Double_t[nBins4[1]+1];
 626     xPt4[0] = 0.;
 627     for(int i = 1; i<=nBins4[1];i++){
 628       if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1;
 629       else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5;
 630       else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] +  1.;
 631       else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] +  2.5;
 632       else xPt4[i] = xPt4[i-1] + 5.;
 633     }
 634     fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi;sign",nBinsSparse4,nBins4,xmin4,xmax4);
 635     fhnTrackPtQA[ij]->SetBinEdges(1,xPt4);
 636     fHistList->Add(fhnTrackPtQA[ij]);
 637     delete [] xPt4;
 638
 639     for(int i = 0;i <= kMaxJets;++i){
 640       fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
 641       fHistList->Add(fh1PtIn[ij][i]);
 642
 643
 644       if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 645       fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
 646                                          Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
 647                                               cAdd.Data()),
 648                                          200,0.,200.,nBinPt,binLimitsPt);
 649       fHistList->Add(fh2LTrackPtJetPt[ij][i]);
 650     }
 651
 652
 653     fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 654     fHistList->Add(fh1DijetMinv[ij]);
 655
 656     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);
 657     fHistList->Add(fh2DijetDeltaPhiPt[ij]);
 658
 659     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);
 660     fHistList->Add(fh2DijetAsymPt[ij]);
 661
 662     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.);
 663     fHistList->Add(fh2DijetPt2vsPt1[ij]);
 664     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.);
 665     fHistList->Add( fh2DijetDifvsSum[ij]);
 666   }
 667   // =========== Switch on Sumw2 for all histos ===========
 668   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 669     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 670     if (h1){
 671       h1->Sumw2();
 672       continue;
 673     }
 674     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 675     if(hn)hn->Sumw2();
 676   }
 677   TH1::AddDirectory(oldStatus);
 678 }
 679
 680 void AliAnalysisTaskJetSpectrum2::Init()
 681 {
 682   //
 683   // Initialization
 684   //
 685
 686   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 687
 688 }
 689
 690 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
 691
 692
 693   Bool_t selected = kTRUE;
 694
 695   if(fUseGlobalSelection&&fEventSelectionMask==0){
 696     selected = AliAnalysisHelperJetTasks::Selected();
 697   }
 698   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 699     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 700   }
 701
 702   if(fEventClass>0){
 703     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 704   }
 705
 706   if(!selected){
 707     // no selection by the service task, we continue
 708     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 709     PostData(1, fHistList);
 710     return;
 711   }
 712
 713
 714   static AliAODEvent* aod = 0;
 715
 716   // take all other information from the aod we take the tracks from
 717   if(!aod){
 718     if(fUseAODTrackInput)aod = fAODIn;
 719     else aod = fAODOut;
 720   }
 721
 722
 723   //
 724   // Execute analysis for current event
 725   //
 726   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 727   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 728   if(!aodH){
 729     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 730     return;
 731   }
 732
 733   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 734   TClonesArray *aodRecJets = 0;
 735
 736   if(fAODOut&&!aodRecJets){
 737     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
 738   }
 739   if(fAODExtension&&!aodRecJets){
 740     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
 741   }
 742   if(fAODIn&&!aodRecJets){
 743     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
 744   }
 745
 746
 747
 748   if(!aodRecJets){
 749     if(fDebug){
 750
 751       Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 752       if(fAODIn){
 753         Printf("Input AOD >>>>");
 754         fAODIn->Print();
 755       }
 756       if(fAODExtension){
 757         Printf("AOD Extension >>>>");
 758         fAODExtension->Print();
 759       }
 760       if(fAODOut){
 761         Printf("Output AOD >>>>");
 762         fAODOut->Print();
 763       }
 764     }
 765     return;
 766   }
 767
 768   TClonesArray *aodGenJets = 0;
 769   if(fBranchGen.Length()>0){
 770     if(fAODOut&&!aodGenJets){
 771       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
 772     }
 773     if(fAODExtension&&!aodGenJets){
 774       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
 775     }
 776     if(fAODIn&&!aodGenJets){
 777       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
 778     }
 779
 780     if(!aodGenJets){
 781       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
 782       return;
 783     }
 784   }
 785
 786   TClonesArray *aodBackRecJets = 0;
 787   if(fBranchBkgRec.Length()>0){
 788     if(fAODOut&&!aodBackRecJets){
 789       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
 790     }
 791     if(fAODExtension&&!aodBackRecJets){
 792       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
 793     }
 794     if(fAODIn&&!aodBackRecJets){
 795       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
 796     }
 797
 798     if(!aodBackRecJets){
 799       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
 800       return;
 801     }
 802   }
 803
 804
 805   TClonesArray *aodBackGenJets = 0;
 806
 807   if(fBranchBkgGen.Length()>0){
 808     if(fAODOut&&!aodBackGenJets){
 809       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
 810     }
 811     if(fAODExtension&&!aodBackGenJets){
 812       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
 813     }
 814     if(fAODIn&&!aodBackGenJets){
 815       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
 816     }
 817
 818     if(!aodBackGenJets){
 819       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
 820       return;
 821     }
 822   }
 823
 824
 825   // new Scheme
 826   // first fill all the pure  histograms, i.e. full jets
 827   // in case of the correaltion limit to the n-leading jets
 828
 829   // reconstructed
 830
 831
 832   // generated
 833
 834
 835   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
 836
 837
 838
 839   TList genJetsList;         // full acceptance
 840   TList genJetsListCut;      // acceptance cut
 841   TList recJetsList;         // full acceptance
 842   TList recJetsListCut;      // acceptance cut
 843
 844
 845   GetListOfJets(&recJetsList,aodRecJets,0);
 846   GetListOfJets(&recJetsListCut,aodRecJets,1);
 847
 848   if(fBranchGen.Length()>0){
 849     GetListOfJets(&genJetsList,aodGenJets,0);
 850     GetListOfJets(&genJetsListCut,aodGenJets,1);
 851   }
 852
 853   Double_t eventW = 1;
 854   Double_t ptHard = 0;
 855   Double_t nTrials = 1; // Trials for MC trigger
 856   fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
 857
 858   // Getting some global properties
 859   fCentrality = GetCentrality();
 860   if(fCentrality<=0)fCentrality = 0;
 861   fh1Centrality->Fill(fCentrality);
 862
 863
 864   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 865     // this is the part we only use when we have MC information
 866     AliMCEvent* mcEvent = MCEvent();
 867     //    AliStack *pStack = 0;
 868     if(!mcEvent){
 869       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 870       return;
 871     }
 872     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 873     if(pythiaGenHeader){
 874       nTrials = pythiaGenHeader->Trials();
 875       ptHard  = pythiaGenHeader->GetPtHard();
 876       int iProcessType = pythiaGenHeader->ProcessType();
 877       // 11 f+f -> f+f
 878       // 12 f+barf -> f+barf
 879       // 13 f+barf -> g+g
 880       // 28 f+g -> f+g
 881       // 53 g+g -> f+barf
 882       // 68 g+g -> g+g
 883       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 884       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 885     }
 886   }// (fAnalysisType&kMCESD)==kMCESD)
 887
 888
 889   // we simply fetch the tracks/mc particles as a list of AliVParticles
 890
 891   TList recParticles;
 892   TList genParticles;
 893
 894   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 895
 896   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 897   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 898   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 899
 900   //  CalculateReactionPlaneAngle(&recParticles);
 901   fRPAngle = 0;
 902
 903   if(fRPMethod==0)fRPAngle = aod->GetHeader()->GetEventplane();
 904   else if(fRPMethod==1||fRPMethod==2){
 905     fRPAngle = aod->GetHeader()->GetQTheta(fRPMethod);
 906   }
 907   fh1RP->Fill(fRPAngle);
 908   fh2RPCentrality->Fill(fCentrality,fRPAngle);
 909   // Event control and counting ...
 910   // MC
 911   fh1PtHard->Fill(ptHard,eventW);
 912   fh1PtHardNoW->Fill(ptHard,1);
 913   fh1PtHardTrials->Fill(ptHard,nTrials);
 914
 915   // Real
 916   if(aod->GetPrimaryVertex()){// No vtx for pure MC
 917     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
 918   }
 919
 920
 921   Int_t recMult1 = recParticles.GetEntries();
 922   Int_t genMult1 = genParticles.GetEntries();
 923
 924   Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
 925   Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
 926
 927   fh2MultRec->Fill(recMult1,recMult2);
 928   fh2MultGen->Fill(genMult1,genMult2);
 929   fMultRec = recMult1;
 930   if(fMultRec<=0)fMultRec = recMult2;
 931   fMultGen = genMult1;
 932   if(fMultGen<=0)fMultGen = genMult2;
 933
 934   Double_t var0[4] = {0,};
 935   var0[0] = fCentrality;
 936   var0[1] = fMultRec;
 937   for(int it=0;it<fNTrigger;it++){
 938     if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 939       var0[2] = it;
 940       var0[3] = GetRho(recJetsList);
 941       fhnEvent->Fill(var0);
 942     }
 943   }
 944   // the loops for rec and gen should be indentical... pass it to a separate
 945   // function ...
 946   // Jet Loop
 947   // Track Loop
 948   // Jet Jet Loop
 949   // Jet Track loop
 950
 951   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 952   FillJetHistos(recJetsList,recParticles,kJetRecFull);
 953   FillTrackHistos(recParticles,kJetRec);
 954
 955   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 956   FillJetHistos(genJetsList,genParticles,kJetGenFull);
 957   FillTrackHistos(genParticles,kJetGen);
 958
 959   // Here follows the jet matching part
 960   // delegate to a separated method?
 961
 962   if(fDoMatching){
 963     FillMatchHistos(recJetsList,genJetsList);
 964   }
 965
 966   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 967   PostData(1, fHistList);
 968 }
 969
 970 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 971
 972   if(iType>=kJetTypes){
 973     return;
 974   }
 975   if(!fFlagJetType[iType])return;
 976
 977   Int_t refMult = fMultRec;
 978   if(iType==kJetGen||iType==kJetGenFull){
 979     refMult = fMultGen;
 980   }
 981
 982   Int_t nJets = jetsList.GetEntries();
 983   fh1NJets[iType]->Fill(nJets);
 984
 985   if(nJets<=0)return;
 986
 987   Float_t ptOld = 1.E+32;
 988   Float_t pT0 = 0;
 989   Float_t pT1 = 0;
 990   Float_t phi0 = 0;
 991   Float_t phi1 = 0;
 992   Int_t ij0 = -1;
 993   Int_t ij1 = -1;
 994
 995   Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
 996   Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
 997
 998   var1[2] = fCentrality;
 999   var1b[1] = fCentrality;
1000   var1[3] = refMult;
1001
1002
1003
1004   Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
1005   var2[2] = fCentrality;
1006
1007   for(int ij = 0;ij < nJets;ij++){
1008     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
1009     Float_t ptJet = jet->Pt();
1010
1011
1012     if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
1013
1014     var1b[0] = ptJet;
1015     if(jet->Trigger()&fJetTriggerBestMask){
1016       fh1PtJetsInBest[iType]->Fill(ptJet);
1017       for(int it = 0;it <fNTrigger;it++){
1018         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1019           var1b[2] = it;
1020           fhnJetPtBest[iType]->Fill(var1b);
1021         }
1022       }
1023     }
1024     if(jet->Trigger()&fJetTriggerExcludeMask){
1025       fh1PtJetsInRej[iType]->Fill(ptJet);
1026       for(int it = 0;it <fNTrigger;it++){
1027         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1028           var1b[2] = it;
1029           fhnJetPtRej[iType]->Fill(var1b);
1030         }
1031       }
1032       continue;
1033     }
1034     fh1PtJetsIn[iType]->Fill(ptJet);
1035     ptOld = ptJet;
1036
1037     // find the dijets assume sorting and acceptance cut...
1038     if(ij==0){
1039       ij0 = ij;
1040       pT0 = ptJet;
1041       phi0 = jet->Phi();
1042       if(phi0<0)phi0+=TMath::Pi()*2.;
1043     }
1044     else if(ptJet>pT1){
1045       // take only the backward hemisphere??
1046       phi1 = jet->Phi();
1047       if(phi1<0)phi1+=TMath::Pi()*2.;
1048       Float_t dPhi = phi1 - phi0;
1049       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1050       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1051       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1052         ij1 = ij;
1053         pT1 = ptJet;
1054       }
1055     }
1056     // fill jet histos for kmax jets
1057     Float_t phiJet = jet->Phi();
1058     Float_t etaJet = jet->Eta();
1059     if(phiJet<0)phiJet+=TMath::Pi()*2.;
1060     fh1TmpRho->Reset();
1061     if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1062
1063     fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1064     // fill leading jets...
1065     AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1066     //AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1067     Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1068     Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1069     AliVParticle *tt = (AliVParticle*)particlesList.At(0);
1070     Float_t ttphi=tt->Phi();
1071     if(ttphi<0)ttphi+=TMath::Pi()*2.;
1072     Float_t ttpt=tt->Pt();
1073     Int_t phiBintt = GetPhiBin(ttphi-fRPAngle);
1074      Double_t dphitrigjet=RelativePhi(ttphi,phiJet);
1075     if(fTRP==1){
1076     if(TMath::Abs(dphitrigjet)<TMath::Pi()-0.6) continue;
1077     var1[1] = ptJet;
1078     var1[4] = phiBintt;
1079     var1[5] = jet->EffectiveAreaCharged();
1080     var1[7] = ttpt;
1081     var1[8] = CheckAcceptance(phiJet,etaJet);}
1082
1083      if(fTRP==0){
1084     var1[1] = ptJet;
1085     var1[4] = phiBin;
1086     var1[5] = jet->EffectiveAreaCharged();
1087     var1[7] = ptLead;
1088     var1[8] = CheckAcceptance(phiJet,etaJet);}
1089
1090
1091     var2[1] = ptJet;
1092     var2[3] = etaJet;
1093     var2[4] = phiJet;
1094     var2[5] = jet->EffectiveAreaCharged();
1095     var2[7] = var1[8];
1096     var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1097
1098     if(ij<kMaxJets){
1099       fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1100       var1[0] = ij;
1101       var2[0] = ij;
1102       for(int it = 0;it <fNTrigger;it++){
1103         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1104           var1[6] = it;
1105           var2[6] = it;
1106           fhnJetPt[iType]->Fill(var1);
1107           fhnJetPtQA[iType]->Fill(var2);
1108         }
1109       }
1110     }
1111     var1[0] = kMaxJets;// fill for all jets
1112     var2[0] = kMaxJets;// fill for all jets
1113     for(int it = 0;it <fNTrigger;it++){
1114       if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1115         var1[6] = it;
1116         var2[6] = it;
1117         fhnJetPt[iType]->Fill(var1);
1118         fhnJetPtQA[iType]->Fill(var2);
1119       }
1120     }
1121
1122     fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
1123
1124     if(particlesList.GetSize()&&ij<kMaxJets){
1125       // Particles... correlated with jets...
1126       for(int it = 0;it<particlesList.GetEntries();++it){
1127         AliVParticle *part = (AliVParticle*)particlesList.At(it);
1128         Float_t deltaR = jet->DeltaR(part);
1129         if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1130       }
1131       // fill the jet shapes
1132     }// if we have particles
1133   }// Jet Loop
1134
1135
1136   // do something with dijets...
1137   if(ij0>=0&&ij1>0){
1138     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1139     Double_t ptJet0 = jet0->Pt();
1140     Double_t phiJet0 = jet0->Phi();
1141     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1142
1143     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1144     Double_t ptJet1 = jet1->Pt();
1145     Double_t phiJet1 = jet1->Phi();
1146     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1147
1148     Float_t deltaPhi = phiJet0 - phiJet1;
1149     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1150     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1151     deltaPhi = TMath::Abs(deltaPhi);
1152     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1153
1154     Float_t asym = 9999;
1155     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1156     fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1157     fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1158     fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1159     Float_t minv = 2.*(jet0->P()*jet1->P()-
1160                        jet0->Px()*jet1->Px()-
1161                        jet0->Py()*jet1->Py()-
1162                        jet0->Pz()*jet1->Pz());    // assume mass == 0;
1163     if(minv<0)minv=0; // prevent numerical instabilities
1164     minv = TMath::Sqrt(minv);
1165     fh1DijetMinv[iType]->Fill(minv);
1166   }
1167
1168
1169
1170   // count down the jets above thrueshold
1171   Int_t nOver = nJets;
1172   if(nOver>0){
1173     TIterator *jetIter = jetsList.MakeIterator();
1174     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1175     if(tmpJet){
1176       Float_t pt = tmpJet->Pt();
1177       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1178         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1179         while(pt<ptCut&&tmpJet){
1180           nOver--;
1181           tmpJet = (AliAODJet*)(jetIter->Next());
1182           if(tmpJet){
1183             pt = tmpJet->Pt();
1184           }
1185         }
1186         if(nOver<=0)break;
1187         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1188       }
1189     }
1190     delete jetIter;
1191   }
1192 }
1193
1194 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1195
1196   if(fFlagJetType[iType]<=0)return;
1197   Int_t refMult = fMultRec;
1198   if(iType==kJetGen||iType==kJetGenFull){
1199     refMult = fMultGen;
1200
1201   }
1202
1203   //
1204   Double_t var3[7]; // track number;p_{T};cent;#tracks;RP;cahrge
1205   var3[2] = fCentrality;
1206   var3[3] = refMult;
1207   Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
1208   var4[2] = fCentrality;
1209   Int_t nTrackOver = particlesList.GetSize();
1210   // do the same for tracks and jets
1211   if(nTrackOver>0){
1212     TIterator *trackIter = particlesList.MakeIterator();
1213     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1214     Float_t pt = tmpTrack->Pt();
1215     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1216       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1217       while(pt<ptCut&&tmpTrack){
1218         nTrackOver--;
1219         tmpTrack = (AliVParticle*)(trackIter->Next());
1220         if(tmpTrack){
1221           pt = tmpTrack->Pt();
1222         }
1223       }
1224       if(nTrackOver<=0)break;
1225       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1226     }
1227
1228
1229     trackIter->Reset();
1230     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1231     Float_t sumPt = 0;
1232
1233     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1234       Float_t tmpPt = tmpTrack->Pt();
1235       fh1PtTracksIn[iType]->Fill(tmpPt);
1236       fh1PtTracksInLow[iType]->Fill(tmpPt);
1237
1238       sumPt += tmpPt;
1239       Float_t tmpPhi = tmpTrack->Phi();
1240       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1241
1242
1243       Float_t phiRP = tmpPhi-fRPAngle;
1244       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1245       if(phiRP<0)phiRP += TMath::Pi();
1246       if(phiRP<0)phiRP += TMath::Pi();
1247       const float allPhi = -1./180.*TMath::Pi();
1248
1249       if(tmpPt<100){
1250         fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1251         fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1252
1253         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1254         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1255       }
1256       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1257       var3[0] = 1;
1258       var3[1] = tmpPt;
1259       var3[4] = phiBin;
1260       var3[6] = tmpTrack->Charge();
1261
1262       var4[0] = 1;
1263       var4[1] = tmpPt;
1264       var4[3] = tmpTrack->Eta();
1265       var4[4] = tmpPhi;
1266       var4[5] = tmpTrack->Charge();
1267
1268
1269       for(int it = 0;it <fNTrigger;it++){
1270         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1271           var3[0] = 1;
1272           var3[5] = it;
1273           fhnTrackPt[iType]->Fill(var3);
1274           if(tmpTrack==leading){
1275             var3[0] = 0;
1276             fhnTrackPt[iType]->Fill(var3);
1277           }
1278         }
1279       }
1280
1281
1282
1283
1284       fhnTrackPtQA[iType]->Fill(var4);
1285       if(tmpTrack==leading){
1286         var4[0] = 0;
1287         fhnTrackPtQA[iType]->Fill(var4);
1288         continue;
1289       }
1290
1291
1292
1293
1294
1295     }
1296     fh1SumPtTrack[iType]->Fill(sumPt);
1297     delete trackIter;
1298   }
1299
1300 }
1301
1302
1303 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1304
1305
1306   // Fill al the matching histos
1307   // It is important that the acceptances for the mathing are as large as possible
1308   // to avoid false matches on the edge of acceptance
1309   // therefore we add some extra matching jets as overhead
1310
1311   static TArrayI aGenIndex(recJetsList.GetEntries());
1312   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1313
1314   static TArrayI aRecIndex(genJetsList.GetEntries());
1315   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1316
1317   if(fDebug){
1318     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1319     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1320   }
1321   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1322                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1323                                             aGenIndex,aRecIndex,fDebug);
1324
1325   if(fDebug){
1326     for(int i = 0;i< aGenIndex.GetSize();++i){
1327       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1328     }
1329     for(int i = 0;i< aRecIndex.GetSize();++i){
1330       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1331     }
1332   }
1333
1334   Double_t container[8];
1335   Double_t containerRec[4];
1336   Double_t containerGen[4];
1337
1338   // loop over generated jets
1339   // consider the
1340   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1341     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1342     Double_t ptGen = genJet->Pt();
1343     Double_t phiGen = genJet->Phi();
1344     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1345     Double_t etaGen = genJet->Eta();
1346     containerGen[0] = ptGen;
1347     containerGen[1] = etaGen;
1348     containerGen[2] = phiGen;
1349     containerGen[3] = genJet->GetPtLeading();
1350
1351     fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1352
1353     if(JetSelected(genJet))
1354       fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1355
1356     Int_t ir = aRecIndex[ig];
1357     if(ir>=0&&ir<recJetsList.GetEntries()){
1358       AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1359
1360       fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1361
1362       if(JetSelected(genJet)){
1363         fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1364
1365         if(JetSelected(recJet)) {
1366           fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1367
1368         }
1369         containerGen[3] = recJet->GetPtLeading();
1370         fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner with leading track on reconstructed level
1371       }
1372     }
1373   }// loop over generated jets used for matching...
1374
1375
1376
1377   // loop over reconstructed jets
1378   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1379     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1380     Double_t etaRec = recJet->Eta();
1381     Double_t ptRec = recJet->Pt();
1382     Double_t phiRec = recJet->Phi();
1383     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1384
1385     containerRec[0] = ptRec;
1386     containerRec[1] = etaRec;
1387     containerRec[2] = phiRec;
1388     containerRec[3] = recJet->GetPtLeading();
1389
1390     container[0] = ptRec;
1391     container[1] = etaRec;
1392     container[2] = phiRec;
1393     container[3] = recJet->GetPtLeading();
1394
1395     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1396
1397     if(JetSelected(recJet)){
1398       fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
1399       // Fill Correlation
1400       Int_t ig = aGenIndex[ir];
1401       if(ig>=0 && ig<genJetsList.GetEntries()){
1402         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1403         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1404         Double_t ptGen  = genJet->Pt();
1405         Double_t phiGen = genJet->Phi();
1406         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1407         Double_t etaGen = genJet->Eta();
1408
1409         containerGen[0] = ptGen;
1410         containerGen[1] = etaGen;
1411         containerGen[2] = phiGen;
1412         containerGen[3] = genJet->GetPtLeading();
1413
1414         container[4] = ptGen;
1415         container[5] = etaGen;
1416         container[6] = phiGen;
1417         container[7] = genJet->GetPtLeading();
1418
1419         fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1420
1421         fhnCorrelation->Fill(container);
1422
1423         Float_t delta = (ptRec-ptGen)/ptGen;
1424         fh2RelPtFGen->Fill(ptGen,delta);
1425         fh3RelPtFGenLeadTrkPt->Fill(ptGen,delta,recJet->GetPtLeading());
1426         fh2PtFGen->Fill(ptGen,ptRec);
1427
1428         fh2deltaPt1Pt2->Fill(ptRec-ptGen,ptRec,ptGen);
1429
1430       }
1431     }// loop over reconstructed jets
1432   }
1433   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1434 }
1435
1436
1437 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1438   //
1439   // Create the particle container for the correction framework manager and
1440   // link it
1441   //
1442   const Int_t kNvar   = 4 ; //number of variables on the grid:pt,eta, phi
1443   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1444   const Double_t kEtamin = -1.0, kEtamax = 1.0;
1445   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1446   const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
1447
1448   // can we neglect migration in eta and phi?
1449   // phi should be no problem since we cover full phi and are phi symmetric
1450   // eta migration is more difficult  due to needed acceptance correction
1451   // in limited eta range
1452
1453   //arrays for the number of bins in each dimension
1454   Int_t iBin[kNvar];
1455   iBin[0] = 125; //bins in pt
1456   iBin[1] =  4; //bins in eta
1457   iBin[2] =  1; // bins in phi
1458   iBin[3] =  10; // bins in leading track Pt
1459
1460   //arrays for lower bounds :
1461   Double_t* binEdges[kNvar];
1462   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1463     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1464
1465   //values for bin lower bounds
1466   //  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);
1467   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1468   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1469   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1470   binEdges[3][0]= kPtLeadingTrackPtMin;
1471   binEdges[3][1]= 1.;
1472   binEdges[3][2]= 2.;
1473   binEdges[3][3]= 3.;
1474   binEdges[3][4]= 4.;
1475   binEdges[3][5]= 5.;
1476   binEdges[3][6]= 6.;
1477   binEdges[3][7]= 8.;
1478   binEdges[3][8]= 10.;
1479   binEdges[3][9]= 12.;
1480   binEdges[3][10]= kPtLeadingTrackPtMax;
1481
1482   fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
1483   for (int k=0; k<kNvar; k++) {
1484     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1485   }
1486
1487   //create correlation matrix for unfolding
1488   Int_t thnDim[2*kNvar];
1489   for (int k=0; k<kNvar; k++) {
1490     //first half  : reconstructed
1491     //second half : MC
1492     thnDim[k]      = iBin[k];
1493     thnDim[k+kNvar] = iBin[k];
1494   }
1495
1496   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
1497   for (int k=0; k<kNvar; k++) {
1498     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1499     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1500   }
1501
1502   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1503     delete [] binEdges[ivar];
1504
1505
1506 }
1507
1508 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1509 {
1510   // Terminate analysis
1511   //
1512   if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1513 }
1514
1515
1516 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1517
1518   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1519
1520   Int_t iCount = 0;
1521   if(type==kTrackAOD){
1522     AliAODEvent *aod = 0;
1523     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1524     else aod = AODEvent();
1525     if(!aod){
1526       return iCount;
1527     }
1528     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1529       AliAODTrack *tr = aod->GetTrack(it);
1530       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1531       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1532       if(tr->Pt()<fMinTrackPt)continue;
1533       if(fDebug>0){
1534         if(tr->Pt()>20){
1535           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1536           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1537           tr->Print();
1538           //    tr->Dump();
1539           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1540           if(fESD){
1541             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1542             esdTr->Print("");
1543             // esdTr->Dump();
1544           }
1545         }
1546       }
1547       list->Add(tr);
1548       iCount++;
1549     }
1550   }
1551   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1552     AliMCEvent* mcEvent = MCEvent();
1553     if(!mcEvent)return iCount;
1554     // we want to have alivpartilces so use get track
1555     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1556       if(!mcEvent->IsPhysicalPrimary(it))continue;
1557       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1558       if(part->Pt()<fMinTrackPt)continue;
1559       if(type == kTrackKineAll){
1560         list->Add(part);
1561         iCount++;
1562       }
1563       else if(type == kTrackKineCharged){
1564         if(part->Particle()->GetPDG()->Charge()==0)continue;
1565         list->Add(part);
1566         iCount++;
1567       }
1568     }
1569   }
1570   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1571     AliAODEvent *aod = 0;
1572     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1573     else aod = AODEvent();
1574     if(!aod)return iCount;
1575     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1576     if(!tca)return iCount;
1577     for(int it = 0;it < tca->GetEntriesFast();++it){
1578       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1579       if(!part)continue;
1580       if(part->Pt()<fMinTrackPt)continue;
1581       if(!part->IsPhysicalPrimary())continue;
1582       if(type == kTrackAODMCAll){
1583         list->Add(part);
1584         iCount++;
1585       }
1586       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1587         if(part->Charge()==0)continue;
1588         if(kTrackAODMCCharged){
1589           list->Add(part);
1590         }
1591         else {
1592           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1593           list->Add(part);
1594         }
1595         iCount++;
1596       }
1597     }
1598   }// AODMCparticle
1599   list->Sort();
1600   return iCount;
1601
1602 }
1603
1604
1605 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1606   AliAODEvent *aod = 0;
1607   if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1608   else aod = AODEvent();
1609   if(!aod){
1610     return 101;
1611   }
1612   return aod->GetHeader()->GetCentrality();
1613 }
1614
1615
1616
1617 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1618   Bool_t selected = false;
1619
1620   if(!jet)return selected;
1621
1622   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1623     selected = kTRUE;
1624   }
1625   return selected;
1626
1627 }
1628
1629 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1630
1631   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1632   Int_t iCount = 0;
1633
1634   if(!jarray){
1635     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1636     return iCount;
1637   }
1638
1639
1640   for(int ij=0;ij<jarray->GetEntries();ij++){
1641     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1642     if(!jet)continue;
1643     if(type==0){
1644       // No cut at all, main purpose here is sorting
1645       list->Add(jet);
1646       iCount++;
1647     }
1648     else if(type == 1){
1649       // eta cut
1650       if(JetSelected(jet)){
1651         list->Add(jet);
1652         iCount++;
1653       }
1654     }
1655   }
1656
1657   list->Sort();
1658   return iCount;
1659
1660 }
1661
1662
1663 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1664   if(!jets)return 0;
1665
1666   Int_t refMult = 0;
1667   for(int ij = 0;ij < jets->GetEntries();++ij){
1668     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1669     if(!jet)continue;
1670     TRefArray *refs = jet->GetRefTracks();
1671     if(!refs)continue;
1672     refMult += refs->GetEntries();
1673   }
1674   return refMult;
1675
1676 }
1677
1678
1679 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1680   if(!jet)return 0;
1681   TRefArray *refs = jet->GetRefTracks();
1682   if(!refs) return 0;
1683   AliVParticle *leading = 0;
1684   Float_t fMaxPt = 0;
1685   for(int i = 0;i<refs->GetEntries();i++){
1686     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1687     if(!tmp)continue;
1688     if(tmp->Pt()>fMaxPt){
1689       leading = tmp;
1690       fMaxPt = tmp->Pt();
1691     }
1692   }
1693   return leading;
1694 }
1695
1696
1697 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1698   if(!jet)return 0;
1699   if(!list) return 0;
1700   AliVParticle *leading = 0;
1701   Float_t fMaxPt = 0;
1702   for(int i = 0;i<list->GetEntries();i++){
1703     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1704     if(!tmp)continue;
1705     if(jet->DeltaR(tmp)>r)continue;
1706     if(tmp->Pt()>fMaxPt){
1707       leading = tmp;
1708       fMaxPt = tmp->Pt();
1709     }
1710   }
1711   return leading;
1712 }
1713
1714
1715 Double_t AliAnalysisTaskJetSpectrum2::RelativePhi(Double_t mphi,Double_t vphi){
1716
1717   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1718   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1719   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1720   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1721   double dphi = mphi-vphi;
1722   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1723   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1724   return dphi;//dphi in [-Pi, Pi]
1725 }
1726
1727
1728
1729
1730
1731
1732 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1733 {
1734   Int_t phibin=-1;
1735   if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1736   Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1737   phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1738   if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1739   return phibin;
1740 }
1741
1742 void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1743   //
1744   // set number of trigger bins
1745   //
1746   if(n>0){
1747     fNTrigger = n;
1748     delete [] fTriggerName;
1749     fTriggerName = new TString [fNTrigger];
1750     delete [] fTriggerBit;fTriggerBit = 0;
1751     fTriggerBit = new UInt_t [fNTrigger];
1752   }
1753   else{
1754     fNTrigger = 0;
1755   }
1756 }
1757
1758
1759 void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1760   //
1761   // set trigger bin
1762   //
1763   if(i<fNTrigger){
1764     fTriggerBit[i] = it;
1765     fTriggerName[i] =  c;
1766   }
1767 }
1768
1769
1770 void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1771   //
1772   // Set number of acceptance bins
1773   //
1774
1775   if(n>0){
1776     fNAcceptance = n;
1777     delete [] fAcceptancePhiMin;
1778     delete [] fAcceptancePhiMax;
1779     delete [] fAcceptanceEtaMin;
1780     delete [] fAcceptanceEtaMax;
1781
1782     fAcceptancePhiMin = new Float_t[fNAcceptance];
1783     fAcceptancePhiMax = new Float_t[fNAcceptance];
1784     fAcceptanceEtaMin = new Float_t[fNAcceptance];
1785     fAcceptanceEtaMax = new Float_t[fNAcceptance];
1786   }
1787   else{
1788     fNTrigger = 0;
1789   }
1790 }
1791
1792 void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1793   //
1794   // Set acceptance bins
1795   //
1796
1797   if(i<fNAcceptance){
1798     Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1799
1800     fAcceptancePhiMin[i] = phiMin;
1801     fAcceptancePhiMax[i] = phiMax;
1802     fAcceptanceEtaMin[i] = etaMin;
1803     fAcceptanceEtaMax[i] = etaMax;
1804   }
1805   else{
1806     fNTrigger = 0;
1807   }
1808 }
1809
1810
1811 Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1812
1813   //
1814   // return aceptnace bin do no use ovelapping bins
1815   //
1816
1817   for(int i = 0;i<fNAcceptance;i++){
1818     if(phi<fAcceptancePhiMin[i])continue;
1819     if(phi>fAcceptancePhiMax[i])continue;
1820     if(eta<fAcceptanceEtaMin[i])continue;
1821     if(eta>fAcceptanceEtaMax[i])continue;
1822     return i;
1823   }
1824   // catch the rest
1825   return fNAcceptance;
1826 }
1827
1828 Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){
1829
1830   // invert the correction
1831   AliAODJet *jet = (AliAODJet*)list.At(0);  // highest pt jet
1832   if(!jet)return -1;
1833   if(jet->EffectiveAreaCharged()<=0)return -1;
1834   Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged();
1835   return rho;
1836
1837 }
1838
1839
1840
1841 AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1842   //
1843   delete [] fTriggerBit;
1844   delete [] fTriggerName;
1845
1846   delete [] fAcceptancePhiMin;
1847   delete [] fAcceptancePhiMax;
1848   delete [] fAcceptanceEtaMin;
1849   delete [] fAcceptanceEtaMax;
1850
1851
1852
1853 }