PWGJE/AliAnalysisTaskJetSpectrum2.cxx

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