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 = 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.05,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.   total bins = 9.72 M
 547     const Int_t nBinsSparse2 = 9;
 548     Int_t nBins2[nBinsSparse2] = {     kMaxJets+1, 25,   8,  18,             180, 10,fNTrigger,fNAcceptance+0.5,10};
 549     if(cJetBranch.Contains("RandomCone")){
 550       nBins2[5] = 1;
 551     }
 552     const Double_t xmin2[nBinsSparse2]  = {        -0.5,  0,   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     fHistList->Add(fhnJetPtQA[ij]);
 556
 557     // Bins:track number  pTtrack, cent, mult, RP.   total bins = 224 k
 558     const Int_t nBinsSparse3 = 6;
 559     const Int_t nRPBinsSparse3 = 3;
 560     const Int_t nBins3[nBinsSparse3] = {       2,    100,     10,   1,    nRPBinsSparse3,fNTrigger};
 561     const Double_t xmin3[nBinsSparse3]  = { -0.5,      0,   0,      0,        -0.5,-0.5};
 562     const Double_t xmax3[nBinsSparse3]  = { 1.5,     200, 100,   4000,nRPBinsSparse3-0.5,fNTrigger-0.5};
 563
 564       // change the binning ot the pT axis:
 565     Double_t *xPt3 = new Double_t[nBins3[1]+1];
 566     xPt3[0] = 0.;
 567     for(int i = 1; i<=nBins3[1];i++){
 568       if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
 569       else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
 570       else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
 571       else if(xPt3[i-1]<15)xPt3[i] = xPt3[i-1] +  1.; // 62 - 67
 572       else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 2.; // 67 - 72
 573       else if(xPt3[i-1]<60)xPt3[i] = xPt3[i-1] + 5; // 72 - 76
 574       else xPt3[i] = xPt3[i-1] + 10; // 76 - 100 = 140
 575     }
 576
 577       fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP;trigger",nBinsSparse3,nBins3,xmin3,xmax3);
 578       fhnTrackPt[ij]->SetBinEdges(1,xPt3);
 579       fHistList->Add(fhnTrackPt[ij]);
 580       delete [] xPt3;
 581
 582       // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M
 583       const Int_t nBinsSparse4 = 6;
 584       const Int_t nBins4[nBinsSparse4] =    {    2, 50,  10,  20, 180,2};
 585       const Double_t xmin4[nBinsSparse4]  = { -0.5,  0,   0, -1.0,   0.,-1.5};
 586       const Double_t xmax4[nBinsSparse4]  = {  1.5,150, 100,  1.0,2.*TMath::Pi(),1.5};
 587
 588       // change the binning ot the pT axis:
 589       Double_t *xPt4 = new Double_t[nBins4[1]+1];
 590       xPt4[0] = 0.;
 591       for(int i = 1; i<=nBins4[1];i++){
 592         if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1;
 593         else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5;
 594         else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] +  1.;
 595         else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] +  2.5;
 596         else xPt4[i] = xPt4[i-1] + 5.;
 597       }
 598       fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi;sign",nBinsSparse4,nBins4,xmin4,xmax4);
 599       fhnTrackPtQA[ij]->SetBinEdges(1,xPt4);
 600       fHistList->Add(fhnTrackPtQA[ij]);
 601       delete [] xPt4;
 602
 603       for(int i = 0;i <= kMaxJets;++i){
 604         fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
 605         fHistList->Add(fh1PtIn[ij][i]);
 606
 607
 608         if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
 609         fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
 610                                            Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
 611                                                 cAdd.Data()),
 612                                            200,0.,200.,nBinPt,binLimitsPt);
 613         fHistList->Add(fh2LTrackPtJetPt[ij][i]);
 614       }
 615
 616
 617       fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
 618       fHistList->Add(fh1DijetMinv[ij]);
 619
 620       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);
 621       fHistList->Add(fh2DijetDeltaPhiPt[ij]);
 622
 623       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);
 624       fHistList->Add(fh2DijetAsymPt[ij]);
 625
 626       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.);
 627       fHistList->Add(fh2DijetPt2vsPt1[ij]);
 628       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.);
 629       fHistList->Add( fh2DijetDifvsSum[ij]);
 630     }
 631   // =========== Switch on Sumw2 for all histos ===========
 632   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 633     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 634     if (h1){
 635       h1->Sumw2();
 636       continue;
 637     }
 638     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 639     if(hn)hn->Sumw2();
 640   }
 641   TH1::AddDirectory(oldStatus);
 642 }
 643
 644 void AliAnalysisTaskJetSpectrum2::Init()
 645 {
 646   //
 647   // Initialization
 648   //
 649
 650   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 651
 652 }
 653
 654 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
 655
 656
 657   Bool_t selected = kTRUE;
 658
 659   if(fUseGlobalSelection&&fEventSelectionMask==0){
 660     selected = AliAnalysisHelperJetTasks::Selected();
 661   }
 662   else if(fUseGlobalSelection&&fEventSelectionMask>0){
 663     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
 664   }
 665
 666   if(fEventClass>0){
 667     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
 668   }
 669
 670   if(!selected){
 671     // no selection by the service task, we continue
 672     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
 673     PostData(1, fHistList);
 674     return;
 675   }
 676
 677
 678   static AliAODEvent* aod = 0;
 679
 680   // take all other information from the aod we take the tracks from
 681   if(!aod){
 682    if(fUseAODTrackInput)aod = fAODIn;
 683    else aod = fAODOut;
 684   }
 685
 686
 687   //
 688   // Execute analysis for current event
 689   //
 690   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
 691   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 692   if(!aodH){
 693     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 694     return;
 695   }
 696
 697   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 698   TClonesArray *aodRecJets = 0;
 699
 700   if(fAODOut&&!aodRecJets){
 701     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
 702   }
 703   if(fAODExtension&&!aodRecJets){
 704     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
 705   }
 706   if(fAODIn&&!aodRecJets){
 707     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
 708   }
 709
 710
 711
 712   if(!aodRecJets){
 713     if(fDebug){
 714
 715         Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 716         if(fAODIn){
 717           Printf("Input AOD >>>>");
 718           fAODIn->Print();
 719         }
 720         if(fAODExtension){
 721           Printf("AOD Extension >>>>");
 722           fAODExtension->Print();
 723         }
 724         if(fAODOut){
 725           Printf("Output AOD >>>>");
 726           fAODOut->Print();
 727         }
 728     }
 729         return;
 730   }
 731
 732   TClonesArray *aodGenJets = 0;
 733   if(fBranchGen.Length()>0){
 734     if(fAODOut&&!aodGenJets){
 735       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
 736     }
 737     if(fAODExtension&&!aodGenJets){
 738       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
 739     }
 740     if(fAODIn&&!aodGenJets){
 741       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
 742     }
 743
 744     if(!aodGenJets){
 745       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
 746       return;
 747     }
 748   }
 749
 750   TClonesArray *aodBackRecJets = 0;
 751   if(fBranchBkgRec.Length()>0){
 752     if(fAODOut&&!aodBackRecJets){
 753       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
 754     }
 755     if(fAODExtension&&!aodBackRecJets){
 756       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
 757     }
 758     if(fAODIn&&!aodBackRecJets){
 759       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
 760     }
 761
 762     if(!aodBackRecJets){
 763       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
 764       return;
 765     }
 766   }
 767
 768
 769   TClonesArray *aodBackGenJets = 0;
 770
 771   if(fBranchBkgGen.Length()>0){
 772     if(fAODOut&&!aodBackGenJets){
 773       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
 774     }
 775     if(fAODExtension&&!aodBackGenJets){
 776       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
 777     }
 778     if(fAODIn&&!aodBackGenJets){
 779       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
 780     }
 781
 782     if(!aodBackGenJets){
 783       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
 784       return;
 785     }
 786   }
 787
 788
 789   // new Scheme
 790   // first fill all the pure  histograms, i.e. full jets
 791   // in case of the correaltion limit to the n-leading jets
 792
 793   // reconstructed
 794
 795
 796   // generated
 797
 798
 799   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
 800
 801
 802
 803   TList genJetsList;         // full acceptance
 804   TList genJetsListCut;      // acceptance cut
 805   TList recJetsList;         // full acceptance
 806   TList recJetsListCut;      // acceptance cut
 807
 808
 809   GetListOfJets(&recJetsList,aodRecJets,0);
 810   GetListOfJets(&recJetsListCut,aodRecJets,1);
 811
 812   if(fBranchGen.Length()>0){
 813     GetListOfJets(&genJetsList,aodGenJets,0);
 814     GetListOfJets(&genJetsListCut,aodGenJets,1);
 815   }
 816
 817   Double_t eventW = 1;
 818   Double_t ptHard = 0;
 819   Double_t nTrials = 1; // Trials for MC trigger
 820   fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
 821
 822   // Getting some global properties
 823   fCentrality = GetCentrality();
 824   if(fCentrality<=0)fCentrality = 0;
 825   fh1Centrality->Fill(fCentrality);
 826
 827
 828   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
 829     // this is the part we only use when we have MC information
 830     AliMCEvent* mcEvent = MCEvent();
 831     //    AliStack *pStack = 0;
 832     if(!mcEvent){
 833       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 834       return;
 835     }
 836     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 837     if(pythiaGenHeader){
 838       nTrials = pythiaGenHeader->Trials();
 839       ptHard  = pythiaGenHeader->GetPtHard();
 840       int iProcessType = pythiaGenHeader->ProcessType();
 841       // 11 f+f -> f+f
 842       // 12 f+barf -> f+barf
 843       // 13 f+barf -> g+g
 844       // 28 f+g -> f+g
 845       // 53 g+g -> f+barf
 846       // 68 g+g -> g+g
 847       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 848       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 849     }
 850   }// (fAnalysisType&kMCESD)==kMCESD)
 851
 852
 853   // we simply fetch the tracks/mc particles as a list of AliVParticles
 854
 855   TList recParticles;
 856   TList genParticles;
 857
 858   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 859   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 860   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 861   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 862
 863   //  CalculateReactionPlaneAngle(&recParticles);
 864   fRPAngle = 0;
 865
 866   if(fRPMethod==0)fRPAngle = aod->GetHeader()->GetEventplane();
 867   else if(fRPMethod==1||fRPMethod==2){
 868     fRPAngle = aod->GetHeader()->GetQTheta(fRPMethod);
 869   }
 870   fh1RP->Fill(fRPAngle);
 871   fh2RPCentrality->Fill(fCentrality,fRPAngle);
 872   // Event control and counting ...
 873   // MC
 874   fh1PtHard->Fill(ptHard,eventW);
 875   fh1PtHardNoW->Fill(ptHard,1);
 876   fh1PtHardTrials->Fill(ptHard,nTrials);
 877
 878   // Real
 879   if(aod->GetPrimaryVertex()){// No vtx for pure MC
 880     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
 881   }
 882
 883
 884   Int_t recMult1 = recParticles.GetEntries();
 885   Int_t genMult1 = genParticles.GetEntries();
 886
 887   Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
 888   Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
 889
 890   fh2MultRec->Fill(recMult1,recMult2);
 891   fh2MultGen->Fill(genMult1,genMult2);
 892   fMultRec = recMult1;
 893   if(fMultRec<=0)fMultRec = recMult2;
 894   fMultGen = genMult1;
 895   if(fMultGen<=0)fMultGen = genMult2;
 896
 897   Double_t var0[3] = {0,};
 898   var0[0] = fCentrality;
 899   var0[1] = fMultRec;
 900   for(int it=0;it<fNTrigger;it++){
 901     if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 902       var0[2] = it;
 903       fhnEvent->Fill(var0);
 904     }
 905   }
 906   // the loops for rec and gen should be indentical... pass it to a separate
 907   // function ...
 908   // Jet Loop
 909   // Track Loop
 910   // Jet Jet Loop
 911   // Jet Track loop
 912
 913   FillJetHistos(recJetsListCut,recParticles,kJetRec);
 914   FillJetHistos(recJetsList,recParticles,kJetRecFull);
 915   FillTrackHistos(recParticles,kJetRec);
 916
 917   FillJetHistos(genJetsListCut,genParticles,kJetGen);
 918   FillJetHistos(genJetsList,genParticles,kJetGenFull);
 919   FillTrackHistos(genParticles,kJetGen);
 920
 921   // Here follows the jet matching part
 922   // delegate to a separated method?
 923
 924   if(fDoMatching){
 925     FillMatchHistos(recJetsList,genJetsList);
 926   }
 927
 928   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 929   PostData(1, fHistList);
 930 }
 931
 932 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
 933
 934   if(iType>=kJetTypes){
 935     return;
 936   }
 937   if(!fFlagJetType[iType])return;
 938
 939   Int_t refMult = fMultRec;
 940   if(iType==kJetGen||iType==kJetGenFull){
 941     refMult = fMultGen;
 942   }
 943
 944   Int_t nJets = jetsList.GetEntries();
 945   fh1NJets[iType]->Fill(nJets);
 946
 947   if(nJets<=0)return;
 948
 949   Float_t ptOld = 1.E+32;
 950   Float_t pT0 = 0;
 951   Float_t pT1 = 0;
 952   Float_t phi0 = 0;
 953   Float_t phi1 = 0;
 954   Int_t ij0 = -1;
 955   Int_t ij1 = -1;
 956
 957   Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
 958   Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
 959
 960   var1[2] = fCentrality;
 961   var1b[1] = fCentrality;
 962   var1[3] = refMult;
 963
 964
 965
 966   Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
 967   var2[2] = fCentrality;
 968
 969   for(int ij = 0;ij < nJets;ij++){
 970     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
 971     Float_t ptJet = jet->Pt();
 972
 973
 974     if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
 975
 976     var1b[0] = ptJet;
 977     if(jet->Trigger()&fJetTriggerBestMask){
 978       fh1PtJetsInBest[iType]->Fill(ptJet);
 979       for(int it = 0;it <fNTrigger;it++){
 980         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 981           var1b[2] = it;
 982           fhnJetPtBest[iType]->Fill(var1b);
 983         }
 984       }
 985     }
 986     if(jet->Trigger()&fJetTriggerExcludeMask){
 987       fh1PtJetsInRej[iType]->Fill(ptJet);
 988       for(int it = 0;it <fNTrigger;it++){
 989         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
 990           var1b[2] = it;
 991           fhnJetPtRej[iType]->Fill(var1b);
 992         }
 993       }
 994       continue;
 995     }
 996     fh1PtJetsIn[iType]->Fill(ptJet);
 997     ptOld = ptJet;
 998
 999     // find the dijets assume sorting and acceptance cut...
