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