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       Error("%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 = 3;
 364   const Int_t nBins0[nBinsSparse0] = {     100, 500,fNTrigger};
 365   const Double_t xmin0[nBinsSparse0]  = {    0,   0, -0.5};
 366   const Double_t xmax0[nBinsSparse0]  = {  100,5000,fNTrigger-0.5};
 367
 368
 369   fhnEvent = new THnSparseF("fhnEvent",";cent;mult",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,10};
 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[3] = {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       fhnEvent->Fill(var0);
 905     }
 906   }
 907   // the loops for rec and gen should be indentical... pass it to a separate
 908   // function ...
 909   // Jet Loop
 910   // Track Loop
 911   // Jet Jet Loop
 912   // Jet Track loop
 913
 914   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 915   FillJetHistos(recJetsList,recParticles,kJetRecFull);
 916   FillTrackHistos(recParticles,kJetRec);
 917
 918   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 919   FillJetHistos(genJetsList,genParticles,kJetGenFull);
 920   FillTrackHistos(genParticles,kJetGen);
 921
 922   // Here follows the jet matching part
 923   // delegate to a separated method?
 924
 925   if(fDoMatching){
 926     FillMatchHistos(recJetsList,genJetsList);
 927   }
 928
 929   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 930   PostData(1, fHistList);
 931 }
 932
 933 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 934
 935   if(iType>=kJetTypes){
 936     return;
 937   }
 938   if(!fFlagJetType[iType])return;
 939
 940   Int_t refMult = fMultRec;
 941   if(iType==kJetGen||iType==kJetGenFull){
 942     refMult = fMultGen;
 943   }
 944
 945   Int_t nJets = jetsList.GetEntries();
 946   fh1NJets[iType]->Fill(nJets);
 947
 948   if(nJets<=0)return;
 949
 950   Float_t ptOld = 1.E+32;
 951   Float_t pT0 = 0;
 952   Float_t pT1 = 0;
 953   Float_t phi0 = 0;
 954   Float_t phi1 = 0;
 955   Int_t ij0 = -1;
 956   Int_t ij1 = -1;
 957
 958   Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
 959   Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
 960
 961   var1[2] = fCentrality;
 962   var1b[1] = fCentrality;
 963   var1[3] = refMult;
 964
 965
 966
 967   Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
 968   var2[2] = fCentrality;
 969
 970   for(int ij = 0;ij < nJets;ij++){
 971     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
 972     Float_t ptJet = jet->Pt();
 973
 974
 975     if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
 976
 977     var1b[0] = ptJet;
 978     if(jet->Trigger()&fJetTriggerBestMask){
 979       fh1PtJetsInBest[iType]->Fill(ptJet);
 980       for(int it = 0;it <fNTrigger;it++){
 981         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 982           var1b[2] = it;
 983           fhnJetPtBest[iType]->Fill(var1b);
 984         }
 985       }
 986     }
 987     if(jet->Trigger()&fJetTriggerExcludeMask){
 988       fh1PtJetsInRej[iType]->Fill(ptJet);
 989       for(int it = 0;it <fNTrigger;it++){
 990         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 991           var1b[2] = it;
 992           fhnJetPtRej[iType]->Fill(var1b);
 993         }
 994       }
 995       continue;
 996     }
 997     fh1PtJetsIn[iType]->Fill(ptJet);
 998     ptOld = ptJet;
 999