1000     if(ij==0){
1001       ij0 = ij;
1002       pT0 = ptJet;
1003       phi0 = jet->Phi();
1004       if(phi0<0)phi0+=TMath::Pi()*2.;
1005     }
1006     else if(ptJet>pT1){
1007       // take only the backward hemisphere??
1008       phi1 = jet->Phi();
1009       if(phi1<0)phi1+=TMath::Pi()*2.;
1010       Float_t dPhi = phi1 - phi0;
1011       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1012       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1013       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1014         ij1 = ij;
1015         pT1 = ptJet;
1016       }
1017     }
1018     // fill jet histos for kmax jets
1019
1020       Float_t phiJet = jet->Phi();
1021       Float_t etaJet = jet->Eta();
1022       if(phiJet<0)phiJet+=TMath::Pi()*2.;
1023       fh1TmpRho->Reset();
1024       if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1025
1026       fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1027       // fill leading jets...
1028       AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1029       //      AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1030       Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1031       Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1032
1033       var1[1] = ptJet;
1034       var1[4] = phiBin;
1035       var1[5] = jet->EffectiveAreaCharged();
1036       var1[7] = ptLead;
1037       var1[8] = CheckAcceptance(phiJet,etaJet);
1038
1039       var2[1] = ptJet;
1040       var2[3] = etaJet;
1041       var2[4] = phiJet;
1042       var2[5] = jet->EffectiveAreaCharged();
1043       var2[7] = var1[8];
1044       var2[8] = (leadTrack?leadTrack->Charge()*leadTrack->Pt():0);//pT of leading jet x charge
1045
1046       if(ij<kMaxJets){
1047         fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1048         var1[0] = ij;
1049         var2[0] = ij;
1050         for(int it = 0;it <fNTrigger;it++){
1051           if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1052             var1[6] = it;
1053             var2[6] = it;
1054             fhnJetPt[iType]->Fill(var1);
1055             fhnJetPtQA[iType]->Fill(var2);
1056           }
1057         }
1058       }
1059       var1[0] = kMaxJets;// fill for all jets
1060       var2[0] = kMaxJets;// fill for all jets
1061       for(int it = 0;it <fNTrigger;it++){
1062         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1063           var1[6] = it;
1064           fhnJetPt[iType]->Fill(var1);
1065           fhnJetPtQA[iType]->Fill(var2);
1066         }
1067       }
1068
1069       fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
1070
1071       if(particlesList.GetSize()&&ij<kMaxJets){
1072         // Particles... correlated with jets...
1073         for(int it = 0;it<particlesList.GetEntries();++it){
1074           AliVParticle *part = (AliVParticle*)particlesList.At(it);
1075           Float_t deltaR = jet->DeltaR(part);
1076           if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1077         }
1078         // fill the jet shapes
1079       }// if we have particles
1080   }// Jet Loop
1081
1082
1083   // do something with dijets...
1084   if(ij0>=0&&ij1>0){
1085     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1086     Double_t ptJet0 = jet0->Pt();
1087     Double_t phiJet0 = jet0->Phi();
1088     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1089
1090     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1091     Double_t ptJet1 = jet1->Pt();
1092     Double_t phiJet1 = jet1->Phi();
1093     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1094
1095     Float_t deltaPhi = phiJet0 - phiJet1;
1096     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1097     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1098     deltaPhi = TMath::Abs(deltaPhi);
1099     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1100
1101     Float_t asym = 9999;
1102     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1103       fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1104       fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1105       fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1106       Float_t minv = 2.*(jet0->P()*jet1->P()-
1107                          jet0->Px()*jet1->Px()-
1108                          jet0->Py()*jet1->Py()-
1109                          jet0->Pz()*jet1->Pz());    // assume mass == 0;
1110       if(minv<0)minv=0; // prevent numerical instabilities
1111       minv = TMath::Sqrt(minv);
1112       fh1DijetMinv[iType]->Fill(minv);
1113   }
1114
1115
1116
1117   // count down the jets above thrueshold
1118   Int_t nOver = nJets;
1119   if(nOver>0){
1120     TIterator *jetIter = jetsList.MakeIterator();
1121     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1122     if(tmpJet){
1123       Float_t pt = tmpJet->Pt();
1124       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1125         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1126         while(pt<ptCut&&tmpJet){
1127           nOver--;
1128           tmpJet = (AliAODJet*)(jetIter->Next());
1129           if(tmpJet){
1130             pt = tmpJet->Pt();
1131           }
1132         }
1133         if(nOver<=0)break;
1134         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1135       }
1136     }
1137     delete jetIter;
1138   }
1139 }
1140
1141 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1142
1143   if(fFlagJetType[iType]<=0)return;
1144   Int_t refMult = fMultRec;
1145   if(iType==kJetGen||iType==kJetGenFull){
1146     refMult = fMultGen;
1147
1148   }
1149
1150   //
1151   Double_t var3[6]; // track number;p_{T};cent;#tracks;RP
1152   var3[2] = fCentrality;
1153   var3[3] = refMult;
1154   Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
1155   var4[2] = fCentrality;
1156   Int_t nTrackOver = particlesList.GetSize();
1157   // do the same for tracks and jets
1158   if(nTrackOver>0){
1159     TIterator *trackIter = particlesList.MakeIterator();
1160     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1161     Float_t pt = tmpTrack->Pt();
1162     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1163       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1164       while(pt<ptCut&&tmpTrack){
1165         nTrackOver--;
1166         tmpTrack = (AliVParticle*)(trackIter->Next());
1167         if(tmpTrack){
1168           pt = tmpTrack->Pt();
1169         }
1170       }
1171       if(nTrackOver<=0)break;
1172       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1173     }
1174
1175
1176     trackIter->Reset();
1177     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1178     Float_t sumPt = 0;
1179
1180     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1181       Float_t tmpPt = tmpTrack->Pt();
1182       fh1PtTracksIn[iType]->Fill(tmpPt);
1183       fh1PtTracksInLow[iType]->Fill(tmpPt);
1184
1185       sumPt += tmpPt;
1186       Float_t tmpPhi = tmpTrack->Phi();
1187       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1188
1189
1190       Float_t phiRP = tmpPhi-fRPAngle;
1191       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1192       if(phiRP<0)phiRP += TMath::Pi();
1193       if(phiRP<0)phiRP += TMath::Pi();
1194       const float allPhi = -1./180.*TMath::Pi();
1195
1196       if(tmpPt<100){
1197         fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1198         fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1199
1200         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1201         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1202       }
1203       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1204       var3[0] = 1;
1205       var3[1] = tmpPt;
1206       var3[4] = phiBin;
1207
1208       var4[0] = 1;
1209       var4[1] = tmpPt;
1210       var4[3] = tmpTrack->Eta();
1211       var4[4] = tmpPhi;
1212       var4[5] = tmpTrack->Charge();
1213
1214
1215       for(int it = 0;it <fNTrigger;it++){
1216         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1217           var3[0] = 1;
1218           var3[5] = it;
1219           fhnTrackPt[iType]->Fill(var3);
1220           if(tmpTrack==leading){
1221             var3[0] = 0;
1222             fhnTrackPt[iType]->Fill(var3);
1223           }
1224         }
1225       }
1226
1227
1228
1229
1230
1231       fhnTrackPtQA[iType]->Fill(var4);
1232
1233       if(tmpTrack==leading){
1234         var4[0] = 0;
1235         fhnTrackPtQA[iType]->Fill(var4);
1236         continue;
1237       }
1238     }
1239     fh1SumPtTrack[iType]->Fill(sumPt);
1240     delete trackIter;
1241   }
1242
1243 }
1244
1245
1246 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1247
1248
1249   // Fill al the matching histos
1250   // It is important that the acceptances for the mathing are as large as possible
1251   // to avoid false matches on the edge of acceptance
1252   // therefore we add some extra matching jets as overhead
1253
1254   static TArrayI aGenIndex(recJetsList.GetEntries());
1255   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1256
1257   static TArrayI aRecIndex(genJetsList.GetEntries());
1258   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1259
1260   if(fDebug){
1261     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1262     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1263   }
1264   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1265                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1266                                             aGenIndex,aRecIndex,fDebug);
1267
1268   if(fDebug){
1269     for(int i = 0;i< aGenIndex.GetSize();++i){
1270       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1271     }
1272     for(int i = 0;i< aRecIndex.GetSize();++i){
1273       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1274     }
1275   }
1276
1277   Double_t container[8];
1278   Double_t containerRec[4];
1279   Double_t containerGen[4];
1280
1281   // loop over generated jets
1282   // consider the
1283   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1284     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1285     Double_t ptGen = genJet->Pt();
1286     Double_t phiGen = genJet->Phi();
1287     if(phiGen<0)phiGen+=TMath::Pi()*2.;
1288     Double_t etaGen = genJet->Eta();
1289     containerGen[0] = ptGen;
1290     containerGen[1] = etaGen;
1291     containerGen[2] = phiGen;
1292     containerGen[3] = genJet->GetPtLeading();
1293
1294     fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1295
1296     Int_t ir = aRecIndex[ig];
1297     if(ir>=0&&ir<recJetsList.GetEntries()){
1298       fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1299
1300       if(JetSelected(genJet)){
1301         fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1302
1303         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1304
1305         fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1306         if(JetSelected(recJet)) {
1307           fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1308
1309           containerGen[3] = recJet->GetPtLeading();
1310           fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner in eta window with leading track on reconstructed level
1311         }
1312       }
1313     }
1314   }// loop over generated jets used for matching...
1315
1316
1317
1318   // loop over reconstructed jets
1319   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1320     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1321     Double_t etaRec = recJet->Eta();
1322     Double_t ptRec = recJet->Pt();
1323     Double_t phiRec = recJet->Phi();
1324     if(phiRec<0)phiRec+=TMath::Pi()*2.;
1325
1326     containerRec[0] = ptRec;
1327     containerRec[1] = etaRec;
1328     containerRec[2] = phiRec;
1329     containerRec[3] = recJet->GetPtLeading();
1330
1331     container[0] = ptRec;
1332     container[1] = etaRec;
1333     container[2] = phiRec;
1334     container[3] = recJet->GetPtLeading();
1335
1336     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1337
1338     if(JetSelected(recJet)){
1339       fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
1340       // Fill Correlation
1341       Int_t ig = aGenIndex[ir];
1342       if(ig>=0 && ig<genJetsList.GetEntries()){
1343         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1344         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1345         Double_t ptGen  = genJet->Pt();
1346         Double_t phiGen = genJet->Phi();
1347         if(phiGen<0)phiGen+=TMath::Pi()*2.;
1348         Double_t etaGen = genJet->Eta();
1349
1350         containerGen[0] = ptGen;
1351         containerGen[1] = etaGen;
1352         containerGen[2] = phiGen;
1353         containerGen[3] = genJet->GetPtLeading();
1354
1355         container[4] = ptGen;
1356         container[5] = etaGen;
1357         container[6] = phiGen;
1358         container[7] = genJet->GetPtLeading();
1359
1360         fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1361
1362         fhnCorrelation->Fill(container);
1363         if(ptGen>0){
1364           Float_t delta = (ptRec-ptGen)/ptGen;
1365           fh2RelPtFGen->Fill(ptGen,delta);
1366           fh2PtFGen->Fill(ptGen,ptRec);
1367         }
1368       }
1369     }// loop over reconstructed jets
1370   }
1371   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1372 }
1373
1374
1375 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1376   //
1377   // Create the particle container for the correction framework manager and
1378   // link it
1379   //
1380   const Int_t kNvar   = 4 ; //number of variables on the grid:pt,eta, phi
1381   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1382   const Double_t kEtamin = -1.0, kEtamax = 1.0;
1383   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1384   const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
1385
1386   // can we neglect migration in eta and phi?
1387   // phi should be no problem since we cover full phi and are phi symmetric
1388   // eta migration is more difficult  due to needed acceptance correction
1389   // in limited eta range
1390
1391   //arrays for the number of bins in each dimension
1392   Int_t iBin[kNvar];
1393   iBin[0] = 125; //bins in pt
1394   iBin[1] =  4; //bins in eta
1395   iBin[2] =  1; // bins in phi
1396   iBin[3] =  10; // bins in leading track Pt
1397
1398   //arrays for lower bounds :
1399   Double_t* binEdges[kNvar];
1400   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1401     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1402
1403   //values for bin lower bounds
1404   //  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);
1405   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1406   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1407   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1408   binEdges[3][0]= kPtLeadingTrackPtMin;
1409   binEdges[3][1]= 1.;
1410   binEdges[3][2]= 2.;
1411   binEdges[3][3]= 3.;
1412   binEdges[3][4]= 4.;
1413   binEdges[3][5]= 5.;
1414   binEdges[3][6]= 6.;
1415   binEdges[3][7]= 8.;
1416   binEdges[3][8]= 10.;
1417   binEdges[3][9]= 12.;
1418   binEdges[3][10]= kPtLeadingTrackPtMax;
1419
1420   fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
1421   for (int k=0; k<kNvar; k++) {
1422     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1423   }
1424
1425   //create correlation matrix for unfolding
1426   Int_t thnDim[2*kNvar];
1427   for (int k=0; k<kNvar; k++) {
1428     //first half  : reconstructed
1429     //second half : MC
1430     thnDim[k]      = iBin[k];
1431     thnDim[k+kNvar] = iBin[k];
1432   }
1433
1434   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
1435   for (int k=0; k<kNvar; k++) {
1436     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1437     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1438   }
1439
1440   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1441     delete [] binEdges[ivar];
1442
1443
1444 }
1445
1446 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1447 {
1448 // Terminate analysis
1449 //
1450     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1451 }
1452
1453
1454 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1455
1456   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1457
1458   Int_t iCount = 0;
1459   if(type==kTrackAOD){
1460     AliAODEvent *aod = 0;
1461     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1462     else aod = AODEvent();
1463     if(!aod){
1464       return iCount;
1465     }
1466     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1467       AliAODTrack *tr = aod->GetTrack(it);
1468       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1469       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1470       if(tr->Pt()<fMinTrackPt)continue;
1471       if(fDebug>0){
1472         if(tr->Pt()>20){
1473           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1474           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1475           tr->Print();
1476           //    tr->Dump();
1477           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1478           if(fESD){
1479             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1480             esdTr->Print("");
1481             // esdTr->Dump();
1482           }
1483         }
1484       }
1485       list->Add(tr);
1486       iCount++;
1487     }
1488   }
1489   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1490     AliMCEvent* mcEvent = MCEvent();
1491     if(!mcEvent)return iCount;
1492     // we want to have alivpartilces so use get track
1493     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1494       if(!mcEvent->IsPhysicalPrimary(it))continue;
1495       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1496       if(part->Pt()<fMinTrackPt)continue;
1497       if(type == kTrackKineAll){
1498         list->Add(part);
1499         iCount++;
1500       }
1501       else if(type == kTrackKineCharged){
1502         if(part->Particle()->GetPDG()->Charge()==0)continue;
1503         list->Add(part);
1504         iCount++;
1505       }
1506     }
1507   }
1508   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1509     AliAODEvent *aod = 0;
1510     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1511     else aod = AODEvent();
1512     if(!aod)return iCount;
1513     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1514     if(!tca)return iCount;
1515     for(int it = 0;it < tca->GetEntriesFast();++it){
1516       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1517       if(!part)continue;
1518       if(part->Pt()<fMinTrackPt)continue;
1519       if(!part->IsPhysicalPrimary())continue;
1520       if(type == kTrackAODMCAll){
1521         list->Add(part);
1522         iCount++;
1523       }
1524       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1525         if(part->Charge()==0)continue;
1526         if(kTrackAODMCCharged){
1527           list->Add(part);
1528         }
1529         else {
1530           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1531           list->Add(part);
1532         }
1533         iCount++;
1534       }
1535     }
1536   }// AODMCparticle
1537   list->Sort();
1538   return iCount;
1539
1540 }
1541
1542
1543 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1544     AliAODEvent *aod = 0;
1545     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1546     else aod = AODEvent();
1547     if(!aod){
1548       return 101;
1549     }
1550     return aod->GetHeader()->GetCentrality();
1551 }
1552
1553
1554
1555 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1556   Bool_t selected = false;
1557
1558   if(!jet)return selected;
1559
1560   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1561     selected = kTRUE;
1562   }
1563   return selected;
1564
1565 }
1566
1567 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1568
1569   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1570   Int_t iCount = 0;
1571
1572   if(!jarray){
1573     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1574     return iCount;
1575   }
1576
1577
1578   for(int ij=0;ij<jarray->GetEntries();ij++){
1579     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1580     if(!jet)continue;
1581     if(type==0){
1582       // No cut at all, main purpose here is sorting
1583       list->Add(jet);
1584       iCount++;
1585     }
1586     else if(type == 1){
1587       // eta cut
1588       if(JetSelected(jet)){
1589         list->Add(jet);
1590         iCount++;
1591       }
1592     }
1593   }
1594
1595   list->Sort();
1596   return iCount;
1597
1598 }
1599
1600
1601 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1602   if(!jets)return 0;
1603
1604   Int_t refMult = 0;
1605   for(int ij = 0;ij < jets->GetEntries();++ij){
1606     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1607     if(!jet)continue;
1608     TRefArray *refs = jet->GetRefTracks();
1609     if(!refs)continue;
1610     refMult += refs->GetEntries();
1611   }
1612   return refMult;
1613
1614 }
1615
1616
1617 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1618   if(!jet)return 0;
1619   TRefArray *refs = jet->GetRefTracks();
1620   if(!refs) return 0;
1621   AliVParticle *leading = 0;
1622   Float_t fMaxPt = 0;
1623   for(int i = 0;i<refs->GetEntries();i++){
1624     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1625     if(!tmp)continue;
1626     if(tmp->Pt()>fMaxPt){
1627       leading = tmp;
1628       fMaxPt = tmp->Pt();
1629     }
1630   }
1631   return leading;
1632 }
1633
1634
1635 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1636   if(!jet)return 0;
1637   if(!list) return 0;
1638   AliVParticle *leading = 0;
1639   Float_t fMaxPt = 0;
1640   for(int i = 0;i<list->GetEntries();i++){
1641     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1642     if(!tmp)continue;
1643     if(jet->DeltaR(tmp)>r)continue;
1644     if(tmp->Pt()>fMaxPt){
1645       leading = tmp;
1646       fMaxPt = tmp->Pt();
1647     }
1648   }
1649   return leading;
1650 }
1651
1652 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1653 {
1654     Int_t phibin=-1;
1655     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1656     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1657     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1658     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1659     return phibin;
1660 }
1661
1662 void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1663   //
1664   // set number of trigger bins
1665   //
1666   if(n>0){
1667       fNTrigger = n;
1668       delete [] fTriggerName;
1669       fTriggerName = new TString [fNTrigger];
1670       delete [] fTriggerBit;fTriggerBit = 0;
1671       fTriggerBit = new UInt_t [fNTrigger];
1672   }
1673   else{
1674     fNTrigger = 0;
1675   }
1676 }
1677
1678
1679 void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1680   //
1681   // set trigger bin
1682   //
1683   if(i<fNTrigger){
1684     fTriggerBit[i] = it;
1685     fTriggerName[i] =  c;
1686   }
1687 }
1688
1689
1690 void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1691   //
1692   // Set number of acceptance bins
1693   //
1694
1695   if(n>0){
1696       fNAcceptance = n;
1697       delete [] fAcceptancePhiMin;
1698       delete [] fAcceptancePhiMax;
1699       delete [] fAcceptanceEtaMin;
1700       delete [] fAcceptanceEtaMax;
1701
1702       fAcceptancePhiMin = new Float_t[fNAcceptance];
1703       fAcceptancePhiMax = new Float_t[fNAcceptance];
1704       fAcceptanceEtaMin = new Float_t[fNAcceptance];
1705       fAcceptanceEtaMax = new Float_t[fNAcceptance];
1706   }
1707   else{
1708     fNTrigger = 0;
1709   }
1710 }
1711
1712 void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1713   //
1714   // Set acceptance bins
1715   //
1716
1717   if(i<fNAcceptance){
1718     Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1719
1720     fAcceptancePhiMin[i] = phiMin;
1721     fAcceptancePhiMax[i] = phiMax;
1722     fAcceptanceEtaMin[i] = etaMin;
1723     fAcceptanceEtaMax[i] = etaMax;
1724   }
1725   else{
1726     fNTrigger = 0;
1727   }
1728 }
1729
1730
1731 Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1732
1733   //
1734   // return aceptnace bin do no use ovelapping bins
1735   //
1736
1737   for(int i = 0;i<fNAcceptance;i++){
1738     if(phi<fAcceptancePhiMin[i])continue;
1739     if(phi>fAcceptancePhiMax[i])continue;
1740     if(eta<fAcceptanceEtaMin[i])continue;
1741     if(eta>fAcceptanceEtaMax[i])continue;
1742     return i;
1743   }
1744   // catch the rest
1745   return fNAcceptance;
1746 }
1747
1748
1749
1750
1751
1752 AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1753   //
1754   delete [] fTriggerBit;
1755   delete [] fTriggerName;
1756
1757   delete [] fAcceptancePhiMin;
1758   delete [] fAcceptancePhiMax;
1759   delete [] fAcceptanceEtaMin;
1760   delete [] fAcceptanceEtaMax;
1761
1762
1763
1764 }