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