1000     // find the dijets assume sorting and acceptance cut...
1001     if(ij==0){
1002       ij0 = ij;
1003       pT0 = ptJet;
1004       phi0 = jet->Phi();
1005       if(phi0<0)phi0+=TMath::Pi()*2.;
1006     }
1007     else if(ptJet>pT1){
1008       // take only the backward hemisphere??
1009       phi1 = jet->Phi();
1010       if(phi1<0)phi1+=TMath::Pi()*2.;
1011       Float_t dPhi = phi1 - phi0;
1012       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1013       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1014       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1015         ij1 = ij;
1016         pT1 = ptJet;
1017       }
1018     }
1019     // fill jet histos for kmax jets
1020
1021     Float_t phiJet = jet->Phi();
1022     Float_t etaJet = jet->Eta();
1023     if(phiJet<0)phiJet+=TMath::Pi()*2.;
1024     fh1TmpRho->Reset();
1025     if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1026
1027     fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1028     // fill leading jets...
1029     AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1030     //      AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1031     Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1032     Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1033
1034     var1[1] = ptJet;
1035     var1[4] = phiBin;
1036     var1[5] = jet->EffectiveAreaCharged();
1037     var1[7] = ptLead;
1038     var1[8] = CheckAcceptance(phiJet,etaJet);
1039
1040     var2[1] = ptJet;
1041     var2[3] = etaJet;
1042     var2[4] = phiJet;
1043     var2[5] = jet->EffectiveAreaCharged();
1044     var2[7] = var1[8];
1045     var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1046
1047     if(ij<kMaxJets){
1048       fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1049       var1[0] = ij;
1050       var2[0] = ij;
1051       for(int it = 0;it <fNTrigger;it++){
1052         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1053           var1[6] = it;
1054           var2[6] = it;
1055           fhnJetPt[iType]->Fill(var1);
1056           fhnJetPtQA[iType]->Fill(var2);
1057         }
1058       }
1059     }
1060     var1[0] = kMaxJets;// fill for all jets
1061     var2[0] = kMaxJets;// fill for all jets
1062     for(int it = 0;it <fNTrigger;it++){
1063       if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1064         var1[6] = it;
1065         fhnJetPt[iType]->Fill(var1);
1066         fhnJetPtQA[iType]->Fill(var2);
1067       }
1068     }
1069
1070     fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
1071
1072     if(particlesList.GetSize()&&ij<kMaxJets){
1073       // Particles... correlated with jets...
1074       for(int it = 0;it<particlesList.GetEntries();++it){
1075         AliVParticle *part = (AliVParticle*)particlesList.At(it);
1076         Float_t deltaR = jet->DeltaR(part);
1077         if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1078       }
1079       // fill the jet shapes
1080     }// if we have particles
1081   }// Jet Loop
1082
1083
1084   // do something with dijets...
1085   if(ij0>=0&&ij1>0){
1086     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1087     Double_t ptJet0 = jet0->Pt();
1088     Double_t phiJet0 = jet0->Phi();
1089     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1090
1091     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1092     Double_t ptJet1 = jet1->Pt();
1093     Double_t phiJet1 = jet1->Phi();
1094     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1095
1096     Float_t deltaPhi = phiJet0 - phiJet1;
1097     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1098     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1099     deltaPhi = TMath::Abs(deltaPhi);
1100     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1101
1102     Float_t asym = 9999;
1103     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1104     fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1105     fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1106     fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1107     Float_t minv = 2.*(jet0->P()*jet1->P()-
1108                        jet0->Px()*jet1->Px()-
1109                        jet0->Py()*jet1->Py()-
1110                        jet0->Pz()*jet1->Pz());    // assume mass == 0;
1111     if(minv<0)minv=0; // prevent numerical instabilities
1112     minv = TMath::Sqrt(minv);
1113     fh1DijetMinv[iType]->Fill(minv);
1114   }
1115
1116
1117
1118   // count down the jets above thrueshold
1119   Int_t nOver = nJets;
1120   if(nOver>0){
1121     TIterator *jetIter = jetsList.MakeIterator();
1122     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1123     if(tmpJet){
1124       Float_t pt = tmpJet->Pt();
1125       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1126         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1127         while(pt<ptCut&&tmpJet){
1128           nOver--;
1129           tmpJet = (AliAODJet*)(jetIter->Next());
1130           if(tmpJet){
1131             pt = tmpJet->Pt();
1132           }
1133         }
1134         if(nOver<=0)break;
1135         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1136       }
1137     }
1138     delete jetIter;
1139   }
1140 }
1141
1142 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1143
1144   if(fFlagJetType[iType]<=0)return;
1145   Int_t refMult = fMultRec;
1146   if(iType==kJetGen||iType==kJetGenFull){
1147     refMult = fMultGen;
1148
1149   }
1150
1151   //
1152   Double_t var3[6]; // track number;p_{T};cent;#tracks;RP
1153   var3[2] = fCentrality;
1154   var3[3] = refMult;
1155   Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
1156   var4[2] = fCentrality;
1157   Int_t nTrackOver = particlesList.GetSize();
1158   // do the same for tracks and jets
1159   if(nTrackOver>0){
1160     TIterator *trackIter = particlesList.MakeIterator();
1161     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1162     Float_t pt = tmpTrack->Pt();
1163     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1164       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1165       while(pt<ptCut&&tmpTrack){
1166         nTrackOver--;
1167         tmpTrack = (AliVParticle*)(trackIter->Next());
1168         if(tmpTrack){
1169           pt = tmpTrack->Pt();
1170         }
1171       }
1172       if(nTrackOver<=0)break;
1173       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1174     }
1175
1176
1177     trackIter->Reset();
1178     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1179     Float_t sumPt = 0;
1180
1181     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1182       Float_t tmpPt = tmpTrack->Pt();
1183       fh1PtTracksIn[iType]->Fill(tmpPt);
1184       fh1PtTracksInLow[iType]->Fill(tmpPt);
1185
1186       sumPt += tmpPt;
1187       Float_t tmpPhi = tmpTrack->Phi();
1188       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1189
1190
1191       Float_t phiRP = tmpPhi-fRPAngle;
1192       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1193       if(phiRP<0)phiRP += TMath::Pi();
1194       if(phiRP<0)phiRP += TMath::Pi();
1195       const float allPhi = -1./180.*TMath::Pi();
1196
1197       if(tmpPt<100){
1198         fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1199         fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1200
1201         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1202         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1203       }
1204       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1205       var3[0] = 1;
1206       var3[1] = tmpPt;
1207       var3[4] = phiBin;
1208
1209       var4[0] = 1;
1210       var4[1] = tmpPt;
1211       var4[3] = tmpTrack->Eta();
1212       var4[4] = tmpPhi;
1213       var4[5] = tmpTrack->Charge();
1214
1215
1216       for(int it = 0;it <fNTrigger;it++){
1217         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1218           var3[0] = 1;
1219           var3[5] = it;
1220           fhnTrackPt[iType]->Fill(var3);
1221           if(tmpTrack==leading){
1222             var3[0] = 0;
1223             fhnTrackPt[iType]->Fill(var3);
1224           }
1225         }
1226       }
1227
1228
1229
1230
1231
1232       fhnTrackPtQA[iType]->Fill(var4);
1233
1234       if(tmpTrack==leading){
1235         var4[0] = 0;
1236         fhnTrackPtQA[iType]->Fill(var4);
1237         continue;
1238       }
1239     }
1240     fh1SumPtTrack[iType]->Fill(sumPt);
1241     delete trackIter;
1242   }
1243
1244 }
1245
1246
1247 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1248
1249
1250   // Fill al the matching histos
1251   // It is important that the acceptances for the mathing are as large as possible
1252   // to avoid false matches on the edge of acceptance
1253   // therefore we add some extra matching jets as overhead
1254
1255   static TArrayI aGenIndex(recJetsList.GetEntries());
1256   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1257
1258   static TArrayI aRecIndex(genJetsList.GetEntries());
1259   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1260
1261   if(fDebug){
1262     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1263     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1264   }
1265   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1266                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1267                                             aGenIndex,aRecIndex,fDebug);
1268
1269   if(fDebug){
1270     for(int i = 0;i< aGenIndex.GetSize();++i){
1271       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1272     }
1273     for(int i = 0;i< aRecIndex.GetSize();++i){
1274       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1275     }
1276   }
1277
1278   Double_t container[8];
1279   Double_t containerRec[4];
1280   Double_t containerGen[4];
1281
1282   // loop over generated jets
1283   // consider the
1284   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1285     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1286     Double_t ptGen = genJet->Pt();
1287     Double_t phiGen = genJet->Phi();
1288     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1289     Double_t etaGen = genJet->Eta();
1290     containerGen[0] = ptGen;
1291     containerGen[1] = etaGen;
1292     containerGen[2] = phiGen;
1293     containerGen[3] = genJet->GetPtLeading();
1294
1295     fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1296
1297     Int_t ir = aRecIndex[ig];
1298     if(ir>=0&&ir<recJetsList.GetEntries()){
1299       fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1300
1301       if(JetSelected(genJet)){
1302         fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1303
1304         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1305
1306         fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1307         if(JetSelected(recJet)) {
1308           fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1309
1310           containerGen[3] = recJet->GetPtLeading();
1311           fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner in eta window with leading track on reconstructed level
1312         }
1313       }
1314     }
1315   }// loop over generated jets used for matching...
1316
1317
1318
1319   // loop over reconstructed jets
1320   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1321     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1322     Double_t etaRec = recJet->Eta();
1323     Double_t ptRec = recJet->Pt();
1324     Double_t phiRec = recJet->Phi();
1325     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1326
1327     containerRec[0] = ptRec;
1328     containerRec[1] = etaRec;
1329     containerRec[2] = phiRec;
1330     containerRec[3] = recJet->GetPtLeading();
1331
1332     container[0] = ptRec;
1333     container[1] = etaRec;
1334     container[2] = phiRec;
1335     container[3] = recJet->GetPtLeading();
1336
1337     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1338
1339     if(JetSelected(recJet)){
1340       fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
1341       // Fill Correlation
1342       Int_t ig = aGenIndex[ir];
1343       if(ig>=0 && ig<genJetsList.GetEntries()){
1344         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1345         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1346         Double_t ptGen  = genJet->Pt();
1347         Double_t phiGen = genJet->Phi();
1348         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1349         Double_t etaGen = genJet->Eta();
1350
1351         containerGen[0] = ptGen;
1352         containerGen[1] = etaGen;
1353         containerGen[2] = phiGen;
1354         containerGen[3] = genJet->GetPtLeading();
1355
1356         container[4] = ptGen;
1357         container[5] = etaGen;
1358         container[6] = phiGen;
1359         container[7] = genJet->GetPtLeading();
1360
1361         fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1362
1363         fhnCorrelation->Fill(container);
1364         if(ptGen>0){
1365           Float_t delta = (ptRec-ptGen)/ptGen;
1366           fh2RelPtFGen->Fill(ptGen,delta);
1367           fh2PtFGen->Fill(ptGen,ptRec);
1368         }
1369       }
1370     }// loop over reconstructed jets
1371   }
1372   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1373 }
1374
1375
1376 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1377   //
1378   // Create the particle container for the correction framework manager and
1379   // link it
1380   //
1381   const Int_t kNvar   = 4 ; //number of variables on the grid:pt,eta, phi
1382   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1383   const Double_t kEtamin = -1.0, kEtamax = 1.0;
1384   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1385   const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
1386
1387   // can we neglect migration in eta and phi?
1388   // phi should be no problem since we cover full phi and are phi symmetric
1389   // eta migration is more difficult  due to needed acceptance correction
1390   // in limited eta range
1391
1392   //arrays for the number of bins in each dimension
1393   Int_t iBin[kNvar];
1394   iBin[0] = 125; //bins in pt
1395   iBin[1] =  4; //bins in eta
1396   iBin[2] =  1; // bins in phi
1397   iBin[3] =  10; // bins in leading track Pt
1398
1399   //arrays for lower bounds :
1400   Double_t* binEdges[kNvar];
1401   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1402     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1403
1404   //values for bin lower bounds
1405   //  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);
1406   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1407   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1408   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1409   binEdges[3][0]= kPtLeadingTrackPtMin;
1410   binEdges[3][1]= 1.;
1411   binEdges[3][2]= 2.;
1412   binEdges[3][3]= 3.;
1413   binEdges[3][4]= 4.;
1414   binEdges[3][5]= 5.;
1415   binEdges[3][6]= 6.;
1416   binEdges[3][7]= 8.;
1417   binEdges[3][8]= 10.;
1418   binEdges[3][9]= 12.;
1419   binEdges[3][10]= kPtLeadingTrackPtMax;
1420
1421   fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
1422   for (int k=0; k<kNvar; k++) {
1423     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1424   }
1425
1426   //create correlation matrix for unfolding
1427   Int_t thnDim[2*kNvar];
1428   for (int k=0; k<kNvar; k++) {
1429     //first half  : reconstructed
1430     //second half : MC
1431     thnDim[k]      = iBin[k];
1432     thnDim[k+kNvar] = iBin[k];
1433   }
1434
1435   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
1436   for (int k=0; k<kNvar; k++) {
1437     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1438     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1439   }
1440
1441   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1442     delete [] binEdges[ivar];
1443
1444
1445 }
1446
1447 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1448 {
1449   // Terminate analysis
1450   //
1451   if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1452 }
1453
1454
1455 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1456
1457   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1458
1459   Int_t iCount = 0;
1460   if(type==kTrackAOD){
1461     AliAODEvent *aod = 0;
1462     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1463     else aod = AODEvent();
1464     if(!aod){
1465       return iCount;
1466     }
1467     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1468       AliAODTrack *tr = aod->GetTrack(it);
1469       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1470       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1471       if(tr->Pt()<fMinTrackPt)continue;
1472       if(fDebug>0){
1473         if(tr->Pt()>20){
1474           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1475           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1476           tr->Print();
1477           //    tr->Dump();
1478           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1479           if(fESD){
1480             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1481             esdTr->Print("");
1482             // esdTr->Dump();
1483           }
1484         }
1485       }
1486       list->Add(tr);
1487       iCount++;
1488     }
1489   }
1490   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1491     AliMCEvent* mcEvent = MCEvent();
1492     if(!mcEvent)return iCount;
1493     // we want to have alivpartilces so use get track
1494     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1495       if(!mcEvent->IsPhysicalPrimary(it))continue;
1496       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1497       if(part->Pt()<fMinTrackPt)continue;
1498       if(type == kTrackKineAll){
1499         list->Add(part);
1500         iCount++;
1501       }
1502       else if(type == kTrackKineCharged){
1503         if(part->Particle()->GetPDG()->Charge()==0)continue;
1504         list->Add(part);
1505         iCount++;
1506       }
1507     }
1508   }
1509   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1510     AliAODEvent *aod = 0;
1511     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1512     else aod = AODEvent();
1513     if(!aod)return iCount;
1514     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1515     if(!tca)return iCount;
1516     for(int it = 0;it < tca->GetEntriesFast();++it){
1517       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1518       if(!part)continue;
1519       if(part->Pt()<fMinTrackPt)continue;
1520       if(!part->IsPhysicalPrimary())continue;
1521       if(type == kTrackAODMCAll){
1522         list->Add(part);
1523         iCount++;
1524       }
1525       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1526         if(part->Charge()==0)continue;
1527         if(kTrackAODMCCharged){
1528           list->Add(part);
1529         }
1530         else {
1531           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1532           list->Add(part);
1533         }
1534         iCount++;
1535       }
1536     }
1537   }// AODMCparticle
1538   list->Sort();
1539   return iCount;
1540
1541 }
1542
1543
1544 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1545   AliAODEvent *aod = 0;
1546   if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1547   else aod = AODEvent();
1548   if(!aod){
1549     return 101;
1550   }
1551   return aod->GetHeader()->GetCentrality();
1552 }
1553
1554
1555
1556 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1557   Bool_t selected = false;
1558
1559   if(!jet)return selected;
1560
1561   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1562     selected = kTRUE;
1563   }
1564   return selected;
1565
1566 }
1567
1568 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1569
1570   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1571   Int_t iCount = 0;
1572
1573   if(!jarray){
1574     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1575     return iCount;
1576   }
1577
1578
1579   for(int ij=0;ij<jarray->GetEntries();ij++){
1580     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1581     if(!jet)continue;
1582     if(type==0){
1583       // No cut at all, main purpose here is sorting
1584       list->Add(jet);
1585       iCount++;
1586     }
1587     else if(type == 1){
1588       // eta cut
1589       if(JetSelected(jet)){
1590         list->Add(jet);
1591         iCount++;
1592       }
1593     }
1594   }
1595
1596   list->Sort();
1597   return iCount;
1598
1599 }
1600
1601
1602 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1603   if(!jets)return 0;
1604
1605   Int_t refMult = 0;
1606   for(int ij = 0;ij < jets->GetEntries();++ij){
1607     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1608     if(!jet)continue;
1609     TRefArray *refs = jet->GetRefTracks();
1610     if(!refs)continue;
1611     refMult += refs->GetEntries();
1612   }
1613   return refMult;
1614
1615 }
1616
1617
1618 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1619   if(!jet)return 0;
1620   TRefArray *refs = jet->GetRefTracks();
1621   if(!refs) return 0;
1622   AliVParticle *leading = 0;
1623   Float_t fMaxPt = 0;
1624   for(int i = 0;i<refs->GetEntries();i++){
1625     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1626     if(!tmp)continue;
1627     if(tmp->Pt()>fMaxPt){
1628       leading = tmp;
1629       fMaxPt = tmp->Pt();
1630     }
1631   }
1632   return leading;
1633 }
1634
1635
1636 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1637   if(!jet)return 0;
1638   if(!list) return 0;
1639   AliVParticle *leading = 0;
1640   Float_t fMaxPt = 0;
1641   for(int i = 0;i<list->GetEntries();i++){
1642     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1643     if(!tmp)continue;
1644     if(jet->DeltaR(tmp)>r)continue;
1645     if(tmp->Pt()>fMaxPt){
1646       leading = tmp;
1647       fMaxPt = tmp->Pt();
1648     }
1649   }
1650   return leading;
1651 }
1652
1653 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1654 {
1655   Int_t phibin=-1;
1656   if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1657   Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1658   phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1659   if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1660   return phibin;
1661 }
1662
1663 void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1664   //
1665   // set number of trigger bins
1666   //
1667   if(n>0){
1668     fNTrigger = n;
1669     delete [] fTriggerName;
1670     fTriggerName = new TString [fNTrigger];
1671     delete [] fTriggerBit;fTriggerBit = 0;
1672     fTriggerBit = new UInt_t [fNTrigger];
1673   }
1674   else{
1675     fNTrigger = 0;
1676   }
1677 }
1678
1679
1680 void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1681   //
1682   // set trigger bin
1683   //
1684   if(i<fNTrigger){
1685     fTriggerBit[i] = it;
1686     fTriggerName[i] =  c;
1687   }
1688 }
1689
1690
1691 void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1692   //
1693   // Set number of acceptance bins
1694   //
1695
1696   if(n>0){
1697     fNAcceptance = n;
1698     delete [] fAcceptancePhiMin;
1699     delete [] fAcceptancePhiMax;
1700     delete [] fAcceptanceEtaMin;
1701     delete [] fAcceptanceEtaMax;
1702
1703     fAcceptancePhiMin = new Float_t[fNAcceptance];
1704     fAcceptancePhiMax = new Float_t[fNAcceptance];
1705     fAcceptanceEtaMin = new Float_t[fNAcceptance];
1706     fAcceptanceEtaMax = new Float_t[fNAcceptance];
1707   }
1708   else{
1709     fNTrigger = 0;
1710   }
1711 }
1712
1713 void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1714   //
1715   // Set acceptance bins
1716   //
1717
1718   if(i<fNAcceptance){
1719     Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1720
1721     fAcceptancePhiMin[i] = phiMin;
1722     fAcceptancePhiMax[i] = phiMax;
1723     fAcceptanceEtaMin[i] = etaMin;
1724     fAcceptanceEtaMax[i] = etaMax;
1725   }
1726   else{
1727     fNTrigger = 0;
1728   }
1729 }
1730
1731
1732 Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1733
1734   //
1735   // return aceptnace bin do no use ovelapping bins
1736   //
1737
1738   for(int i = 0;i<fNAcceptance;i++){
1739     if(phi<fAcceptancePhiMin[i])continue;
1740     if(phi>fAcceptancePhiMax[i])continue;
1741     if(eta<fAcceptanceEtaMin[i])continue;
1742     if(eta>fAcceptanceEtaMax[i])continue;
1743     return i;
1744   }
1745   // catch the rest
1746   return fNAcceptance;
1747 }
1748
1749
1750
1751
1752
1753 AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1754   //
1755   delete [] fTriggerBit;
1756   delete [] fTriggerName;
1757
1758   delete [] fAcceptancePhiMin;
1759   delete [] fAcceptancePhiMax;
1760   delete [] fAcceptanceEtaMin;
1761   delete [] fAcceptanceEtaMax;
1762
1763
1764
1765 }