PWGJE/AliAnalysisTaskJetSpectrum2.cxx

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