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