JETAN/AliSISConeJetFinder.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 //---------------------------------------------------------------------
  18 // FastJet v2.3.4 finder algorithm interface
  19 //
  20 // Author: swensy.jangal@ires.in2p3.fr
  21 //
  22 // Last modification: Neutral cell energy included in the jet reconstruction
  23 //
  24 // Author: Magali.estienne@subatech.in2p3.fr
  25 //---------------------------------------------------------------------
  26
  27
  28 #include <Riostream.h>
  29 #include <TArrayF.h>
  30 #include <TFile.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TLorentzVector.h>
  34 #include <TRandom.h>
  35 #include <TClonesArray.h>
  36
  37 #include "AliHeader.h"
  38 #include "AliJet.h"
  39 #include "AliJetKineReader.h"
  40 #include "AliJetReader.h"
  41 #include "AliJetReaderHeader.h"
  42 #include "AliJetUnitArray.h"
  43 #include "AliSISConeJetFinder.h"
  44 #include "AliSISConeJetHeader.h"
  45
  46 #include "fastjet/AreaDefinition.hh"
  47 #include "fastjet/ClusterSequenceArea.hh"
  48 #include "fastjet/JetDefinition.hh"
  49 #include "fastjet/PseudoJet.hh"
  50
  51 // get info on how fastjet was configured
  52 #include "fastjet/config.h"
  53
  54 #ifdef ENABLE_PLUGIN_SISCONE
  55 #include "fastjet/SISConePlugin.hh"
  56 #endif
  57
  58 #include<sstream>  // needed for internal io
  59 #include<vector>
  60 #include <cmath>
  61
  62 using namespace std;
  63
  64
  65 ClassImp(AliSISConeJetFinder)
  66
  67
  68 //____________________________________________________________________________
  69
  70 AliSISConeJetFinder::AliSISConeJetFinder():
  71 AliJetFinder()
  72 {
  73   // Constructor
  74 }
  75
  76 //____________________________________________________________________________
  77
  78 AliSISConeJetFinder::~AliSISConeJetFinder()
  79 {
  80   // destructor
  81 }
  82
  83 //______________________________________________________________________________
  84 void AliSISConeJetFinder::FindJets()
  85 {
  86
  87   // Pick up siscone header
  88   AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader;
  89   Bool_t debug  = header->GetDebug();     // debug option
  90   Int_t fOpt    = fReader->GetReaderHeader()->GetDetector();
  91
  92   // Check if we are reading AOD jets
  93   TRefArray *refs = 0;
  94   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
  95   if (fromAod) { refs = fReader->GetReferences(); }
  96
  97
  98 //******************************** SISCONE PLUGIN CONFIGURATION
  99 // Here we look for SISCone parameters in the header and we define our plugin.
 100
 101   Double_t coneRadius                = header->GetConeRadius();                 // cone radius
 102   Double_t overlapThreshold          = header->GetOverlapThreshold();           // overlap parameter
 103   Int_t    nPassMax                  = header->GetNPassMax();                   // maximum number of passes
 104   Double_t ptProtoJetMin             = header->GetPtProtojetMin();              // pT min of protojets
 105   Double_t caching                   = header->GetCaching();                    // do we record found cones for this set of data?
 106
 107 //  if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
 108 //  Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
 109
 110   fastjet::JetDefinition::Plugin * plugin;
 111   plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
 112
 113   vector<fastjet::PseudoJet> inputParticles;
 114 //********************************  READING OF INPUT PARTICLES
 115 // Here we look for px, py pz and energy of each particle that we gather in a PseudoJet object, and we put all these PseudoJet in a vector of PseudoJets : input_particles.
 116
 117     if(fOpt==0)
 118       {
 119         TClonesArray *lvArray = fReader->GetMomentumArray();
 120         Int_t nIn = lvArray->GetEntries();
 121
 122         // We check if lvArray is ok
 123         if(lvArray == 0)
 124           {
 125             cout << "Could not get the momentum array" << endl;
 126             return;
 127           }
 128
 129         if(nIn == 0)// nIn = Number of particles in the event
 130           {
 131             if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
 132             return;
 133           }
 134
 135         fJets->SetNinput(nIn) ; // fJets = AliJet number of input objects
 136         Float_t px,py,pz,en;
 137
 138         // Load input vectors
 139         for(Int_t i = 0; i < nIn; i++)
 140           {
 141             TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
 142             px = lv->Px();
 143             py = lv->Py();
 144             pz = lv->Pz();
 145             en = lv->Energy();
 146
 147             fastjet::PseudoJet inputPart(px,py,pz,en);
 148             inputPart.set_user_index(i);
 149             inputParticles.push_back(inputPart);
 150
 151           }
 152       }
 153     else {
 154       TClonesArray* fUnit = fReader->GetUnitArray();
 155       if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
 156       Int_t         nCandidate = fReader->GetNumCandidate();
 157       Int_t         nIn = fUnit->GetEntries();
 158       if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
 159       fJets->SetNinput(nCandidate); // number of input objects // ME
 160       // Information extracted from fUnitArray
 161       // load input vectors and calculate total energy in array
 162       Float_t pt,eta,phi,theta,px,py,pz,en;
 163       Int_t ipart = 0;
 164       for(Int_t i=0; i<nIn; i++)
 165         {
 166           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
 167
 168           if(uArray->GetUnitEnergy()>0.){
 169
 170             // It is not necessary anymore to cut on particle pt
 171             pt    = uArray->GetUnitEnergy();
 172             eta   = uArray->GetUnitEta();
 173             phi   = uArray->GetUnitPhi();
 174             theta = EtaToTheta(eta);
 175             en    = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
 176             px    = TMath::Cos(phi)*pt;
 177             py    = TMath::Sin(phi)*pt;
 178             pz    = en*TMath::TanH(eta);
 179             if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
 180
 181             fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
 182             inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
 183             inputParticles.push_back(inputPart);  // back of the input_particles vector
 184             ipart++;
 185           }
 186         } // End loop on UnitArray
 187     }
 188
 189 //******************************** CHOICE OF JET AREA
 190 // Here we determine jets area for subtracting background later
 191 // For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
 192
 193   Double_t ghostEtamax        = header->GetGhostEtaMax();       // maximum eta in which a ghost can be generated
 194   Double_t ghostArea          = header->GetGhostArea();         // area of a ghost
 195   Int_t    activeAreaRepeats  = header->GetActiveAreaRepeats(); // do we repeat area calculation?
 196   Double_t gridScatter        = header->GetGridScatter();       // fractional random fluctuations of the position of the ghosts on the y-phi grid
 197   Double_t ktScatter          = header->GetKtScatter();         // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
 198   Double_t meanGhostKt        = header->GetMeanGhostKt();       // average transverse momentum of the ghosts.
 199
 200   Double_t areaTypeNumber  = header->GetAreaTypeNumber();       // the number determines jet area type
 201   fastjet::AreaType areaType = fastjet::active_area;
 202   if (areaTypeNumber == 1)  areaType = fastjet::active_area;
 203   if (areaTypeNumber == 2)  areaType = fastjet::active_area_explicit_ghosts;
 204   if (areaTypeNumber == 3)  areaType = fastjet::one_ghost_passive_area;
 205   if (areaTypeNumber == 4)  areaType = fastjet::passive_area;
 206   if (areaTypeNumber == 5)  areaType = fastjet::voronoi_area;
 207
 208   fastjet::AreaDefinition areaDef;
 209
 210   if (areaTypeNumber < 5)
 211   {
 212     fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
 213     areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 214   }
 215
 216   if (areaTypeNumber == 5)
 217   {
 218     Double_t effectiveRFact = header->GetEffectiveRFact();
 219     fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
 220     areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 221   }
 222
 223 //********************************
 224 //********************************
 225
 226   Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
 227
 228 //********************************
 229 //********************************
 230
 231 //********************************
 232 //******************************** BG SUBTRACTION
 233   if (bgMode == 1)// BG subtraction
 234   {
 235     //******************************** JETS FINDING AND EXTRACTION
 236     fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
 237     // Here we extract inclusive jets with pt > ptmin, sorted by pt
 238     Double_t ptMin = header->GetMinJetPt();
 239     vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
 240     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
 241
 242     //***************************** BACKGROUND SUBTRACTION
 243
 244     // Set the rapidity-azimuth range within which to study background
 245     Double_t rapMin = header->GetRapMin();
 246     Double_t rapMax = header->GetRapMax();
 247     Double_t phiMin = header->GetPhiMin();
 248     Double_t phiMax = header->GetPhiMax();
 249     fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
 250
 251     // rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas)
 252     Int_t algo = header->GetBGAlgorithm();
 253     fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
 254     if (algo == 0) algorithm = fastjet::kt_algorithm;
 255     if (algo == 1) algorithm = fastjet::cambridge_algorithm;
 256     fastjet::JetDefinition jetDefForRho(algorithm, 0.5);
 257     fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef);
 258     Double_t rho = csForRho.median_pt_per_unit_area_4vector(range);
 259     cout<<"rho = "<<rho<<endl;
 260
 261     // Vector of corrected jets
 262     vector<fastjet::PseudoJet> corrJets = clustSeq.subtracted_jets(rho, ptMin);
 263
 264     //***************************** JETS DISPLAY
 265
 266     for (size_t j = 0; j < jets.size(); j++)
 267     {
 268       // If the jet is only made of ghosts, continue.
 269       if (clustSeq.is_pure_ghost(jets[j]) == 1) continue;
 270
 271       // If the correction is > jet energy px = py = pz = e = 0
 272       if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;
 273
 274       cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
 275       cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
 276 //     cout<<"px = "<<jets[j].px()<<endl;
 277 //     cout<<"py = "<<jets[j].py()<<endl;
 278 //     cout<<"pz = "<<jets[j].pz()<<endl;
 279       cout<<"e = "<<jets[j].E()<<endl;
 280
 281       cout<<"********************************** Corrected jet(s)"<<endl;
 282       cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
 283 //     cout<<"px = "<<corrJets[j].px()<<endl;
 284 //     cout<<"py = "<<corrJets[j].py()<<endl;
 285 //     cout<<"pz = "<<corrJets[j].pz()<<endl;
 286       cout<<"e = "<<corrJets[j].E()<<endl;
 287
 288       // Go to write AOD info
 289       AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
 290       if(debug) aodjet.Print("");
 291       AddJet(aodjet);
 292     }
 293   }
 294
 295 //********************************
 296 //******************************** NO BG SUBTRACTION
 297
 298   if (bgMode == 0)// No BG subtraction
 299   {
 300     //******************************** JETS FINDING AND EXTRACTION
 301     fastjet::ClusterSequence clustSeq(inputParticles, plugin);
 302     // Here we extract inclusive jets with pt > ptmin, sorted by pt
 303     Double_t ptMin = header->GetMinJetPt();
 304     vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
 305     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
 306
 307     //***************************** JETS DISPLAY
 308
 309     for (size_t k = 0; k < jets.size(); k++)
 310     {
 311       cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
 312       cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
 313 //     cout<<"px = "<<jets[k].px()<<endl;
 314 //     cout<<"py = "<<jets[k].py()<<endl;
 315 //     cout<<"pz = "<<jets[k].pz()<<endl;
 316       cout<<"e = "<<jets[k].E()<<endl;
 317
 318       // Go to write AOD info
 319       AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
 320       if(debug) aodjet.Print("");
 321       AddJet(aodjet);
 322     }
 323   }
 324 }
 325
 326 //____________________________________________________________________________
 327
 328 Float_t  AliSISConeJetFinder::EtaToTheta(Float_t arg)
 329 {
 330   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
 331   return 2.*atan(exp(-arg));
 332
 333
 334 }
 335
 336 //____________________________________________________________________________
 337
 338 void AliSISConeJetFinder::InitTask(TChain *tree)
 339 {
 340
 341   printf("SISCone jet finder initialization ******************");
 342   fReader->CreateTasks(tree);
 343
 344 }
 345
 346
 347 //____________________________________________________________________________
 348
 349 void AliSISConeJetFinder::WriteJHeaderToFile()
 350 {
 351   fHeader->Write();
 352 }
 353
 354 //____________________________________________________________________________
 355