JETAN/AliAnalysisTaskJetCluster.cxx

   1 // **************************************
   2 // Task 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 <TRandom3.h>
  25 #include <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TRefArray.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 <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetCluster.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliJetFinder.h"
  43 #include "AliJetHeader.h"
  44 #include "AliJetReader.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliAODHandler.h"
  48 #include "AliAODTrack.h"
  49 #include "AliAODJet.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliGenPythiaEventHeader.h"
  55 #include "AliJetKineReaderHeader.h"
  56 #include "AliGenCocktailEventHeader.h"
  57 #include "AliInputEventHandler.h"
  58 #include "AliAODJetEventBackground.h"
  59
  60 #include "fastjet/PseudoJet.hh"
  61 #include "fastjet/ClusterSequenceArea.hh"
  62 #include "fastjet/AreaDefinition.hh"
  63 #include "fastjet/JetDefinition.hh"
  64 // get info on how fastjet was configured
  65 #include "fastjet/config.h"
  66
  67
  68 ClassImp(AliAnalysisTaskJetCluster)
  69
  70 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
  71   delete fRef;
  72   delete fRandom;
  73 }
  74
  75 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
  76   fAOD(0x0),
  77   fAODExtension(0x0),
  78   fRef(new TRefArray),
  79   fUseAODTrackInput(kFALSE),
  80   fUseAODMCInput(kFALSE),
  81   fUseGlobalSelection(kFALSE),
  82   fUseBackgroundCalc(kFALSE),
  83   fFilterMask(0),
  84   fTrackTypeRec(kTrackUndef),
  85   fTrackTypeGen(kTrackUndef),
  86   fNSkipLeadingRan(0),
  87   fNRandomCones(10),
  88   fAvgTrials(1),
  89   fExternalWeight(1),
  90   fRecEtaWindow(0.5),
  91   fTrackPtCut(0.),
  92   fJetOutputMinPt(1),
  93   fJetTriggerPtCut(0),
  94   fCentCutUp(0),
  95   fCentCutLo(0),
  96   fNonStdBranch(""),
  97   fBackgroundBranch(""),
  98   fNonStdFile(""),
  99   fRparam(1.0),
 100   fAlgorithm(fastjet::kt_algorithm),
 101   fStrategy(fastjet::Best),
 102   fRecombScheme(fastjet::BIpt_scheme),
 103   fAreaType(fastjet::active_area),
 104   fGhostArea(0.01),
 105   fActiveAreaRepeats(1),
 106   fGhostEtamax(1.5),
 107   fRandom(0),
 108   fh1Xsec(0x0),
 109   fh1Trials(0x0),
 110   fh1PtHard(0x0),
 111   fh1PtHardNoW(0x0),
 112   fh1PtHardTrials(0x0),
 113   fh1NJetsRec(0x0),
 114   fh1NConstRec(0x0),
 115   fh1NConstLeadingRec(0x0),
 116   fh1PtJetsRecIn(0x0),
 117   fh1PtJetsLeadingRecIn(0x0),
 118   fh1PtJetConstRec(0x0),
 119   fh1PtJetConstLeadingRec(0x0),
 120   fh1PtTracksRecIn(0x0),
 121   fh1PtTracksLeadingRecIn(0x0),
 122   fh1NJetsRecRan(0x0),
 123   fh1NConstRecRan(0x0),
 124   fh1PtJetsLeadingRecInRan(0x0),
 125   fh1NConstLeadingRecRan(0x0),
 126   fh1PtJetsRecInRan(0x0),
 127   fh1PtTracksGenIn(0x0),
 128   fh1Nch(0x0),
 129   fh1CentralityPhySel(0x0),
 130   fh1Centrality(0x0),
 131   fh1CentralitySelect(0x0),
 132   fh1ZPhySel(0x0),
 133   fh1Z(0x0),
 134   fh1ZSelect(0x0),
 135   fh2NRecJetsPt(0x0),
 136   fh2NRecTracksPt(0x0),
 137   fh2NConstPt(0x0),
 138   fh2NConstLeadingPt(0x0),
 139   fh2JetPhiEta(0x0),
 140   fh2LeadingJetPhiEta(0x0),
 141   fh2JetEtaPt(0x0),
 142   fh2LeadingJetEtaPt(0x0),
 143   fh2TrackEtaPt(0x0),
 144   fh2LeadingTrackEtaPt(0x0),
 145   fh2JetsLeadingPhiEta(0x0),
 146   fh2JetsLeadingPhiPt(0x0),
 147   fh2TracksLeadingPhiEta(0x0),
 148   fh2TracksLeadingPhiPt(0x0),
 149   fh2TracksLeadingJetPhiPt(0x0),
 150   fh2JetsLeadingPhiPtW(0x0),
 151   fh2TracksLeadingPhiPtW(0x0),
 152   fh2TracksLeadingJetPhiPtW(0x0),
 153   fh2NRecJetsPtRan(0x0),
 154   fh2NConstPtRan(0x0),
 155   fh2NConstLeadingPtRan(0x0),
 156   fh2PtNch(0x0),
 157   fh2PtNchRan(0x0),
 158   fh2PtNchN(0x0),
 159   fh2PtNchNRan(0x0),
 160   fh2TracksLeadingJetPhiPtRan(0x0),
 161   fh2TracksLeadingJetPhiPtWRan(0x0),
 162   fHistList(0x0)
 163 {
 164   for(int i = 0;i<3;i++){
 165     fh1BiARandomCones[i] = 0;
 166     fh1BiARandomConesRan[i] = 0;
 167   }
 168   for(int i = 0;i<kMaxCent;i++){
 169     fh2JetsLeadingPhiPtC[i] = 0;
 170     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 171     fh2TracksLeadingJetPhiPtC[i] = 0;
 172     fh2TracksLeadingJetPhiPtWC[i] = 0;
 173   }
 174 }
 175
 176 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
 177   AliAnalysisTaskSE(name),
 178   fAOD(0x0),
 179   fAODExtension(0x0),
 180   fRef(new TRefArray),
 181   fUseAODTrackInput(kFALSE),
 182   fUseAODMCInput(kFALSE),
 183   fUseGlobalSelection(kFALSE),
 184   fUseBackgroundCalc(kFALSE),
 185   fFilterMask(0),
 186   fTrackTypeRec(kTrackUndef),
 187   fTrackTypeGen(kTrackUndef),
 188   fNSkipLeadingRan(0),
 189   fNRandomCones(5),
 190   fAvgTrials(1),
 191   fExternalWeight(1),
 192   fRecEtaWindow(0.5),
 193   fTrackPtCut(0.),
 194   fJetOutputMinPt(1),
 195   fJetTriggerPtCut(0),
 196   fCentCutUp(0),
 197   fCentCutLo(0),
 198   fNonStdBranch(""),
 199   fBackgroundBranch(""),
 200   fNonStdFile(""),
 201   fRparam(1.0),
 202   fAlgorithm(fastjet::kt_algorithm),
 203   fStrategy(fastjet::Best),
 204   fRecombScheme(fastjet::BIpt_scheme),
 205   fAreaType(fastjet::active_area),
 206   fGhostArea(0.01),
 207   fActiveAreaRepeats(1),
 208   fGhostEtamax(1.5),
 209   fRandom(0),
 210   fh1Xsec(0x0),
 211   fh1Trials(0x0),
 212   fh1PtHard(0x0),
 213   fh1PtHardNoW(0x0),
 214   fh1PtHardTrials(0x0),
 215   fh1NJetsRec(0x0),
 216   fh1NConstRec(0x0),
 217   fh1NConstLeadingRec(0x0),
 218   fh1PtJetsRecIn(0x0),
 219   fh1PtJetsLeadingRecIn(0x0),
 220   fh1PtJetConstRec(0x0),
 221   fh1PtJetConstLeadingRec(0x0),
 222   fh1PtTracksRecIn(0x0),
 223   fh1PtTracksLeadingRecIn(0x0),
 224   fh1NJetsRecRan(0x0),
 225   fh1NConstRecRan(0x0),
 226   fh1PtJetsLeadingRecInRan(0x0),
 227   fh1NConstLeadingRecRan(0x0),
 228   fh1PtJetsRecInRan(0x0),
 229   fh1PtTracksGenIn(0x0),
 230   fh1Nch(0x0),
 231   fh1CentralityPhySel(0x0),
 232   fh1Centrality(0x0),
 233   fh1CentralitySelect(0x0),
 234   fh1ZPhySel(0x0),
 235   fh1Z(0x0),
 236   fh1ZSelect(0x0),
 237   fh2NRecJetsPt(0x0),
 238   fh2NRecTracksPt(0x0),
 239   fh2NConstPt(0x0),
 240   fh2NConstLeadingPt(0x0),
 241   fh2JetPhiEta(0x0),
 242   fh2LeadingJetPhiEta(0x0),
 243   fh2JetEtaPt(0x0),
 244   fh2LeadingJetEtaPt(0x0),
 245   fh2TrackEtaPt(0x0),
 246   fh2LeadingTrackEtaPt(0x0),
 247   fh2JetsLeadingPhiEta(0x0),
 248   fh2JetsLeadingPhiPt(0x0),
 249   fh2TracksLeadingPhiEta(0x0),
 250   fh2TracksLeadingPhiPt(0x0),
 251   fh2TracksLeadingJetPhiPt(0x0),
 252   fh2JetsLeadingPhiPtW(0x0),
 253   fh2TracksLeadingPhiPtW(0x0),
 254   fh2TracksLeadingJetPhiPtW(0x0),
 255   fh2NRecJetsPtRan(0x0),
 256   fh2NConstPtRan(0x0),
 257   fh2NConstLeadingPtRan(0x0),
 258   fh2PtNch(0x0),
 259   fh2PtNchRan(0x0),
 260   fh2PtNchN(0x0),
 261   fh2PtNchNRan(0x0),
 262   fh2TracksLeadingJetPhiPtRan(0x0),
 263   fh2TracksLeadingJetPhiPtWRan(0x0),
 264   fHistList(0x0)
 265 {
 266   for(int i = 0;i<3;i++){
 267     fh1BiARandomCones[i] = 0;
 268     fh1BiARandomConesRan[i] = 0;
 269   }
 270   for(int i = 0;i<kMaxCent;i++){
 271     fh2JetsLeadingPhiPtC[i] = 0;
 272     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 273     fh2TracksLeadingJetPhiPtC[i] = 0;
 274     fh2TracksLeadingJetPhiPtWC[i] = 0;
 275   }
 276  DefineOutput(1, TList::Class());
 277 }
 278
 279
 280
 281 Bool_t AliAnalysisTaskJetCluster::Notify()
 282 {
 283   //
 284   // Implemented Notify() to read the cross sections
 285   // and number of trials from pyxsec.root
 286   //
 287   return kTRUE;
 288 }
 289
 290 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 291 {
 292
 293   //
 294   // Create the output container
 295   //
 296
 297   fRandom = new TRandom3(0);
 298
 299
 300   // Connect the AOD
 301
 302
 303   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 304
 305
 306
 307   if(fNonStdBranch.Length()!=0)
 308     {
 309       // only create the output branch if we have a name
 310       // Create a new branch for jets...
 311       //  -> cleared in the UserExec....
 312       // here we can also have the case that the brnaches are written to a separate file
 313
 314       TClonesArray *tca = new TClonesArray("AliAODJet", 0);
 315       tca->SetName(fNonStdBranch.Data());
 316       AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
 317
 318
 319       TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
 320       tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
 321       AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
 322
 323
 324       if(fUseBackgroundCalc){
 325         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
 326           AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
 327           evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
 328           AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
 329         }
 330         // create the branch for the random cones with the same R
 331         TString cName = Form("%sRandomCone",fNonStdBranch.Data());
 332         if(!AODEvent()->FindListObject(cName.Data())){
 333           TClonesArray *tcaC = new TClonesArray("AliAODJet", 0);
 334           tcaC->SetName(cName.Data());
 335           AddAODBranch("TClonesArray",&tcaC,fNonStdFile.Data());
 336         }
 337
 338         // create the branch with the random for the random cones on the random event
 339         cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 340         if(!AODEvent()->FindListObject(cName.Data())){
 341           TClonesArray *tcaCran = new TClonesArray("AliAODJet", 0);
 342           tcaCran->SetName(cName.Data());
 343           AddAODBranch("TClonesArray",&tcaCran,fNonStdFile.Data());
 344         }
 345       }
 346
 347       if(fNonStdFile.Length()!=0){
 348         //
 349         // case that we have an AOD extension we need to fetch the jets from the extended output
 350         // we identifay the extension aod event by looking for the branchname
 351         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 352         TObjArray* extArray = aodH->GetExtensions();
 353         if (extArray) {
 354           TIter next(extArray);
 355           while ((fAODExtension=(AliAODExtension*)next())){
 356             TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
 357             if(fDebug>10){
 358               Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
 359               fAODExtension->GetAOD()->Dump();
 360             }
 361             if(obj){
 362               if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
 363               break;
 364             }
 365             fAODExtension = 0;
 366           }
 367         }
 368       }
 369     }
 370
 371
 372   if(!fHistList)fHistList = new TList();
 373   fHistList->SetOwner();
 374
 375   Bool_t oldStatus = TH1::AddDirectoryStatus();
 376   TH1::AddDirectory(kFALSE);
 377
 378   //
 379   //  Histogram
 380
 381   const Int_t nBinPt = 200;
 382   Double_t binLimitsPt[nBinPt+1];
 383   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 384     if(iPt == 0){
 385       binLimitsPt[iPt] = 0.0;
 386     }
 387     else {// 1.0
 388       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 389     }
 390   }
 391
 392   const Int_t nBinPhi = 90;
 393   Double_t binLimitsPhi[nBinPhi+1];
 394   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 395     if(iPhi==0){
 396       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 397     }
 398     else{
 399       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 400     }
 401   }
 402
 403
 404
 405   const Int_t nBinEta = 40;
 406   Double_t binLimitsEta[nBinEta+1];
 407   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 408     if(iEta==0){
 409       binLimitsEta[iEta] = -2.0;
 410     }
 411     else{
 412       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 413     }
 414   }
 415
 416   const int nChMax = 4000;
 417
 418   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 419   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 420
 421   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 422   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 423
 424
 425   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 426   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 427
 428   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 429   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 430   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 431   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 432
 433
 434   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 435   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 436   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 437
 438   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 439   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 440   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 441   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 442   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 443   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 444   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 445   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 446   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 447   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 448
 449   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 450   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
 451   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 452
 453   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 454   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
 455   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 456
 457   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 458   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 459   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 460   //
 461
 462
 463   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 464   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 465   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 466   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 467
 468   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 469   fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 470   fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 471   fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 472
 473
 474
 475   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 476                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 477   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 478                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 479
 480   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 481                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 482   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 483                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 484
 485   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 486                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 487   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 488                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 489
 490
 491
 492   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 493                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 494   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 495                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 496   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 497                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 498   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 499                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 500   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 501                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 502   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 503                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 504
 505   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 506                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 507   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 508                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 509
 510   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 511                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 512   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 513                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 514
 515
 516   if(fUseBackgroundCalc){
 517     for(int i = 0;i<3;i++){
 518       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 519       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 520     }
 521   }
 522   for(int i = 0;i < kMaxCent;i++){
 523     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
 524     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
 525     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
 526     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
 527   }
 528
 529   const Int_t saveLevel = 3; // large save level more histos
 530   if(saveLevel>0){
 531     fHistList->Add(fh1Xsec);
 532     fHistList->Add(fh1Trials);
 533
 534     fHistList->Add(fh1NJetsRec);
 535     fHistList->Add(fh1NConstRec);
 536     fHistList->Add(fh1NConstLeadingRec);
 537     fHistList->Add(fh1PtJetsRecIn);
 538     fHistList->Add(fh1PtJetsLeadingRecIn);
 539     fHistList->Add(fh1PtTracksRecIn);
 540     fHistList->Add(fh1PtTracksLeadingRecIn);
 541     fHistList->Add(fh1PtJetConstRec);
 542     fHistList->Add(fh1PtJetConstLeadingRec);
 543     fHistList->Add(fh1NJetsRecRan);
 544     fHistList->Add(fh1NConstRecRan);
 545     fHistList->Add(fh1PtJetsLeadingRecInRan);
 546     fHistList->Add(fh1NConstLeadingRecRan);
 547     fHistList->Add(fh1PtJetsRecInRan);
 548     fHistList->Add(fh1Nch);
 549     fHistList->Add(fh1Centrality);
 550     fHistList->Add(fh1CentralitySelect);
 551     fHistList->Add(fh1CentralityPhySel);
 552     fHistList->Add(fh1Z);
 553     fHistList->Add(fh1ZSelect);
 554     fHistList->Add(fh1ZPhySel);
 555     if(fUseBackgroundCalc){
 556       for(int i = 0;i<3;i++){
 557         fHistList->Add(fh1BiARandomCones[i]);
 558         fHistList->Add(fh1BiARandomConesRan[i]);
 559       }
 560     }
 561   for(int i = 0;i < kMaxCent;i++){
 562     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
 563     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
 564     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
 565     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
 566   }
 567
 568     fHistList->Add(fh2NRecJetsPt);
 569     fHistList->Add(fh2NRecTracksPt);
 570     fHistList->Add(fh2NConstPt);
 571     fHistList->Add(fh2NConstLeadingPt);
 572     fHistList->Add(fh2PtNch);
 573     fHistList->Add(fh2PtNchRan);
 574     fHistList->Add(fh2PtNchN);
 575     fHistList->Add(fh2PtNchNRan);
 576     fHistList->Add(fh2JetPhiEta);
 577     fHistList->Add(fh2LeadingJetPhiEta);
 578     fHistList->Add(fh2JetEtaPt);
 579     fHistList->Add(fh2LeadingJetEtaPt);
 580     fHistList->Add(fh2TrackEtaPt);
 581     fHistList->Add(fh2LeadingTrackEtaPt);
 582     fHistList->Add(fh2JetsLeadingPhiEta );
 583     fHistList->Add(fh2JetsLeadingPhiPt);
 584     fHistList->Add(fh2TracksLeadingPhiEta);
 585     fHistList->Add(fh2TracksLeadingPhiPt);
 586     fHistList->Add(fh2TracksLeadingJetPhiPt);
 587     fHistList->Add(fh2JetsLeadingPhiPtW);
 588     fHistList->Add(fh2TracksLeadingPhiPtW);
 589     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 590     fHistList->Add(fh2NRecJetsPtRan);
 591     fHistList->Add(fh2NConstPtRan);
 592     fHistList->Add(fh2NConstLeadingPtRan);
 593     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 594     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 595     }
 596
 597   // =========== Switch on Sumw2 for all histos ===========
 598   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 599     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 600     if (h1){
 601       h1->Sumw2();
 602       continue;
 603     }
 604     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 605     if(hn)hn->Sumw2();
 606   }
 607   TH1::AddDirectory(oldStatus);
 608 }
 609
 610 void AliAnalysisTaskJetCluster::Init()
 611 {
 612   //
 613   // Initialization
 614   //
 615
 616   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 617
 618 }
 619
 620 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 621 {
 622
 623   if(fUseGlobalSelection){
 624     // no selection by the service task, we continue
 625     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 626     PostData(1, fHistList);
 627     return;
 628   }
 629
 630
 631
 632   // handle and reset the output jet branch
 633   // only need this once
 634   TClonesArray* jarray = 0;
 635   TClonesArray* jarrayran = 0;
 636   TClonesArray* rConeArray = 0;
 637   TClonesArray* rConeArrayRan = 0;
 638   AliAODJetEventBackground*  evBkg = 0;
 639   if(fNonStdBranch.Length()!=0) {
 640     if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
 641     if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
 642     if(jarray)jarray->Delete();    // this is our responsibility, clear before filling again
 643     if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
 644     if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
 645     if(jarrayran)jarrayran->Delete();    // this is our responsibility, clear before filling again
 646
 647     if(fUseBackgroundCalc){
 648       if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
 649       if(!evBkg)  evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
 650       if(evBkg)evBkg->Reset();
 651
 652       TString cName = Form("%sRandomCone",fNonStdBranch.Data());
 653       if(AODEvent())rConeArray = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
 654       if(!rConeArray)rConeArray =  (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
 655       if(rConeArray)rConeArray->Delete();
 656
 657       cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 658       if(AODEvent())rConeArrayRan = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
 659       if(!rConeArrayRan)rConeArrayRan =  (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
 660       if(rConeArrayRan)rConeArrayRan->Delete();
 661     }
 662   }
 663
 664   static AliAODJetEventBackground* externalBackground = 0;
 665   if(!externalBackground&&fBackgroundBranch.Length()){
 666     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
 667     Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 668   }
 669   //
 670   // Execute analysis for current event
 671   //
 672   AliESDEvent *fESD = 0;
 673   if(fUseAODTrackInput){
 674     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 675     if(!fAOD){
 676       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 677       return;
 678     }
 679     // fethc the header
 680   }
 681   else{
 682     //  assume that the AOD is in the general output...
 683     fAOD  = AODEvent();
 684     if(!fAOD){
 685       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 686       return;
 687     }
 688     if(fDebug>0){
 689       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 690     }
 691   }
 692
 693   Bool_t selectEvent =  false;
 694   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 695
 696   Float_t cent = 0;
 697   Float_t zVtx  = 0;
 698   Int_t cenClass = -1;
 699   if(fAOD){
 700     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 701     TString vtxTitle(vtxAOD->GetTitle());
 702     zVtx = vtxAOD->GetZ();
 703
 704     cent = fAOD->GetHeader()->GetCentrality();
 705     if(cent<10)cenClass = 0;
 706     else if(cent<30)cenClass = 1;
 707     else if(cent<50)cenClass = 2;
 708     else if(cent<80)cenClass = 3;
 709     if(physicsSelection){
 710       fh1CentralityPhySel->Fill(cent);
 711       fh1ZPhySel->Fill(zVtx);
 712     }
 713
 714
 715     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 716         Float_t yvtx = vtxAOD->GetY();
 717         Float_t xvtx = vtxAOD->GetX();
 718         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 719         if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
 720           if(physicsSelection){
 721             selectEvent = true;
 722           }
 723         }
 724     }
 725     if(fCentCutUp>0){
 726       if(cent<fCentCutLo||cent>fCentCutUp){
 727         selectEvent = false;
 728       }
 729     }
 730
 731   }
 732   if(!selectEvent){
 733     PostData(1, fHistList);
 734     return;
 735   }
 736   fh1Centrality->Fill(cent);
 737   fh1Z->Fill(zVtx);
 738   fh1Trials->Fill("#sum{ntrials}",1);
 739
 740
 741   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 742
 743   // ==== General variables needed
 744
 745
 746
 747   // we simply fetch the tracks/mc particles as a list of AliVParticles
 748
 749   TList recParticles;
 750   TList genParticles;
 751
 752   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 753   Float_t nCh = recParticles.GetEntries();
 754   fh1Nch->Fill(nCh);
 755   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 756   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 757   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 758
 759   // find the jets....
 760
 761   vector<fastjet::PseudoJet> inputParticlesRec;
 762   vector<fastjet::PseudoJet> inputParticlesRecRan;
 763
 764   // Generate the random cones
 765
 766   AliAODJet vTmpRan(1,0,0,1);
 767   for(int i = 0; i < recParticles.GetEntries(); i++){
 768     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 769     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 770     // we take total momentum here
 771     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 772     jInp.set_user_index(i);
 773     inputParticlesRec.push_back(jInp);
 774
 775     // the randomized input changes eta and phi, but keeps the p_T
 776     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 777       Double_t pT = vp->Pt();
 778       Double_t eta = 1.8 * fRandom->Rndm() - 0.9;
 779       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
 780
 781       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
 782       Double_t pZ = pT/TMath::Tan(theta);
 783
 784       Double_t pX = pT * TMath::Cos(phi);
 785       Double_t pY = pT * TMath::Sin(phi);
 786       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 787       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 788
 789       jInpRan.set_user_index(i);
 790       inputParticlesRecRan.push_back(jInpRan);
 791       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
 792     }
 793
 794     // fill the tref array, only needed when we write out jets
 795     if(jarray){
 796       if(i == 0){
 797         fRef->Delete(); // make sure to delete before placement new...
 798         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
 799       }
 800       fRef->Add(vp);
 801     }
 802   }// recparticles
 803
 804   if(inputParticlesRec.size()==0){
 805     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 806     PostData(1, fHistList);
 807     return;
 808   }
 809
 810   // run fast jet
 811   // employ setters for these...
 812
 813
 814   // now create the object that holds info about ghosts
 815   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
 816     // reduce CPU time...
 817     fGhostArea = 0.5;
 818     fActiveAreaRepeats = 0;
 819   }
 820
 821  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 822   fastjet::AreaType areaType =   fastjet::active_area;
 823   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 824   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 825   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 826   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
 827
 828   //range where to compute background
 829   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 830   phiMin = 0;
 831   phiMax = 2*TMath::Pi();
 832   rapMax = fGhostEtamax - fRparam;
 833   rapMin = - fGhostEtamax + fRparam;
 834   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 835
 836
 837
 838   inclusiveJets = clustSeq.inclusive_jets();
 839   sortedJets    = sorted_by_pt(inclusiveJets);
 840
 841   fh1NJetsRec->Fill(sortedJets.size());
 842
 843  // loop over all jets an fill information, first on is the leading jet
 844
 845   Int_t nRecOver = inclusiveJets.size();
 846   Int_t nRec     = inclusiveJets.size();
 847   if(inclusiveJets.size()>0){
 848     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 849     Double_t area = clustSeq.area(sortedJets[0]);
 850     leadingJet.SetEffArea(area,0);
 851     Float_t pt = leadingJet.Pt();
 852     Int_t nAodOutJets = 0;
 853     Int_t nAodOutTracks = 0;
 854     AliAODJet *aodOutJet = 0;
 855
 856     Int_t iCount = 0;
 857     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 858       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 859       while(pt<ptCut&&iCount<nRec){
 860         nRecOver--;
 861         iCount++;
 862         if(iCount<nRec){
 863           pt = sortedJets[iCount].perp();
 864         }
 865       }
 866       if(nRecOver<=0)break;
 867       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 868     }
 869     Float_t phi = leadingJet.Phi();
 870     if(phi<0)phi+=TMath::Pi()*2.;
 871     Float_t eta = leadingJet.Eta();
 872     Float_t pTback = 0;
 873     if(externalBackground){
 874       // carefull has to be filled in a task before
 875       // todo, ReArrange to the botom
 876       pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
 877     }
 878     pt = leadingJet.Pt() - pTback;
 879     // correlation of leading jet with tracks
 880     TIterator *recIter = recParticles.MakeIterator();
 881     recIter->Reset();
 882     AliVParticle *tmpRecTrack = 0;
 883     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 884       Float_t tmpPt = tmpRecTrack->Pt();
 885       // correlation
 886       Float_t tmpPhi =  tmpRecTrack->Phi();
 887       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 888       Float_t dPhi = phi - tmpPhi;
 889       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 890       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 891       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 892       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 893       if(cenClass>=0){
 894         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
 895         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
 896       }
 897
 898     }
 899
 900
 901
 902     for(int j = 0; j < nRec;j++){
 903       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 904       aodOutJet = 0;
 905       nAodOutTracks = 0;
 906       Float_t tmpPt = tmpRec.Pt();
 907       Float_t tmpPtBack = 0;
 908       if(externalBackground){
 909         // carefull has to be filled in a task before
 910        // todo, ReArrange to the botom
 911         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
 912       }
 913       tmpPt = tmpPt - tmpPtBack;
 914       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
 915
 916       fh1PtJetsRecIn->Fill(tmpPt);
 917       // Fill Spectra with constituents
 918       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 919       fh1NConstRec->Fill(constituents.size());
 920       fh2PtNch->Fill(nCh,tmpPt);
 921       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 922       fh2NConstPt->Fill(tmpPt,constituents.size());
 923       // loop over constiutents and fill spectrum
 924
 925       // Add the jet information and the track references to the output AOD
 926
 927       if(tmpPt>fJetOutputMinPt&&jarray){
 928         aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
 929         Double_t area1 = clustSeq.area(sortedJets[j]);
 930         aodOutJet->SetEffArea(area1,0);
 931       }
 932
 933
 934       if(fUseBackgroundCalc){
 935
 936         // create a random jet within the acceptance
 937         Double_t etaMax = 0.9 - fRparam;
 938         Int_t nCone = 0;
 939         Int_t nConeRan = 0;
 940         Double_t pTC = 1; // small number
 941         for(int ir = 0;ir < fNRandomCones;ir++){
 942           Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
 943           Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
 944           // massless jet
 945           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 946           Double_t pZC = pTC/TMath::Tan(thetaC);
 947           Double_t pXC = pTC * TMath::Cos(phiC);
 948           Double_t pYC = pTC * TMath::Sin(phiC);
 949           Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 950           AliAODJet tmpRecC (pXC,pYC,pZC, pC);
 951           bool skip = false;
 952           for(int jj = 0; jj < TMath::Min(nRec,2);jj++){
 953             AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
 954             if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
 955               skip = true;
 956               break;
 957             }
 958           }
 959           if(skip)continue;
 960           tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
 961           if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
 962           if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
 963         }// random cones
 964
 965
 966         // loop over the reconstructed particles and add up the pT in the random cones
 967         // maybe better to loop over randomized particles not in the real jets...
 968         // but this by definition brings dow average energy in the whole  event
 969         AliAODJet vTmpRanR(1,0,0,1);
 970         for(int i = 0; i < recParticles.GetEntries(); i++){
 971           AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 972           if(rConeArray){
 973             for(int ir = 0;ir < fNRandomCones;ir++){
 974               AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
 975               if(jC&&jC->DeltaR(vp)<fRparam){
 976                 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
 977               }
 978             }
 979           }// add up energy in cone
 980
 981           // the randomized input changes eta and phi, but keeps the p_T
 982           if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 983             Double_t pTR = vp->Pt();
 984             Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
 985             Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
 986
 987             Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
 988             Double_t pZR = pTR/TMath::Tan(thetaR);
 989
 990             Double_t pXR = pTR * TMath::Cos(phiR);
 991             Double_t pYR = pTR * TMath::Sin(phiR);
 992             Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
 993             vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
 994             if(rConeArrayRan){
 995               for(int ir = 0;ir < fNRandomCones;ir++){
 996                 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
 997                 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
 998                   jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
 999                 }
1000               }
1001             }
1002           }
1003         }// loop over recparticles
1004
1005         Float_t jetArea = fRparam*fRparam*TMath::Pi();
1006         for(int ir = 0;ir < fNRandomCones;ir++){
1007           // rescale the momntum vectors for the random cones
1008           if(!rConeArray)continue;
1009           AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
1010           if(rC){
1011             Double_t etaC = rC->Eta();
1012             Double_t phiC = rC->Phi();
1013             // massless jet, unit vector
1014             pTC = rC->ChargedBgEnergy();
1015             if(pTC<=0)pTC = 0.1; // for almost empty events
1016             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1017             Double_t pZC = pTC/TMath::Tan(thetaC);
1018             Double_t pXC = pTC * TMath::Cos(phiC);
1019             Double_t pYC = pTC * TMath::Sin(phiC);
1020             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1021             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1022             rC->SetBgEnergy(0,0);
1023             rC->SetEffArea(jetArea,0);
1024           }
1025           rC = (AliAODJet*)rConeArrayRan->At(ir);
1026           // same wit random
1027           if(rC){
1028             Double_t etaC = rC->Eta();
1029             Double_t phiC = rC->Phi();
1030             // massless jet, unit vector
1031             pTC = rC->ChargedBgEnergy();
1032             if(pTC<=0)pTC = 0.1;// for almost empty events
1033             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1034             Double_t pZC = pTC/TMath::Tan(thetaC);
1035             Double_t pXC = pTC * TMath::Cos(phiC);
1036             Double_t pYC = pTC * TMath::Sin(phiC);
1037             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1038             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1039             rC->SetBgEnergy(0,0);
1040             rC->SetEffArea(jetArea,0);
1041           }
1042         }
1043       }// if(fUseBackgroundCalc
1044
1045      for(unsigned int ic = 0; ic < constituents.size();ic++){
1046        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1047        fh1PtJetConstRec->Fill(part->Pt());
1048        if(aodOutJet){
1049          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1050        }
1051        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1052      }
1053
1054      // correlation
1055      Float_t tmpPhi =  tmpRec.Phi();
1056      Float_t tmpEta =  tmpRec.Eta();
1057      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1058
1059
1060
1061      if(j==0){
1062        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1063        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1064        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1065        fh1NConstLeadingRec->Fill(constituents.size());
1066        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1067        continue;
1068      }
1069      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1070      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1071      Float_t dPhi = phi - tmpPhi;
1072      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1073      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1074      Float_t dEta = eta - tmpRec.Eta();
1075      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1076      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1077      if(cenClass>=0){
1078        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1079        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1080      }
1081      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1082     }
1083    delete recIter;
1084
1085    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1086
1087
1088    if(evBkg){
1089      vector<fastjet::PseudoJet> jets2=sortedJets;
1090      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1091      Double_t bkg1=0;
1092      Double_t sigma1=0.;
1093      Double_t meanarea1=0.;
1094      Double_t bkg2=0;
1095      Double_t sigma2=0.;
1096      Double_t meanarea2=0.;
1097
1098      clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
1099      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
1100
1101      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1102      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1103
1104      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1105      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
1106      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1107      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1108
1109    }
1110   }
1111
1112
1113
1114
1115
1116   // fill track information
1117   Int_t nTrackOver = recParticles.GetSize();
1118   // do the same for tracks and jets
1119
1120   if(nTrackOver>0){
1121    TIterator *recIter = recParticles.MakeIterator();
1122    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1123    Float_t pt = tmpRec->Pt();
1124
1125    //    Printf("Leading track p_t %3.3E",pt);
1126    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1127      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1128      while(pt<ptCut&&tmpRec){
1129        nTrackOver--;
1130        tmpRec = (AliVParticle*)(recIter->Next());
1131        if(tmpRec){
1132          pt = tmpRec->Pt();
1133        }
1134      }
1135      if(nTrackOver<=0)break;
1136      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1137    }
1138
1139    recIter->Reset();
1140    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1141    Float_t phi = leading->Phi();
1142    if(phi<0)phi+=TMath::Pi()*2.;
1143    Float_t eta = leading->Eta();
1144    pt  = leading->Pt();
1145    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1146      Float_t tmpPt = tmpRec->Pt();
1147      Float_t tmpEta = tmpRec->Eta();
1148      fh1PtTracksRecIn->Fill(tmpPt);
1149      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1150      if(tmpRec==leading){
1151        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1152        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1153        continue;
1154      }
1155       // correlation
1156      Float_t tmpPhi =  tmpRec->Phi();
1157
1158      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1159      Float_t dPhi = phi - tmpPhi;
1160      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1161      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1162      Float_t dEta = eta - tmpRec->Eta();
1163      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1164      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1165      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1166    }
1167    delete recIter;
1168  }
1169
1170  // find the random jets
1171  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1172  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1173
1174  inclusiveJetsRan = clustSeqRan.inclusive_jets();
1175  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1176
1177  // fill the jet information from random track
1178
1179   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1180
1181  // loop over all jets an fill information, first on is the leading jet
1182
1183  Int_t nRecOverRan = inclusiveJetsRan.size();
1184  Int_t nRecRan     = inclusiveJetsRan.size();
1185  if(inclusiveJetsRan.size()>0){
1186    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1187    Float_t pt = leadingJet.Pt();
1188
1189    Int_t iCount = 0;
1190    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1191      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1192       while(pt<ptCut&&iCount<nRecRan){
1193         nRecOverRan--;
1194         iCount++;
1195         if(iCount<nRecRan){
1196           pt = sortedJetsRan[iCount].perp();
1197         }
1198       }
1199       if(nRecOverRan<=0)break;
1200       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1201     }
1202     Float_t phi = leadingJet.Phi();
1203     if(phi<0)phi+=TMath::Pi()*2.;
1204     pt = leadingJet.Pt();
1205
1206     // correlation of leading jet with random tracks
1207
1208     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1209       {
1210         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1211         // correlation
1212         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1213         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1214         Float_t dPhi = phi - tmpPhi;
1215         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1216         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1217         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1218         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1219       }
1220
1221     Int_t nAodOutJetsRan = 0;
1222      AliAODJet *aodOutJetRan = 0;
1223     for(int j = 0; j < nRecRan;j++){
1224       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1225       Float_t tmpPt = tmpRec.Pt();
1226       fh1PtJetsRecInRan->Fill(tmpPt);
1227       // Fill Spectra with constituents
1228       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1229       fh1NConstRecRan->Fill(constituents.size());
1230       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1231       fh2PtNchRan->Fill(nCh,tmpPt);
1232       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1233
1234
1235      if(tmpPt>fJetOutputMinPt&&jarrayran){
1236        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1237        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1238
1239        aodOutJetRan->SetEffArea(arearan,0);    }
1240
1241
1242
1243
1244      for(unsigned int ic = 0; ic < constituents.size();ic++){
1245        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1246        // fh1PtJetConstRec->Fill(part->Pt());
1247        if(aodOutJetRan){
1248          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1249        }
1250      }
1251
1252
1253
1254
1255       // correlation
1256       Float_t tmpPhi =  tmpRec.Phi();
1257       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1258
1259       if(j==0){
1260         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1261         fh1NConstLeadingRecRan->Fill(constituents.size());
1262         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1263         continue;
1264       }
1265     }
1266
1267
1268     if(evBkg){
1269      Double_t bkg3=0.;
1270      Double_t sigma3=0.;
1271      Double_t meanarea3=0.;
1272      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1273      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1274      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1275      /*
1276      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1277      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1278      */
1279     }
1280
1281
1282
1283  }
1284
1285
1286  // do the event selection if activated
1287  if(fJetTriggerPtCut>0){
1288    bool select = false;
1289    Float_t minPt = fJetTriggerPtCut;
1290    /*
1291    // hard coded for now ...
1292    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1293    if(cent<10)minPt = 50;
1294    else if(cent<30)minPt = 42;
1295    else if(cent<50)minPt = 28;
1296    else if(cent<80)minPt = 18;
1297    */
1298    float rho = 0;
1299    if(externalBackground)rho = externalBackground->GetBackground(2);
1300    if(jarray){
1301      for(int i = 0;i < jarray->GetEntriesFast();i++){
1302        AliAODJet *jet = (AliAODJet*)jarray->At(i);
1303        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1304        if(ptSub>=minPt){
1305          select = true;
1306          break;
1307        }
1308      }
1309    }
1310
1311    if(select){
1312      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1313      fh1CentralitySelect->Fill(cent);
1314      fh1ZSelect->Fill(zVtx);
1315      aodH->SetFillAOD(kTRUE);
1316    }
1317  }
1318
1319  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1320  PostData(1, fHistList);
1321 }
1322
1323 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1324 {
1325 // Terminate analysis
1326 //
1327     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1328 }
1329
1330
1331 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1332
1333   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1334
1335   Int_t iCount = 0;
1336   if(type==kTrackAOD){
1337     AliAODEvent *aod = 0;
1338     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1339     else aod = AODEvent();
1340     if(!aod){
1341       return iCount;
1342     }
1343     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1344       AliAODTrack *tr = aod->GetTrack(it);
1345       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1346       if(TMath::Abs(tr->Eta())>0.9)continue;
1347       if(tr->Pt()<fTrackPtCut)continue;
1348       list->Add(tr);
1349       iCount++;
1350     }
1351   }
1352   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1353     AliMCEvent* mcEvent = MCEvent();
1354     if(!mcEvent)return iCount;
1355     // we want to have alivpartilces so use get track
1356     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1357       if(!mcEvent->IsPhysicalPrimary(it))continue;
1358       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1359       if(type == kTrackKineAll){
1360         if(part->Pt()<fTrackPtCut)continue;
1361         list->Add(part);
1362         iCount++;
1363       }
1364       else if(type == kTrackKineCharged){
1365         if(part->Particle()->GetPDG()->Charge()==0)continue;
1366         if(part->Pt()<fTrackPtCut)continue;
1367         list->Add(part);
1368         iCount++;
1369       }
1370     }
1371   }
1372   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1373     AliAODEvent *aod = 0;
1374     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1375     else aod = AODEvent();
1376     if(!aod)return iCount;
1377     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1378     if(!tca)return iCount;
1379     for(int it = 0;it < tca->GetEntriesFast();++it){
1380       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1381       if(!part->IsPhysicalPrimary())continue;
1382       if(type == kTrackAODMCAll){
1383         if(part->Pt()<fTrackPtCut)continue;
1384         list->Add(part);
1385         iCount++;
1386       }
1387       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1388         if(part->Charge()==0)continue;
1389         if(part->Pt()<fTrackPtCut)continue;
1390         if(kTrackAODMCCharged){
1391           list->Add(part);
1392         }
1393         else {
1394           if(TMath::Abs(part->Eta())>0.9)continue;
1395           list->Add(part);
1396         }
1397         iCount++;
1398       }
1399     }
1400   }// AODMCparticle
1401   list->Sort();
1402   return iCount;
1403 }
1404
1405 /*
1406 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1407   for(int i = 0; i < particles.GetEntries(); i++){
1408     AliVParticle *vp = (AliVParticle*)particles.At(i);
1409     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1410     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1411     jInp.set_user_index(i);
1412     inputParticles.push_back(jInp);
1413   }
1414
1415   return 0;
1416
1417 }
1418 */