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