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
 908       if(tmpPt>fJetOutputMinPt&&jarray){// cut on the non-background subtracted...
 909         aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
 910         Double_t area1 = clustSeq.area(sortedJets[j]);
 911         aodOutJet->SetEffArea(area1,0);
 912       }
 913
 914
 915       Float_t tmpPtBack = 0;
 916       if(externalBackground){
 917         // carefull has to be filled in a task before
 918        // todo, ReArrange to the botom
 919         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
 920       }
 921       tmpPt = tmpPt - tmpPtBack;
 922       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
 923
 924       fh1PtJetsRecIn->Fill(tmpPt);
 925       // Fill Spectra with constituents
 926       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 927
 928       fh1NConstRec->Fill(constituents.size());
 929       fh2PtNch->Fill(nCh,tmpPt);
 930       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 931       fh2NConstPt->Fill(tmpPt,constituents.size());
 932       // loop over constiutents and fill spectrum
 933
 934       // Add the jet information and the track references to the output AOD
 935
 936
 937       if(fUseBackgroundCalc){
 938
 939         // create a random jet within the acceptance
 940         Double_t etaMax = 0.9 - fRparam;
 941         Int_t nCone = 0;
 942         Int_t nConeRan = 0;
 943         Double_t pTC = 1; // small number
 944         for(int ir = 0;ir < fNRandomCones;ir++){
 945           Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
 946           Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
 947           // massless jet
 948           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 949           Double_t pZC = pTC/TMath::Tan(thetaC);
 950           Double_t pXC = pTC * TMath::Cos(phiC);
 951           Double_t pYC = pTC * TMath::Sin(phiC);
 952           Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 953           AliAODJet tmpRecC (pXC,pYC,pZC, pC);
 954           bool skip = false;
 955           for(int jj = 0; jj < TMath::Min(nRec,2);jj++){
 956             AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
 957             if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
 958               skip = true;
 959               break;
 960             }
 961           }
 962           if(skip)continue;
 963           tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
 964           if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
 965           if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
 966         }// random cones
 967
 968
 969         // loop over the reconstructed particles and add up the pT in the random cones
 970         // maybe better to loop over randomized particles not in the real jets...
 971         // but this by definition brings dow average energy in the whole  event
 972         AliAODJet vTmpRanR(1,0,0,1);
 973         for(int i = 0; i < recParticles.GetEntries(); i++){
 974           AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 975           if(rConeArray){
 976             for(int ir = 0;ir < fNRandomCones;ir++){
 977               AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
 978               if(jC&&jC->DeltaR(vp)<fRparam){
 979                 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
 980               }
 981             }
 982           }// add up energy in cone
 983
 984           // the randomized input changes eta and phi, but keeps the p_T
 985           if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 986             Double_t pTR = vp->Pt();
 987             Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
 988             Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
 989
 990             Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
 991             Double_t pZR = pTR/TMath::Tan(thetaR);
 992
 993             Double_t pXR = pTR * TMath::Cos(phiR);
 994             Double_t pYR = pTR * TMath::Sin(phiR);
 995             Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
 996             vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
 997             if(rConeArrayRan){
 998               for(int ir = 0;ir < fNRandomCones;ir++){
 999                 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
1000                 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1001                   jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1002                 }
1003               }
1004             }
1005           }
1006         }// loop over recparticles
1007
1008         Float_t jetArea = fRparam*fRparam*TMath::Pi();
1009         for(int ir = 0;ir < fNRandomCones;ir++){
1010           // rescale the momntum vectors for the random cones
1011           if(!rConeArray)continue;
1012           AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
1013           if(rC){
1014             Double_t etaC = rC->Eta();
1015             Double_t phiC = rC->Phi();
1016             // massless jet, unit vector
1017             pTC = rC->ChargedBgEnergy();
1018             if(pTC<=0)pTC = 0.1; // for almost empty events
1019             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1020             Double_t pZC = pTC/TMath::Tan(thetaC);
1021             Double_t pXC = pTC * TMath::Cos(phiC);
1022             Double_t pYC = pTC * TMath::Sin(phiC);
1023             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1024             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1025             rC->SetBgEnergy(0,0);
1026             rC->SetEffArea(jetArea,0);
1027           }
1028           rC = (AliAODJet*)rConeArrayRan->At(ir);
1029           // same wit random
1030           if(rC){
1031             Double_t etaC = rC->Eta();
1032             Double_t phiC = rC->Phi();
1033             // massless jet, unit vector
1034             pTC = rC->ChargedBgEnergy();
1035             if(pTC<=0)pTC = 0.1;// for almost empty events
1036             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1037             Double_t pZC = pTC/TMath::Tan(thetaC);
1038             Double_t pXC = pTC * TMath::Cos(phiC);
1039             Double_t pYC = pTC * TMath::Sin(phiC);
1040             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1041             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1042             rC->SetBgEnergy(0,0);
1043             rC->SetEffArea(jetArea,0);
1044           }
1045         }
1046       }// if(fUseBackgroundCalc
1047
1048      for(unsigned int ic = 0; ic < constituents.size();ic++){
1049        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1050        fh1PtJetConstRec->Fill(part->Pt());
1051        if(aodOutJet){
1052          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1053        }
1054        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1055      }
1056
1057      // correlation
1058      Float_t tmpPhi =  tmpRec.Phi();
1059      Float_t tmpEta =  tmpRec.Eta();
1060      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1061
1062
1063
1064      if(j==0){
1065        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1066        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1067        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1068        fh1NConstLeadingRec->Fill(constituents.size());
1069        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1070        continue;
1071      }
1072      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1073      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1074      Float_t dPhi = phi - tmpPhi;
1075      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1076      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1077      Float_t dEta = eta - tmpRec.Eta();
1078      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1079      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1080      if(cenClass>=0){
1081        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1082        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1083      }
1084      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1085     }
1086    delete recIter;
1087
1088    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1089
1090
1091    if(evBkg){
1092      vector<fastjet::PseudoJet> jets2=sortedJets;
1093      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1094      Double_t bkg1=0;
1095      Double_t sigma1=0.;
1096      Double_t meanarea1=0.;
1097      Double_t bkg2=0;
1098      Double_t sigma2=0.;
1099      Double_t meanarea2=0.;
1100
1101      clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
1102      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
1103
1104      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1105      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1106
1107      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1108      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
1109      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1110      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1111
1112    }
1113   }
1114
1115
1116
1117
1118
1119   // fill track information
1120   Int_t nTrackOver = recParticles.GetSize();
1121   // do the same for tracks and jets
1122
1123   if(nTrackOver>0){
1124    TIterator *recIter = recParticles.MakeIterator();
1125    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1126    Float_t pt = tmpRec->Pt();
1127
1128    //    Printf("Leading track p_t %3.3E",pt);
1129    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1130      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1131      while(pt<ptCut&&tmpRec){
1132        nTrackOver--;
1133        tmpRec = (AliVParticle*)(recIter->Next());
1134        if(tmpRec){
1135          pt = tmpRec->Pt();
1136        }
1137      }
1138      if(nTrackOver<=0)break;
1139      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1140    }
1141
1142    recIter->Reset();
1143    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1144    Float_t phi = leading->Phi();
1145    if(phi<0)phi+=TMath::Pi()*2.;
1146    Float_t eta = leading->Eta();
1147    pt  = leading->Pt();
1148    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1149      Float_t tmpPt = tmpRec->Pt();
1150      Float_t tmpEta = tmpRec->Eta();
1151      fh1PtTracksRecIn->Fill(tmpPt);
1152      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1153      if(tmpRec==leading){
1154        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1155        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1156        continue;
1157      }
1158       // correlation
1159      Float_t tmpPhi =  tmpRec->Phi();
1160
1161      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1162      Float_t dPhi = phi - tmpPhi;
1163      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1164      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1165      Float_t dEta = eta - tmpRec->Eta();
1166      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1167      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1168      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1169    }
1170    delete recIter;
1171  }
1172
1173  // find the random jets
1174  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1175  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1176
1177  inclusiveJetsRan = clustSeqRan.inclusive_jets();
1178  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1179
1180  // fill the jet information from random track
1181
1182   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1183
1184  // loop over all jets an fill information, first on is the leading jet
1185
1186  Int_t nRecOverRan = inclusiveJetsRan.size();
1187  Int_t nRecRan     = inclusiveJetsRan.size();
1188  if(inclusiveJetsRan.size()>0){
1189    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1190    Float_t pt = leadingJet.Pt();
1191
1192    Int_t iCount = 0;
1193    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1194      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1195       while(pt<ptCut&&iCount<nRecRan){
1196         nRecOverRan--;
1197         iCount++;
1198         if(iCount<nRecRan){
1199           pt = sortedJetsRan[iCount].perp();
1200         }
1201       }
1202       if(nRecOverRan<=0)break;
1203       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1204     }
1205     Float_t phi = leadingJet.Phi();
1206     if(phi<0)phi+=TMath::Pi()*2.;
1207     pt = leadingJet.Pt();
1208
1209     // correlation of leading jet with random tracks
1210
1211     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1212       {
1213         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1214         // correlation
1215         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1216         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1217         Float_t dPhi = phi - tmpPhi;
1218         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1219         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1220         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1221         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1222       }
1223
1224     Int_t nAodOutJetsRan = 0;
1225      AliAODJet *aodOutJetRan = 0;
1226     for(int j = 0; j < nRecRan;j++){
1227       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1228       Float_t tmpPt = tmpRec.Pt();
1229       fh1PtJetsRecInRan->Fill(tmpPt);
1230       // Fill Spectra with constituents
1231       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1232       fh1NConstRecRan->Fill(constituents.size());
1233       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1234       fh2PtNchRan->Fill(nCh,tmpPt);
1235       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1236
1237
1238      if(tmpPt>fJetOutputMinPt&&jarrayran){
1239        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1240        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1241
1242        aodOutJetRan->SetEffArea(arearan,0);    }
1243
1244
1245
1246
1247      for(unsigned int ic = 0; ic < constituents.size();ic++){
1248        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1249        // fh1PtJetConstRec->Fill(part->Pt());
1250        if(aodOutJetRan){
1251          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1252        }
1253      }
1254
1255
1256
1257
1258       // correlation
1259       Float_t tmpPhi =  tmpRec.Phi();
1260       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1261
1262       if(j==0){
1263         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1264         fh1NConstLeadingRecRan->Fill(constituents.size());
1265         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1266         continue;
1267       }
1268     }
1269
1270
1271     if(evBkg){
1272      Double_t bkg3=0.;
1273      Double_t sigma3=0.;
1274      Double_t meanarea3=0.;
1275      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1276      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1277      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1278      /*
1279      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1280      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1281      */
1282     }
1283
1284
1285
1286  }
1287
1288
1289  // do the event selection if activated
1290  if(fJetTriggerPtCut>0){
1291    bool select = false;
1292    Float_t minPt = fJetTriggerPtCut;
1293    /*
1294    // hard coded for now ...
1295    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1296    if(cent<10)minPt = 50;
1297    else if(cent<30)minPt = 42;
1298    else if(cent<50)minPt = 28;
1299    else if(cent<80)minPt = 18;
1300    */
1301    float rho = 0;
1302    if(externalBackground)rho = externalBackground->GetBackground(2);
1303    if(jarray){
1304      for(int i = 0;i < jarray->GetEntriesFast();i++){
1305        AliAODJet *jet = (AliAODJet*)jarray->At(i);
1306        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1307        if(ptSub>=minPt){
1308          select = true;
1309          break;
1310        }
1311      }
1312    }
1313
1314    if(select){
1315      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1316      fh1CentralitySelect->Fill(cent);
1317      fh1ZSelect->Fill(zVtx);
1318      aodH->SetFillAOD(kTRUE);
1319    }
1320  }
1321
1322  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1323  PostData(1, fHistList);
1324 }
1325
1326 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1327 {
1328 // Terminate analysis
1329 //
1330     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1331 }
1332
1333
1334 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1335
1336   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1337
1338   Int_t iCount = 0;
1339   if(type==kTrackAOD){
1340     AliAODEvent *aod = 0;
1341     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1342     else aod = AODEvent();
1343     if(!aod){
1344       return iCount;
1345     }
1346     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1347       AliAODTrack *tr = aod->GetTrack(it);
1348       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1349       if(TMath::Abs(tr->Eta())>0.9)continue;
1350       if(tr->Pt()<fTrackPtCut)continue;
1351       list->Add(tr);
1352       iCount++;
1353     }
1354   }
1355   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1356     AliMCEvent* mcEvent = MCEvent();
1357     if(!mcEvent)return iCount;
1358     // we want to have alivpartilces so use get track
1359     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1360       if(!mcEvent->IsPhysicalPrimary(it))continue;
1361       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1362       if(type == kTrackKineAll){
1363         if(part->Pt()<fTrackPtCut)continue;
1364         list->Add(part);
1365         iCount++;
1366       }
1367       else if(type == kTrackKineCharged){
1368         if(part->Particle()->GetPDG()->Charge()==0)continue;
1369         if(part->Pt()<fTrackPtCut)continue;
1370         list->Add(part);
1371         iCount++;
1372       }
1373     }
1374   }
1375   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1376     AliAODEvent *aod = 0;
1377     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1378     else aod = AODEvent();
1379     if(!aod)return iCount;
1380     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1381     if(!tca)return iCount;
1382     for(int it = 0;it < tca->GetEntriesFast();++it){
1383       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1384       if(!part->IsPhysicalPrimary())continue;
1385       if(type == kTrackAODMCAll){
1386         if(part->Pt()<fTrackPtCut)continue;
1387         list->Add(part);
1388         iCount++;
1389       }
1390       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1391         if(part->Charge()==0)continue;
1392         if(part->Pt()<fTrackPtCut)continue;
1393         if(kTrackAODMCCharged){
1394           list->Add(part);
1395         }
1396         else {
1397           if(TMath::Abs(part->Eta())>0.9)continue;
1398           list->Add(part);
1399         }
1400         iCount++;
1401       }
1402     }
1403   }// AODMCparticle
1404   list->Sort();
1405   return iCount;
1406 }
1407
1408 /*
1409 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1410   for(int i = 0; i < particles.GetEntries(); i++){
1411     AliVParticle *vp = (AliVParticle*)particles.At(i);
1412     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1413     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1414     jInp.set_user_index(i);
1415     inputParticles.push_back(jInp);
1416   }
1417
1418   return 0;
1419
1420 }
1421 */