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