JETAN/AliAnalysisTaskJetClusterKine.cxx

   1 // **************************************
   2 // Task used for jet finding in the Kine Train (generation and analysis on the fly, no detector effects)
   3 // Output is stored in an exchance container
   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 <TF1.h>
  36 #include <TList.h>
  37 #include <TLorentzVector.h>
  38 #include <TClonesArray.h>
  39 #include  "TDatabasePDG.h"
  40 #include <TGrid.h>
  41
  42 #include "AliAnalysisTaskJetClusterKine.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliJetFinder.h"
  45 #include "AliJetHeader.h"
  46 #include "AliJetReader.h"
  47 #include "AliESDEvent.h"
  48 #include "AliAODEvent.h"
  49 #include "AliAODHandler.h"
  50 #include "AliAODExtension.h"
  51 #include "AliAODTrack.h"
  52 #include "AliAODJet.h"
  53 #include "AliVParticle.h" //FK//
  54 #include "AliAODMCParticle.h"
  55 #include "AliMCEventHandler.h"
  56 #include "AliMCEvent.h"
  57 #include "AliStack.h"
  58 #include "AliGenEventHeader.h" //FK//
  59 #include "AliGenPythiaEventHeader.h"
  60 #include "AliJetKineReaderHeader.h"
  61 #include "AliGenCocktailEventHeader.h"
  62 #include "AliInputEventHandler.h"
  63 #include "AliAODJetEventBackground.h"
  64
  65 #include "fastjet/PseudoJet.hh"
  66 #include "fastjet/ClusterSequenceArea.hh"
  67 #include "fastjet/AreaDefinition.hh"
  68 #include "fastjet/JetDefinition.hh"
  69 // get info on how fastjet was configured
  70 #include "fastjet/config.h"
  71
  72 using std::vector;
  73
  74 ClassImp(AliAnalysisTaskJetClusterKine)
  75
  76 AliAnalysisTaskJetClusterKine::~AliAnalysisTaskJetClusterKine(){
  77   //
  78   // Destructor
  79   //
  80
  81   delete fRef;
  82
  83   if(fTCAJetsOut)fTCAJetsOut->Delete();
  84   delete fTCAJetsOut;
  85
  86 }
  87
  88 //_____________________________________________________________________
  89
  90 AliAnalysisTaskJetClusterKine::AliAnalysisTaskJetClusterKine():
  91   AliAnalysisTaskSE(),
  92   fMcEvent(0x0),
  93   fMcHandler(0x0),
  94   fRef(new TRefArray),
  95   fTrackTypeGen(kTrackKineCharged), //Kine charged?
  96   fAvgTrials(1),
  97   fTrackEtaWindow(0.9),
  98   fTrackPtCut(0.),
  99   fJetOutputMinPt(0.150),
 100   fMaxTrackPtInJet(100.),
 101   fVtxZCut(10.0),
 102   fNonStdBranch(""),
 103   fOutContainer(kNoOutput), //FF//
 104   fNonStdFile(""),
 105   fRparam(1.0),
 106   fAlgorithm(fastjet::kt_algorithm),
 107   fStrategy(fastjet::Best),
 108   fRecombScheme(fastjet::BIpt_scheme),
 109   fAreaType(fastjet::active_area),
 110   fGhostArea(0.01),
 111   fActiveAreaRepeats(1),
 112   fGhostEtamax(1.5),
 113   fTCAJetsOut(0x0),
 114   fh1Xsec(0x0),
 115   fh1Trials(0x0),
 116   fh1PtHard(0x0),
 117   fh1PtHardNoW(0x0),
 118   fh1PtHardTrials(0x0),
 119   fh1NJetsGen(0x0),
 120   fh1NConstGen(0x0),
 121   fh1NConstLeadingGen(0x0),
 122   fh1PtJetsGenIn(0x0),
 123   fh1PtJetsLeadingGenIn(0x0),
 124   fh1PtJetConstGen(0x0),
 125   fh1PtJetConstLeadingGen(0x0),
 126   fh1PtTracksGenIn(0x0),
 127   fh1Nch(0x0),
 128   fh1Z(0x0),
 129   fh2NConstPt(0x0),
 130   fh2NConstLeadingPt(0x0),
 131   fh2JetPhiEta(0x0),
 132   fh2LeadingJetPhiEta(0x0),
 133   fh2JetEtaPt(0x0),
 134   fh2LeadingJetEtaPt(0x0),
 135   fh2TrackEtaPt(0x0),
 136   fh2JetsLeadingPhiEta(0x0),
 137   fh2JetsLeadingPhiPt(0x0),
 138   fh2JetsLeadingPhiPtW(0x0),
 139   fHistList(0x0)
 140 {
 141   //
 142   // Constructor
 143   //
 144 }
 145
 146 //_____________________________________________________________________
 147
 148 AliAnalysisTaskJetClusterKine::AliAnalysisTaskJetClusterKine(const char* name):
 149   AliAnalysisTaskSE(name),
 150   fMcEvent(0x0),
 151   fMcHandler(0x0),
 152   fRef(new TRefArray),
 153   fTrackTypeGen(kTrackKineCharged),  //kine charged?
 154   fAvgTrials(1),
 155   fTrackEtaWindow(0.9),
 156   fTrackPtCut(0.),
 157   fJetOutputMinPt(0.150),
 158   fMaxTrackPtInJet(100.),
 159   fVtxZCut(10.0),
 160   fNonStdBranch(""),
 161   fOutContainer(kNoOutput),//FF//
 162   fNonStdFile(""),
 163   fRparam(1.0),
 164   fAlgorithm(fastjet::kt_algorithm),
 165   fStrategy(fastjet::Best),
 166   fRecombScheme(fastjet::BIpt_scheme),
 167   fAreaType(fastjet::active_area),
 168   fGhostArea(0.01),
 169   fActiveAreaRepeats(1),
 170   fGhostEtamax(1.5),
 171   fTCAJetsOut(0x0),
 172   fh1Xsec(0x0),
 173   fh1Trials(0x0),
 174   fh1PtHard(0x0),
 175   fh1PtHardNoW(0x0),
 176   fh1PtHardTrials(0x0),
 177   fh1NJetsGen(0x0),
 178   fh1NConstGen(0x0),
 179   fh1NConstLeadingGen(0x0),
 180   fh1PtJetsGenIn(0x0),
 181   fh1PtJetsLeadingGenIn(0x0),
 182   fh1PtJetConstGen(0x0),
 183   fh1PtJetConstLeadingGen(0x0),
 184   fh1PtTracksGenIn(0x0),
 185   fh1Nch(0x0),
 186   fh1Z(0x0),
 187   fh2NConstPt(0x0),
 188   fh2NConstLeadingPt(0x0),
 189   fh2JetPhiEta(0x0),
 190   fh2LeadingJetPhiEta(0x0),
 191   fh2JetEtaPt(0x0),
 192   fh2LeadingJetEtaPt(0x0),
 193   fh2TrackEtaPt(0x0),
 194   fh2JetsLeadingPhiEta(0x0),
 195   fh2JetsLeadingPhiPt(0x0),
 196   fh2JetsLeadingPhiPtW(0x0),
 197   fHistList(0x0)
 198 {
 199   //
 200   // named ctor
 201   //
 202   DefineOutput(1, TList::Class());
 203   DefineOutput(2, TClonesArray::Class());
 204 }
 205
 206
 207 //_____________________________________________________________________
 208
 209 Bool_t AliAnalysisTaskJetClusterKine::Notify(){
 210   //
 211   // Implemented Notify() to read the cross sections
 212   // and number of trials from pyxsec.root
 213   //
 214   return kTRUE;
 215 }
 216
 217 //_____________________________________________________________________
 218
 219 void AliAnalysisTaskJetClusterKine::UserCreateOutputObjects(){
 220
 221    //
 222    // Create the output container
 223    //
 224
 225    // Connect the AOD
 226
 227
 228    if(fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 229
 230
 231    if(fNonStdBranch.Length()!=0){
 232       // only create the output branch if we have a name
 233       // Create a new branch for jets...
 234       //  -> cleared in the UserExec....
 235       // here we can also have the case that the brnaches are written to a separate file
 236       fTCAJetsOut = new TClonesArray("AliAODJet", 0);
 237       fTCAJetsOut->SetName(fNonStdBranch.Data());
 238       if(fOutContainer==kAODBranch){   //FF//
 239          AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
 240       }
 241
 242
 243       //if(fNonStdFile.Length()!=0){
 244         //
 245         // case that we have an AOD extension we need to fetch the jets from the extended output
 246         // we identify the extension aod event by looking for the branchname
 247         //AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 248         // case that we have an AOD extension we need can fetch the background maybe from the extended output
 249         //fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 250       //}
 251    }
 252
 253
 254    if(!fHistList) fHistList = new TList();
 255    fHistList->SetOwner(kTRUE);
 256    PostData(1, fHistList); // post data in any case once
 257
 258
 259    if(fOutContainer==kExchCont){
 260       fTCAJetsOut->SetOwner();
 261       PostData(2, fTCAJetsOut); //FF// post data in any case once
 262    }
 263
 264    Bool_t oldStatus = TH1::AddDirectoryStatus();
 265    TH1::AddDirectory(kFALSE);
 266
 267
 268    //
 269    //  Histogram
 270
 271    const Int_t nBinPt = 100;
 272    Double_t binLimitsPt[nBinPt+1];
 273    for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 274       if(iPt == 0){
 275          binLimitsPt[iPt] = 0.0;
 276       }else {// 1.0
 277          binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
 278       }
 279    }
 280
 281    const Int_t nBinPhi = 90;
 282    Double_t binLimitsPhi[nBinPhi+1];
 283    for(Int_t iPhi = 0; iPhi<=nBinPhi; iPhi++){
 284        if(iPhi==0){
 285           binLimitsPhi[iPhi] = -1.*TMath::Pi();
 286        }else{
 287           binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 288        }
 289     }
 290
 291     const Int_t nBinEta = 40;
 292     Double_t binLimitsEta[nBinEta+1];
 293     for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 294        if(iEta==0){
 295           binLimitsEta[iEta] = -2.0;
 296        }else{
 297           binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 298      }
 299    }
 300
 301    const int nChMax = 5000;
 302
 303    fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 304    fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 305
 306    fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 307    fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 308
 309
 310    fh1NJetsGen = new TH1F("fh1NJetsGen","N reconstructed jets",120,-0.5,119.5);
 311
 312    fh1NConstGen = new TH1F("fh1NConstGen","# jet constituents",120,-0.5,119.5);
 313    fh1NConstLeadingGen = new TH1F("fh1NConstLeadingGen","jet constituents",120,-0.5,119.5);
 314
 315
 316    fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 317    fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 318    fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 319
 320    fh1PtJetsGenIn  = new TH1F("fh1PtJetsGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 321    fh1PtJetsLeadingGenIn = new TH1F("fh1PtJetsLeadingGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 322    fh1PtJetConstGen = new TH1F("fh1PtJetsConstGen","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 323    fh1PtJetConstLeadingGen = new TH1F("fh1PtJetsConstLeadingGen","Gen jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 324    fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("Gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
 325    fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 326
 327
 328    fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 329
 330
 331
 332    fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 333    fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 334
 335
 336    fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 337                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 338    fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 339                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 340
 341    fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 342                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 343    fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 344                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 345
 346    fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 347                             nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 348
 349
 350
 351    fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 352                                   nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 353    fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 354                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 355
 356    fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 357                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 358
 359
 360
 361
 362    const Int_t saveLevel = 3; // large save level more histos
 363    if(saveLevel>0){
 364       fHistList->Add(fh1Xsec);
 365       fHistList->Add(fh1Trials);
 366
 367       fHistList->Add(fh1NJetsGen);
 368       fHistList->Add(fh1NConstGen);
 369       fHistList->Add(fh1NConstLeadingGen);
 370       fHistList->Add(fh1PtJetsGenIn);
 371       fHistList->Add(fh1PtJetsLeadingGenIn);
 372       fHistList->Add(fh1PtTracksGenIn);
 373       fHistList->Add(fh1PtJetConstGen);
 374       fHistList->Add(fh1PtJetConstLeadingGen);
 375       fHistList->Add(fh1Nch);
 376       fHistList->Add(fh1Z);
 377
 378       fHistList->Add(fh2NConstPt);
 379       fHistList->Add(fh2NConstLeadingPt);
 380       fHistList->Add(fh2JetPhiEta);
 381       fHistList->Add(fh2LeadingJetPhiEta);
 382       fHistList->Add(fh2JetEtaPt);
 383       fHistList->Add(fh2LeadingJetEtaPt);
 384       fHistList->Add(fh2TrackEtaPt);
 385       fHistList->Add(fh2JetsLeadingPhiEta);
 386       fHistList->Add(fh2JetsLeadingPhiPt);
 387       fHistList->Add(fh2JetsLeadingPhiPtW);
 388    }
 389
 390    // =========== Switch on Sumw2 for all histos ===========
 391    for(Int_t i=0; i<fHistList->GetEntries(); ++i){
 392       TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 393       if(h1){
 394          h1->Sumw2();
 395          continue;
 396       }
 397       THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 398       if(hn) hn->Sumw2();
 399    }
 400    TH1::AddDirectory(oldStatus);
 401 }
 402 //_________________________________________________________________________
 403
 404 void AliAnalysisTaskJetClusterKine::Init(){
 405    // MC handler
 406    if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Init() \n");
 407    fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); //FK//
 408 }
 409
 410 //_________________________________________________________________________
 411 void AliAnalysisTaskJetClusterKine::LocalInit(){
 412    // MC handler
 413    //
 414    // Initialization
 415    //
 416
 417    if(fDebug > 1) printf("AnalysisTaskJetClusterKine::LocalInit() \n");
 418    Init();
 419 }
 420
 421 //_________________________________________________________________________
 422 void AliAnalysisTaskJetClusterKine::UserExec(Option_t* /*option*/){
 423    // handle and reset the output jet branch
 424
 425    if(fDebug > 1) printf("AliAnalysisTaskJetClusterKine::UserExec \n");
 426    if(fTCAJetsOut) fTCAJetsOut->Delete(); //clean TClonesArray
 427
 428    //
 429    // Execute analysis for current event
 430    //
 431    Init();
 432    if(fMcHandler){
 433       fMcEvent = fMcHandler->MCEvent();
 434    }else{
 435       if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Handler() fMcHandler=NULL\n");
 436       PostData(1, fHistList);
 437
 438       if(fOutContainer==kExchCont){
 439          PostData(2, fTCAJetsOut); //FF//
 440       }
 441
 442       return;
 443    }
 444    if(!fMcEvent){
 445       if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Exec()   fMcEvent=NULL \n");
 446       PostData(1, fHistList);
 447
 448       if(fOutContainer==kExchCont){
 449          PostData(2, fTCAJetsOut); //FF//
 450       }
 451
 452       return;
 453    }
 454
 455    const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
 456    Float_t zVtx = vtxMC->GetZ();
 457    if(TMath::Abs(zVtx) > fVtxZCut){ //vertex cut
 458       PostData(1, fHistList);
 459
 460       if(fOutContainer==kExchCont){
 461          PostData(2, fTCAJetsOut); //FF//
 462       }
 463
 464       return;
 465    }
 466
 467    fh1Z->Fill(zVtx);
 468
 469    fh1Trials->Fill("#sum{ntrials}",1);
 470    if(fDebug > 10) Printf("%s:%d",(char*)__FILE__,__LINE__);
 471
 472    // we simply fetch the mc particles as a list of AliVParticles
 473
 474    TList genParticles; //list of particles  with above min pT and  good eta range
 475
 476    Int_t nParticles = GetListOfTracks(&genParticles, fTrackTypeGen); //fill the list
 477    fh1Nch->Fill(nParticles);
 478
 479    if(fDebug>2) Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__, nParticles, genParticles.GetEntries());
 480
 481    // find the jets....
 482
 483    vector<fastjet::PseudoJet> inputParticlesKine;
 484
 485    for(int i = 0; i < genParticles.GetEntries(); i++){  //loop over generated particles
 486       AliVParticle *vp = (AliVParticle*) genParticles.At(i);
 487       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 488       jInp.set_user_index(i);
 489       inputParticlesKine.push_back(jInp);
 490
 491       if(fTCAJetsOut){
 492          if(i == 0){
 493             fRef->Delete(); // make sure to delete before placement new...
 494             new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
 495          }
 496          fRef->Add(vp);  //TRefArray does not work with toy model ...
 497       }
 498    }//generated particles
 499
 500    if(inputParticlesKine.size()==0){ //FK//
 501       if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 502       PostData(1, fHistList);
 503
 504       if(fOutContainer==kExchCont){
 505          PostData(2, fTCAJetsOut); //FF//
 506       }
 507
 508       return;
 509    } //FK//
 510
 511    // run fast jet
 512    // employ setters for these...
 513    // now create the object that holds info about ghosts
 514
 515    fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 516    fastjet::AreaType areaType =   fastjet::active_area;
 517    fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 518    fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 519    fastjet::ClusterSequenceArea clustSeq(inputParticlesKine, jetDef, areaDef); //FK//
 520
 521    //range where to compute background
 522    Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 523    phiMin = 0;
 524    phiMax = 2*TMath::Pi();
 525    rapMax = fGhostEtamax - fRparam;
 526    rapMin = - fGhostEtamax + fRparam;
 527    fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 528
 529
 530    const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
 531    const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
 532
 533
 534    fh1NJetsGen->Fill(sortedJets.size()); //the number of found jets
 535
 536    // loop over all jets an fill information, first on is the leading jet
 537
 538    Int_t nGen     = inclusiveJets.size();
 539    if(inclusiveJets.size()>0){
 540
 541       //leading Jet
 542       AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 543       Double_t area = clustSeq.area(sortedJets[0]);
 544       leadingJet.SetEffArea(area,0);
 545       Float_t pt = leadingJet.Pt();
 546       Float_t phi = leadingJet.Phi();
 547       if(phi<0) phi += TMath::Pi()*2.;
 548       Float_t eta = leadingJet.Eta();
 549
 550       //inclusive jet
 551       Int_t nAodOutJets = 0;
 552       AliAODJet *aodOutJet = NULL;
 553
 554
 555       // correlation of leading jet with tracks
 556   //FK//    TIterator *particleIter = genParticles.MakeIterator();//FK//
 557   //FK//    particleIter->Reset();
 558   //FK//    AliVParticle *tmpParticle = NULL;
 559    //FK//   while((tmpParticle = (AliVParticle*)(particleIter->Next()))){
 560    //FK//      Float_t tmpPt = tmpParticle->Pt();
 561     //FK//     // correlation
 562     //FK//     Float_t tmpPhi =  tmpParticle->Phi();
 563    //FK//      if(tmpPhi<0) tmpPhi+= TMath::Pi()*2.;
 564    //FK//      Float_t dPhi = phi - tmpPhi;
 565      //FK//    if(dPhi > TMath::Pi())    dPhi = dPhi - 2.*TMath::Pi(); //-pi,pi
 566     //FK//     if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();  //-pi,pi
 567      //FK//    fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 568      //FK//    fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 569     //FK//  }
 570
 571       TLorentzVector vecareab;
 572       for(int j = 0; j < nGen;j++){ //loop over inclusive jets
 573          AliAODJet tmpGenJet (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 574          aodOutJet = NULL;
 575          Float_t tmpPt = tmpGenJet.Pt(); //incl jet pt
 576
 577          if((tmpPt > fJetOutputMinPt) && fTCAJetsOut){// cut on the non-background subtracted...
 578             aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpGenJet);
 579             aodOutJet->GetRefTracks()->Clear();
 580             Double_t area1 = clustSeq.area(sortedJets[j]);
 581             aodOutJet->SetEffArea(area1,0);
 582             fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
 583             vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
 584             aodOutJet->SetVectorAreaCharged(&vecareab);
 585          }
 586
 587          fh1PtJetsGenIn->Fill(tmpPt); //incl jet Pt
 588          // Fill Spectra with constituentsemacs
 589          const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
 590
 591          fh1NConstGen->Fill(constituents.size());  //number of constituents
 592          fh2NConstPt->Fill(tmpPt,constituents.size()); //number of constituents vs jet pt
 593
 594          // loop over constiutents and fill spectrum
 595          AliVParticle *partLead = 0;
 596          Float_t ptLead = -1;
 597
 598          for(unsigned int ic = 0; ic < constituents.size(); ic++){
 599             AliVParticle *part = (AliVParticle*) genParticles.At(constituents[ic].user_index());
 600             if(!part) continue;
 601             fh1PtJetConstGen->Fill(part->Pt()); //pt of constituents
 602
 603             if(aodOutJet){
 604                aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));//FK//
 605
 606                if(part->Pt()>fMaxTrackPtInJet){
 607                   aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
 608                }
 609             }
 610
 611             if(part->Pt()>ptLead){
 612                ptLead = part->Pt();
 613                partLead = part;
 614             }
 615
 616             if(j==0) fh1PtJetConstLeadingGen->Fill(part->Pt()); //pt of leading jet constituents
 617          }
 618
 619          if(partLead){
 620             if(aodOutJet){
 621                //set pT of leading constituent of jet
 622                aodOutJet->SetPtLeading(partLead->Pt());
 623             }
 624          }
 625
 626          // correlation
 627          Float_t tmpPhi =  tmpGenJet.Phi(); //incl jet phi
 628          Float_t tmpEta =  tmpGenJet.Eta(); //incl jet eta
 629
 630          if(tmpPhi<0) tmpPhi += TMath::Pi()*2.;
 631          fh2JetPhiEta->Fill(tmpGenJet.Phi(),tmpEta);
 632          fh2JetEtaPt->Fill(tmpEta,tmpPt);
 633
 634          if(j==0){   //leading jet
 635             fh1PtJetsLeadingGenIn->Fill(tmpPt); //leading jet pt
 636             fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
 637             fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
 638             fh1NConstLeadingGen->Fill(constituents.size());  //number of constituents in leading jet
 639             fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
 640             continue;
 641          }
 642
 643          Float_t dPhi = phi - tmpPhi;
 644          if(dPhi > TMath::Pi())     dPhi = dPhi - 2.*TMath::Pi(); //-pi,pi
 645          if(dPhi<(-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi(); //-pi.pi
 646          Float_t dEta = eta - tmpGenJet.Eta();
 647          fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
 648          fh2JetsLeadingPhiPt->Fill(dPhi,pt);
 649          fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
 650       }// loop over reconstructed jets
 651     //delete particleIter;
 652    }//number of jets>0
 653
 654
 655    // fill track information
 656    Int_t nTrackOver = genParticles.GetSize();
 657    // do the same for tracks and jets
 658
 659    if(nTrackOver>0){
 660       TIterator *particleIter = genParticles.MakeIterator();
 661       AliVParticle *tmpGen = 0;
 662       particleIter->Reset();
 663
 664       while((tmpGen = (AliVParticle*)(particleIter->Next()))){
 665          Float_t tmpPt = tmpGen->Pt();
 666          Float_t tmpEta = tmpGen->Eta();
 667          fh1PtTracksGenIn->Fill(tmpPt);
 668          fh2TrackEtaPt->Fill(tmpEta,tmpPt);
 669       }
 670       delete particleIter;
 671    }
 672
 673
 674
 675
 676
 677    if(fDebug > 2){
 678       if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
 679    }
 680    PostData(1, fHistList);
 681
 682   if(fOutContainer==kExchCont){
 683       PostData(2, fTCAJetsOut); //FF//
 684    }
 685
 686
 687 }
 688 //____________________________________________________________________________
 689
 690 void AliAnalysisTaskJetClusterKine::Terminate(Option_t* /*option*/){
 691   //
 692   // Terminate analysis
 693   //
 694   if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
 695
 696
 697 }
 698
 699 //_____________________________________________________________________________________
 700
 701 Int_t  AliAnalysisTaskJetClusterKine::GetListOfTracks(TList *list, Int_t type){
 702
 703   //
 704   // get list of tracks/particles for different types
 705   //
 706
 707    if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
 708
 709    Int_t iCount = 0; //number of particles
 710
 711    if(type ==  kTrackKineAll || type == kTrackKineCharged){
 712       if(! fMcEvent) return iCount;
 713
 714       //we want to have alivpartilces so use get track
 715       for(int it = 0; it < fMcEvent->GetNumberOfTracks(); ++it){
 716          if(!fMcEvent->IsPhysicalPrimary(it)) continue;
 717
 718          AliMCParticle* part = (AliMCParticle*)fMcEvent->GetTrack(it);
 719          if(TMath::Abs(part->Eta()) > fTrackEtaWindow) continue;   //apply eta cut
 720          if(part->Pt() < fTrackPtCut) continue;  //apply pT cut
 721
 722          if(type == kTrackKineAll){  // full jets
 723             list->Add(part);
 724             iCount++;
 725
 726          }else if(type == kTrackKineCharged){  //charged jets
 727             if(part->Particle()->GetPDG()->Charge()==0) continue;
 728             list->Add(part);
 729             iCount++;
 730          }
 731       }
 732    }
 733
 734    list->Sort();  //?//
 735
 736    return iCount; //number of particles in the list
 737 }
 738
 739