JETAN/AliFastJetFinder.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 // Last modification: Neutral cell energy included in the jet reconstruction
  20 //
  21 // Authors: Rafael.Diaz.Valdes@cern.ch
  22 //          Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option)
  23 //
  24 //---------------------------------------------------------------------
  25
  26
  27 #include <Riostream.h>
  28 #include <TLorentzVector.h>
  29 #include <TFile.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TArrayF.h>
  33 #include <TClonesArray.h>
  34
  35 #include "AliFastJetFinder.h"
  36 #include "AliFastJetHeaderV1.h"
  37 #include "AliJetReaderHeader.h"
  38 #include "AliJetReader.h"
  39 #include "AliJetUnitArray.h"
  40 #include "AliFastJetInput.h"
  41 #include "AliJetBkg.h"
  42 #include "AliAODJetEventBackground.h"
  43 #include "AliAODTrack.h"
  44
  45 #include "fastjet/PseudoJet.hh"
  46 #include "fastjet/ClusterSequenceArea.hh"
  47 #include "fastjet/AreaDefinition.hh"
  48 #include "fastjet/JetDefinition.hh"
  49 // get info on how fastjet was configured
  50 #include "fastjet/config.h"
  51
  52 #ifdef ENABLE_PLUGIN_SISCONE
  53 #include "fastjet/SISConePlugin.hh"
  54 #endif
  55
  56 #include<sstream>  // needed for internal io
  57 #include<vector>
  58 #include <cmath>
  59
  60 using namespace std;
  61
  62
  63 ClassImp(AliFastJetFinder)
  64
  65
  66 //____________________________________________________________________________
  67
  68 AliFastJetFinder::AliFastJetFinder():
  69   AliJetFinder(),
  70   fInputFJ(new AliFastJetInput()),
  71   fJetBkg(new  AliJetBkg())
  72 {
  73   // Constructor
  74 }
  75
  76 //____________________________________________________________________________
  77
  78 AliFastJetFinder::~AliFastJetFinder()
  79 {
  80   // destructor
  81   delete  fInputFJ; fInputFJ = 0;
  82   delete  fJetBkg; fJetBkg = 0;
  83 }
  84
  85 //______________________________________________________________________________
  86 void AliFastJetFinder::FindJets()
  87 {
  88
  89   //pick up fastjet header
  90   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  91   Int_t  debug  = header->GetDebug();     // debug option
  92   Bool_t bgMode = header->GetBGMode();    // choose to subtract BG or not
  93   if(debug)cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
  94
  95   // check if we are reading AOD jets
  96   /* NOT used
  97   TRefArray *refs = 0;
  98   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
  99   if (fromAod) { refs = fReader->GetReferences(); }
 100   */
 101   // RUN ALGORITHM
 102   // read input particles -----------------------------
 103
 104   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
 105   if(inputParticles.size()==0){
 106     if(debug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 107     return;
 108   }
 109
 110   // create an object that represents your choice of jet algorithm, and
 111   // the associated parameters
 112   double rParam = header->GetRparam();
 113   fastjet::Strategy strategy = header->GetStrategy();
 114   fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
 115   fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
 116   fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
 117
 118   // create an object that specifies how we to define the area
 119   fastjet::AreaDefinition areaDef;
 120   double ghostEtamax = header->GetGhostEtaMax();
 121   double ghostArea   = header->GetGhostArea();
 122   int    activeAreaRepeats = header->GetActiveAreaRepeats();
 123
 124   // now create the object that holds info about ghosts
 125   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
 126   // and from that get an area definition
 127   fastjet::AreaType areaType = header->GetAreaType();
 128   areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 129
 130
 131
 132
 133
 134   if(bgMode) // BG subtraction!!!!!!!!!! Not entering here
 135     {
 136       //***************************** JETS FINDING AND EXTRACTION
 137       // run the jet clustering with the above jet definition
 138       fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 139
 140       // save a comment in the header
 141
 142       TString comment = "Running FastJet algorithm with the following setup. ";
 143       comment+= "Jet definition: ";
 144       comment+= TString(jetDef.description());
 145       comment+= ". Area definition: ";
 146       comment+= TString(areaDef.description());
 147       comment+= ". Strategy adopted by FastJet: ";
 148       comment+= TString(clust_seq.strategy_string());
 149       header->SetComment(comment);
 150       if(debug){
 151         cout << "--------------------------------------------------------" << endl;
 152         cout << comment << endl;
 153         cout << "--------------------------------------------------------" << endl;
 154       }
 155       //header->PrintParameters();
 156
 157
 158       // extract the inclusive jets with pt > ptmin, sorted by pt
 159       double ptmin = header->GetPtMin();
 160       vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 161
 162       //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 163
 164
 165       //subtract background // ===========================================
 166       // set the rapididty , phi range within which to study the background
 167       double rapMax = header->GetRapMax();
 168       double rapMin = header->GetRapMin();
 169       double phiMax = header->GetPhiMax();
 170       double phiMin = header->GetPhiMin();
 171       fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
 172
 173       // subtract background
 174       vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);
 175
 176       // print out
 177       //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n";
 178       //cout << "---------------------------------------\n";
 179       //cout << endl;
 180       //printf(" ijet   rap      phi        Pt         area  +-   err\n");
 181
 182       // sort jets into increasing pt
 183       vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
 184       for (size_t j = 0; j < jets.size(); j++) { // loop for jets
 185
 186         double area     = clust_seq.area(jets[j]);
 187         double areaError = clust_seq.area_error(jets[j]);
 188
 189         if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError);
 190
 191         // go to write AOD  info
 192         AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
 193         //cout << "Printing jet " << endl;
 194         if(debug) aodjet.Print("");
 195         //cout << "Adding jet ... " ;
 196         AddJet(aodjet);
 197         //cout << "added \n" << endl;
 198
 199       }
 200     }
 201
 202
 203
 204
 205
 206
 207   else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
 208
 209     TClonesArray* fUnit = fReader->GetUnitArray(); //Big mmentum array
 210     if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
 211     Int_t         nIn = fUnit->GetEntries();
 212
 213
 214     //fastjet::ClusterSequence clust_seq(inputParticles, jetDef);
 215     fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 216
 217     // save a comment in the header
 218
 219     TString comment = "Running FastJet algorithm with the following setup. ";
 220     comment+= "Jet definition: ";
 221     comment+= TString(jetDef.description());
 222     comment+= ". Strategy adopted by FastJet: ";
 223     comment+= TString(clust_seq.strategy_string());
 224     header->SetComment(comment);
 225     if(debug){
 226       cout << "--------------------------------------------------------" << endl;
 227       cout << comment << endl;
 228       cout << "--------------------------------------------------------" << endl;
 229     }
 230     //header->PrintParameters();
 231
 232       // extract the inclusive jets with pt > ptmin, sorted by pt
 233     double ptmin = header->GetPtMin();
 234     vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 235
 236     //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 237
 238     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); // Added by me
 239     for (size_t j = 0; j < jets.size(); j++) { // loop for jets     // Added by me
 240
 241       if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp());
 242
 243       vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
 244       int nCon= constituents.size();
 245       TArrayI ind(nCon);
 246       Double_t area=clust_seq.area(jets[j]);
 247       // go to write AOD  info
 248       AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
 249       aodjet.SetEffArea(area,0);
 250       //cout << "Printing jet " << endl;
 251       if(debug) aodjet.Print("");
 252       Int_t count=0;
 253       for (int i=0; i < nCon; i++)
 254         {
 255         fastjet::PseudoJet mPart=constituents[i];
 256         ind[i]=mPart.user_index();
 257         //      cout<<i<<"  index="<<ind[i]<<endl;
 258
 259         //internal loop over all the unit cells
 260         Int_t ipart = 0;
 261         //      count=0;
 262         for(Int_t ii=0; ii<nIn; ii++)
 263           {
 264             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(ii);
 265             if(uArray->GetUnitEnergy()>0.){
 266               uArray->SetUnitTrackID(ipart);//used to have the index available in AliJEtBkg
 267               if(ipart==ind[i]){
 268                 TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef();
 269                 Int_t tracksInCell = trackArray->GetEntries();
 270                 for(int ji = 0; ji < tracksInCell; ji++){
 271                   AliAODTrack * track = (AliAODTrack*)trackArray->At(ji);
 272                   aodjet.AddTrack(track);
 273                 }
 274
 275                 count++;
 276               }
 277               ipart++;
 278             }
 279           }
 280       }
 281
 282
 283
 284       AddJet(aodjet);
 285
 286     } // end loop for jets
 287   }
 288
 289 }
 290
 291 //____________________________________________________________________________
 292 void AliFastJetFinder::RunTest(const char* datafile)
 293
 294 {
 295
 296    // This simple test run the kt algorithm for an ascii file testdata.dat
 297    // read input particles -----------------------------
 298   vector<fastjet::PseudoJet> inputParticles;
 299   Float_t px,py,pz,en;
 300   ifstream in;
 301   Int_t nlines = 0;
 302   // we assume a file basic.dat in the current directory
 303   // this file has 3 columns of float data
 304   in.open(datafile);
 305   while (1) {
 306       in >> px >> py >> pz >> en;
 307       if (!in.good()) break;
 308       //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz);
 309       nlines++;
 310       inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en));
 311    }
 312    //printf(" found %d pointsn",nlines);
 313    in.close();
 314    //////////////////////////////////////////////////
 315
 316   // create an object that represents your choice of jet algorithm, and
 317   // the associated parameters
 318   double rParam = 1.0;
 319   fastjet::Strategy strategy = fastjet::Best;
 320   fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme;
 321   fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy);
 322
 323
 324
 325   // create an object that specifies how we to define the area
 326   fastjet::AreaDefinition areaDef;
 327   double ghostEtamax = 7.0;
 328   double ghostArea    = 0.05;
 329   int    activeAreaRepeats = 1;
 330
 331
 332   // now create the object that holds info about ghosts
 333   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
 334   // and from that get an area definition
 335   areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec);
 336
 337
 338   // run the jet clustering with the above jet definition
 339   fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 340
 341
 342   // tell the user what was done
 343   cout << "--------------------------------------------------------" << endl;
 344   cout << "Jet definition was: " << jetDef.description() << endl;
 345   cout << "Area definition was: " << areaDef.description() << endl;
 346   cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl;
 347   cout << "--------------------------------------------------------" << endl;
 348
 349
 350   // extract the inclusive jets with pt > 5 GeV, sorted by pt
 351   double ptmin = 5.0;
 352   vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 353
 354   cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 355
 356
 357   //subtract background // ===========================================
 358   // set the rapididty range within which to study the background
 359   double rapMax = ghostEtamax - rParam;
 360   fastjet::RangeDefinition range(rapMax);
 361   // subtract background
 362   vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);
 363
 364   // print them out //================================================
 365   cout << "Printing inclusive jets  after background subtraction \n";
 366   cout << "------------------------------------------------------\n";
 367   // sort jets into increasing pt
 368   vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
 369
 370   printf(" ijet   rap      phi        Pt         area  +-   err\n");
 371   for (size_t j = 0; j < jets.size(); j++) {
 372
 373     double area     = clust_seq.area(jets[j]);
 374     double areaError = clust_seq.area_error(jets[j]);
 375
 376     printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(),
 377            jets[j].phi(),jets[j].perp(), area, areaError);
 378   }
 379   cout << endl;
 380   // ================================================================
 381
 382
 383
 384 }
 385
 386 //____________________________________________________________________________
 387
 388 void AliFastJetFinder::WriteJHeaderToFile() const
 389 {
 390   fHeader->Write();
 391 }
 392
 393 //____________________________________________________________________________
 394
 395 Float_t  AliFastJetFinder::EtaToTheta(Float_t arg)
 396 {
 397   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
 398   return 2.*atan(exp(-arg));
 399
 400
 401 }
 402
 403 //____________________________________________________________________________
 404
 405 void AliFastJetFinder::InitTask(TChain *tree)
 406 {
 407
 408   printf("Fast jet finder initialization ******************");
 409   fReader->CreateTasks(tree);
 410
 411 }
 412
 413
 414 Bool_t AliFastJetFinder::ProcessEvent()
 415 {
 416   //
 417   // Process one event
 418   // from meomntum array
 419
 420   Bool_t ok = fReader->FillMomentumArray();
 421
 422   if (!ok) return kFALSE;
 423   fInputFJ->SetHeader(fHeader);
 424   fInputFJ->SetReader(fReader);
 425   fInputFJ->FillInput();
 426   // Jets
 427   FindJets();
 428
 429
 430   if( fAODEvBkg){
 431     fJetBkg->SetHeader(fHeader);
 432     fJetBkg->SetReader(fReader);
 433     Double_t sigma1 = 0,meanarea1= 0,sigma2 = 0,meanarea2 = 0;
 434     Double_t bkg1 = 0,bkg2 = 0;
 435
 436     fJetBkg->SetFastJetInput(fInputFJ);
 437     fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
 438     fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
 439     fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
 440     fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);
 441   }
 442
 443
 444
 445
 446   /*
 447   fJetBkg->SetFastJetInput(fInputFJ);
 448   Double_t bkg1=fJetBkg->BkgFastJet();
 449   Double_t bkg2=fJetBkg->BkgChargedFastJet();
 450   Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets);
 451   Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets);
 452
 453   fAODEvBkg->SetBackground(0,bkg1);
 454   fAODEvBkg->SetBackground(1,bkg2);
 455   fAODEvBkg->SetBackground(2,bkg3);
 456   fAODEvBkg->SetBackground(3,bkg4);
 457   */
 458   Reset();
 459   return kTRUE;
 460
 461 }
 462
 463 Bool_t AliFastJetFinder::ProcessEvent2()
 464 {
 465   //
 466   // Process one event
 467   // Charged only or charged+neutral jets
 468   //
 469
 470   TRefArray* ref = new TRefArray();
 471   Bool_t procid = kFALSE;
 472   Bool_t ok = fReader->ExecTasks(procid,ref);
 473
 474   // Delete reference pointer
 475   if (!ok) {delete ref; return kFALSE;}
 476
 477   // Leading particles
 478   fInputFJ->SetHeader(fHeader);
 479   fInputFJ->SetReader(fReader);
 480   fInputFJ->FillInput();
 481
 482   // Jets
 483   FindJets();
 484
 485   if( fAODEvBkg){
 486     fJetBkg->SetHeader(fHeader);
 487     fJetBkg->SetReader(fReader);
 488     fJetBkg->SetFastJetInput(fInputFJ);
 489     Double_t sigma1,meanarea1,sigma2,meanarea2;
 490     Double_t bkg1,bkg2;
 491     fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
 492     fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
 493     fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
 494     fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);
 495   }
 496
 497
 498 //  Double_t bkg1=fJetBkg->BkgFastJet();
 499 //  Double_t bkg2=fJetBkg->BkgChargedFastJet();
 500 //  Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets);
 501 //  Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets);
 502
 503 //  fAODEvBkg->SetBackground(0,bkg1);
 504 //  fAODEvBkg->SetBackground(1,bkg2);
 505 //  fAODEvBkg->SetBackground(2,bkg3);
 506 //  fAODEvBkg->SetBackground(3,bkg4);
 507
 508   Int_t nEntRef    = ref->GetEntries();
 509
 510   for(Int_t i=0; i<nEntRef; i++)
 511     {
 512       // Reset the UnitArray content which were referenced
 513       ((AliJetUnitArray*)ref->At(i))->SetUnitTrackID(0);
 514       ((AliJetUnitArray*)ref->At(i))->SetUnitEnergy(0.);
 515       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller);
 516       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller);
 517       ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad);
 518       ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlagC(kTRUE,kBad);
 519       ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc);
 520       ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet);
 521       ((AliJetUnitArray*)ref->At(i))->ClearUnitTrackRef();
 522
 523       // Reset process ID
 524       AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i);
 525       uA->ResetBit(kIsReferenced);
 526       uA->SetUniqueID(0);
 527     }
 528
 529   // Delete the reference pointer
 530   ref->Delete();
 531   delete ref;
 532
 533   Reset();
 534
 535   return kTRUE;
 536 }