adding ReverseAlignment function, calculation of adjugate rotation; deleting implemen...
[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;
cc6a2227 91 Int_t debug = header->GetDebug(); // debug option
8838ab7a 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
8ddcf605 96 /* NOT used
392c9b47 97 TRefArray *refs = 0;
98 Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
99 if (fromAod) { refs = fReader->GetReferences(); }
8ddcf605 100 */
392c9b47 101 // RUN ALGORITHM
102 // read input particles -----------------------------
d1993270 103
104 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
03ecc336 105 if(inputParticles.size()==0){
0f5f5d71 106 if(debug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
03ecc336 107 return;
108 }
d1993270 109
392c9b47 110 // create an object that represents your choice of jet algorithm, and
111 // the associated parameters
3dda898d 112 double rParam = header->GetRparam();
392c9b47 113 fastjet::Strategy strategy = header->GetStrategy();
3dda898d 114 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
392c9b47 115 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
d1993270 116 fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
8838ab7a 117
392c9b47 118 // create an object that specifies how we to define the area
3dda898d 119 fastjet::AreaDefinition areaDef;
120 double ghostEtamax = header->GetGhostEtaMax();
121 double ghostArea = header->GetGhostArea();
122 int activeAreaRepeats = header->GetActiveAreaRepeats();
392c9b47 123
124 // now create the object that holds info about ghosts
d1993270 125 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
392c9b47 126 // and from that get an area definition
3dda898d 127 fastjet::AreaType areaType = header->GetAreaType();
d1993270 128 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
392c9b47 129
a38beabd 130
131
132
133
134 if(bgMode) // BG subtraction!!!!!!!!!! Not entering here
8838ab7a 135 {
136 //***************************** JETS FINDING AND EXTRACTION
137 // run the jet clustering with the above jet definition
d1993270 138 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
392c9b47 139
8838ab7a 140 // save a comment in the header
141
142 TString comment = "Running FastJet algorithm with the following setup. ";
143 comment+= "Jet definition: ";
d1993270 144 comment+= TString(jetDef.description());
8838ab7a 145 comment+= ". Area definition: ";
3dda898d 146 comment+= TString(areaDef.description());
8838ab7a 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();
3dda898d 160 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
8838ab7a 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
3dda898d 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);
8838ab7a 172
173 // subtract background
3dda898d 174 vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(range,ptmin);
8838ab7a 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
3dda898d 183 vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
8838ab7a 184 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
185
186 double area = clust_seq.area(jets[j]);
3dda898d 187 double areaError = clust_seq.area_error(jets[j]);
8838ab7a 188
a38beabd 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);
392c9b47 190
191 // go to write AOD info
8838ab7a 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 }
a38beabd 201
202
203
204
205
206
207 else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
d1993270 208
a38beabd 209 TClonesArray* fUnit = fReader->GetUnitArray(); //Big mmentum array
d1993270 210 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
211 Int_t nIn = fUnit->GetEntries();
392c9b47 212
a38beabd 213
d1993270 214 //fastjet::ClusterSequence clust_seq(inputParticles, jetDef);
215 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
216
8838ab7a 217 // save a comment in the header
218
219 TString comment = "Running FastJet algorithm with the following setup. ";
220 comment+= "Jet definition: ";
d1993270 221 comment+= TString(jetDef.description());
8838ab7a 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();
3dda898d 234 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
8838ab7a 235
236 //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
392c9b47 237
3dda898d 238 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); // Added by me
8838ab7a 239 for (size_t j = 0; j < jets.size(); j++) { // loop for jets // Added by me
240
a38beabd 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());
d1993270 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]);
8838ab7a 247 // go to write AOD info
248 AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
d1993270 249 aodjet.SetEffArea(area,0);
8838ab7a 250 //cout << "Printing jet " << endl;
251 if(debug) aodjet.Print("");
a38beabd 252 Int_t count=0;
d1993270 253 for (int i=0; i < nCon; i++)
254 {
255 fastjet::PseudoJet mPart=constituents[i];
256 ind[i]=mPart.user_index();
a38beabd 257 // cout<<i<<" index="<<ind[i]<<endl;
d1993270 258
a38beabd 259 //internal loop over all the unit cells
d1993270 260 Int_t ipart = 0;
a38beabd 261 // count=0;
d1993270 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]){
a38beabd 268 TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef();
269 Int_t tracksInCell = trackArray->GetEntries();
08c28025 270 for(int ji = 0; ji < tracksInCell; ji++){
271 AliAODTrack * track = (AliAODTrack*)trackArray->At(ji);
a38beabd 272 aodjet.AddTrack(track);
273 }
274
275 count++;
d1993270 276 }
277 ipart++;
278 }
279 }
280 }
281
a38beabd 282
d1993270 283
a38beabd 284 AddJet(aodjet);
285
8838ab7a 286 } // end loop for jets
287 }
288
a17e6965 289}
290
392c9b47 291//____________________________________________________________________________
292void AliFastJetFinder::RunTest(const char* datafile)
a17e6965 293
a17e6965 294{
a17e6965 295
392c9b47 296 // This simple test run the kt algorithm for an ascii file testdata.dat
297 // read input particles -----------------------------
3dda898d 298 vector<fastjet::PseudoJet> inputParticles;
392c9b47 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++;
3dda898d 310 inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en));
a17e6965 311 }
392c9b47 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
3dda898d 318 double rParam = 1.0;
392c9b47 319 fastjet::Strategy strategy = fastjet::Best;
3dda898d 320 fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme;
d1993270 321 fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy);
392c9b47 322
323
324
325 // create an object that specifies how we to define the area
3dda898d 326 fastjet::AreaDefinition areaDef;
327 double ghostEtamax = 7.0;
328 double ghostArea = 0.05;
329 int activeAreaRepeats = 1;
392c9b47 330
331
332 // now create the object that holds info about ghosts
d1993270 333 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
392c9b47 334 // and from that get an area definition
d1993270 335 areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec);
392c9b47 336
337
338 // run the jet clustering with the above jet definition
d1993270 339 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
392c9b47 340
341
342 // tell the user what was done
343 cout << "--------------------------------------------------------" << endl;
d1993270 344 cout << "Jet definition was: " << jetDef.description() << endl;
3dda898d 345 cout << "Area definition was: " << areaDef.description() << endl;
392c9b47 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;
3dda898d 352 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
392c9b47 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
3dda898d 359 double rapMax = ghostEtamax - rParam;
360 fastjet::RangeDefinition range(rapMax);
392c9b47 361 // subtract background
3dda898d 362 vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(range,ptmin);
392c9b47 363
364 // print them out //================================================
365 cout << "Printing inclusive jets after background subtraction \n";
366 cout << "------------------------------------------------------\n";
367 // sort jets into increasing pt
3dda898d 368 vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
392c9b47 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]);
3dda898d 374 double areaError = clust_seq.area_error(jets[j]);
392c9b47 375
8838ab7a 376 printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(),
3dda898d 377 jets[j].phi(),jets[j].perp(), area, areaError);
392c9b47 378 }
379 cout << endl;
380 // ================================================================
a17e6965 381
392c9b47 382
383
a17e6965 384}
385
392c9b47 386//____________________________________________________________________________
a17e6965 387
b430dff0 388void AliFastJetFinder::WriteJHeaderToFile() const
a17e6965 389{
a17e6965 390 fHeader->Write();
391}
392
9157b9c9 393//____________________________________________________________________________
8838ab7a 394
395Float_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
405void AliFastJetFinder::InitTask(TChain *tree)
406{
407
408 printf("Fast jet finder initialization ******************");
409 fReader->CreateTasks(tree);
410
411}
d1993270 412
287697fc 413
414Bool_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
294056ab 429
21ead515 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 }
294056ab 442
443
444
445
287697fc 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
d1993270 463Bool_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
21ead515 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 }
294056ab 496
497
f9b8e0ed 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);
d1993270 502
f9b8e0ed 503// fAODEvBkg->SetBackground(0,bkg1);
504// fAODEvBkg->SetBackground(1,bkg2);
505// fAODEvBkg->SetBackground(2,bkg3);
506// fAODEvBkg->SetBackground(3,bkg4);
d1993270 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}