adding JETAN and FASTJETAN development libs for new i/o of tracks/particles for the...
[u/mrichter/AliRoot.git] / JETAN / DEV / 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 /* $Id$ */
17
18 //---------------------------------------------------------------------
19 // FastJet v2.3.4 finder algorithm interface
20 // Last modification: Neutral cell energy included in the jet reconstruction
21 //
22 // Authors: Rafael.Diaz.Valdes@cern.ch
23 //          Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option)
24 //
25 // ** 2011 magali.estienne@subatech.in2p3.fr &  alexandre.shabetai@cern.ch
26 // new implementation of background subtraction
27 // allowing to subtract bkg using a different algo than the one used for signal jets
28 //---------------------------------------------------------------------
29
30 #include <Riostream.h>
31
32 #include "AliFastJetFinder.h"
33 #include "AliFastJetHeaderV1.h"
34 #include "AliFastJetInput.h"
35 #include "AliFastJetBkg.h"
36 #include "AliAODJetEventBackground.h"
37 #include "AliAODJet.h"
38
39 #include "fastjet/PseudoJet.hh"
40 #include "fastjet/ClusterSequenceArea.hh"
41 #include "fastjet/AreaDefinition.hh"
42 #include "fastjet/JetDefinition.hh"
43
44 #include<vector> 
45
46 using namespace std;
47
48 ClassImp(AliFastJetFinder)
49
50 ////////////////////////////////////////////////////////////////////////
51
52 AliFastJetFinder::AliFastJetFinder():
53   AliJetFinder(),
54   fInputFJ(new AliFastJetInput()),
55   fJetBkg(new  AliFastJetBkg())
56 {
57   // Constructor
58 }
59
60 //____________________________________________________________________________
61 AliFastJetFinder::~AliFastJetFinder()
62 {
63   // destructor
64   delete  fInputFJ;
65   delete  fJetBkg;
66
67 }
68
69 //______________________________________________________________________________
70 void AliFastJetFinder::FindJets()
71 {
72   // runs a FASTJET based algo
73
74   //pick up fastjet header
75   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
76   Int_t  debug  = header->GetDebug();     // debug option
77   Bool_t bgMode = header->GetBGMode();    // choose to subtract BG or not
78   if(debug>0) cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
79
80   // RUN ALGORITHM  
81   // read input particles -----------------------------
82
83   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
84   if(inputParticles.size()==0){
85     if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
86     return;
87   }
88
89   // create an object that represents your choice of jet algorithm, and 
90   // the associated parameters
91   double rParam = header->GetRparam();
92   double rBkgParam = header->GetRparamBkg();
93   fastjet::Strategy strategy = header->GetStrategy();
94   fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
95   fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); 
96   fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
97
98   // create an object that specifies how we to define the area
99   fastjet::AreaDefinition areaDef;
100   double ghostEtamax       = header->GetGhostEtaMax(); 
101   double ghostArea         = header->GetGhostArea(); 
102   int    activeAreaRepeats = header->GetActiveAreaRepeats(); 
103   
104   // now create the object that holds info about ghosts
105   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
106   // and from that get an area definition
107   fastjet::AreaType areaType = header->GetAreaType();
108   areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
109   
110   //***************************** JETS FINDING
111   // run the jet clustering with the above jet definition
112   fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
113
114   vector<fastjet::PseudoJet> jets;
115
116   if(bgMode) // Do BG subtraction directly with the same algorithm (cambridge or kt) for jet signal and background
117     {
118       //***************************** JETS FINDING FOR RHO ESTIMATION
119       // run the jet clustering with the above jet definition
120       fastjet::JetAlgorithm algorithmBkg = header->GetBGAlgorithm();
121       fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy);
122       fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef);
123
124       // save a comment in the header
125       TString comment = "Running FastJet algorithm with the following setup. ";
126       comment+= "Jet definition: ";
127       comment+= TString(jetDef.description());
128       comment+= "Jet bckg definition: ";
129       comment+= TString(jetDefBkg.description());
130       comment+= ". Area definition: ";
131       comment+= TString(areaDef.description());
132       comment+= ". Strategy adopted by FastJet: ";
133       comment+= TString(clust_seq.strategy_string());
134       header->SetComment(comment);
135       if(debug>0){
136         cout << "--------------------------------------------------------" << endl;
137         cout << comment << endl;
138         cout << "--------------------------------------------------------" << endl;
139       }
140       
141       // extract the inclusive jets sorted by pt
142       double ptmin = header->GetPtMin(); 
143       vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets();
144       
145       //subtract background // ===========================================
146       // set the rapididty , phi range within which to study the background 
147       double rapMax = header->GetRapMax(); 
148       double rapMin = header->GetRapMin();
149       double phiMax = header->GetPhiMax();
150       double phiMin = header->GetPhiMin();
151       fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
152      
153       // Extract rho and sigma
154       Double_t rho = 0.;
155       Double_t sigma = 0.;
156       Double_t meanarea = 0.;
157       Bool_t kUse4VectorArea = header->Use4VectorArea();
158       vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets();
159       clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, kUse4VectorArea, rho, sigma, meanarea, false);
160
161       // subtract background and extract jets bkg subtracted
162       vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptmin);
163       
164       // sort jets into increasing pt
165       jets = sorted_by_pt(subJets);  
166
167     }
168
169   else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
170     
171     // save a comment in the header
172     TString comment = "Running FastJet algorithm with the following setup. ";
173     comment+= "Jet definition: ";
174     comment+= TString(jetDef.description());
175     comment+= ". Strategy adopted by FastJet: ";
176     comment+= TString(clust_seq.strategy_string());
177     header->SetComment(comment);
178     if(debug>0){
179       cout << "--------------------------------------------------------" << endl;
180       cout << comment << endl;
181       cout << "--------------------------------------------------------" << endl;
182     }
183   
184     // extract the inclusive jets with pt > ptmin, sorted by pt
185     double ptmin = header->GetPtMin(); 
186     vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
187     
188     jets = sorted_by_pt(inclusiveJets); 
189   
190   }
191    
192   for (size_t j = 0; j < jets.size(); j++) { // loop for jets     
193
194     double area      = clust_seq.area(jets[j]);
195     double areaError = clust_seq.area_error(jets[j]);
196
197     if(debug>0) printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp());
198
199     vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
200     int nCon= constituents.size();
201     TArrayI ind(nCon);
202       
203     if ((jets[j].eta() > (header->GetJetEtaMax())) ||
204         (jets[j].eta() < (header->GetJetEtaMin())) ||
205         (jets[j].phi() > (header->GetJetPhiMax())) ||
206         (jets[j].phi() < (header->GetJetPhiMin())) ||
207         (jets[j].perp() < header->GetPtMin())) continue; // acceptance eta range and etmin
208
209     // go to write AOD  info
210     AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
211     aodjet.SetEffArea(area,areaError);
212     //cout << "Printing jet " << endl;
213     if(debug>0) aodjet.Print("");
214       
215     for (int i=0; i < nCon; i++)
216       {
217         fastjet::PseudoJet mPart=constituents[i];
218         ind[i]=mPart.user_index();
219
220         // Jet constituents (charged tracks) added to the AliAODJet
221         AliJetCalTrkEvent* calEvt  = GetCalTrkEvent();
222         for(Int_t itrack=0; itrack<calEvt->GetNCalTrkTracks(); itrack++)
223           {
224             if(itrack==ind[i])
225               {
226                 TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject();
227                 aodjet.AddTrack(track);
228               }
229           }
230       } // End loop on Constituents
231
232     AddJet(aodjet);
233            
234   } // end loop for jets
235
236 }
237
238 //____________________________________________________________________________
239 void AliFastJetFinder::RunTest(const char* datafile)
240 {
241   // This simple test run the kt algorithm for an ascii file testdata.dat
242
243   // read input particles -----------------------------
244   vector<fastjet::PseudoJet> inputParticles;
245   Float_t px,py,pz,en;
246   ifstream in;
247   Int_t nlines = 0;
248   // we assume a file basic.dat in the current directory
249   // this file has 3 columns of float data
250   in.open(datafile);
251   while (1) {
252     in >> px >> py >> pz >> en;
253     if (!in.good()) break;
254     //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz);
255     nlines++;
256     inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en)); 
257   }
258   //printf(" found %d pointsn",nlines);
259   in.close();
260   //////////////////////////////////////////////////
261  
262   // create an object that represents your choice of jet algorithm, and 
263   // the associated parameters
264   double rParam = 1.0;
265   fastjet::Strategy strategy = fastjet::Best;
266   fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme;
267   fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy);
268   
269   // create an object that specifies how we to define the area
270   fastjet::AreaDefinition areaDef;
271   double ghostEtamax = 7.0;
272   double ghostArea    = 0.05;
273   int    activeAreaRepeats = 1;
274   
275   // now create the object that holds info about ghosts
276   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
277   // and from that get an area definition
278   areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec);
279
280   // run the jet clustering with the above jet definition
281   fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
282   
283   // tell the user what was done
284   cout << "--------------------------------------------------------" << endl;
285   cout << "Jet definition was: " << jetDef.description() << endl;
286   cout << "Area definition was: " << areaDef.description() << endl;
287   cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl;
288   cout << "--------------------------------------------------------" << endl;
289   
290   // extract the inclusive jets with pt > 5 GeV, sorted by pt
291   double ptmin = 5.0;
292   vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
293   
294   cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
295  
296   //subtract background // ===========================================
297   // set the rapididty range within which to study the background 
298   double rapMax = ghostEtamax - rParam;
299   fastjet::RangeDefinition range(rapMax);
300   // subtract background
301   vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);  
302   
303   // print them out //================================================
304   cout << "Printing inclusive jets  after background subtraction \n";
305   cout << "------------------------------------------------------\n";
306   // sort jets into increasing pt
307   vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);  
308
309   printf(" ijet   rap      phi        Pt         area  +-   err\n");
310   for (size_t j = 0; j < jets.size(); j++) {
311     double area      = clust_seq.area(jets[j]);
312     double areaError = clust_seq.area_error(jets[j]);
313     printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(),
314            jets[j].phi(),jets[j].perp(), area, areaError);
315   }
316   cout << endl;
317   // ================================================================
318
319 }
320
321 //____________________________________________________________________________
322 void AliFastJetFinder::WriteJHeaderToFile() const
323 {
324   // Write Jet Header to file
325   fHeader->Write();
326
327 }
328
329 //____________________________________________________________________________
330 Bool_t AliFastJetFinder::ProcessEvent()
331 {
332   // Process one event
333   // Charged only or charged+neutral jets
334
335   fInputFJ->SetHeader(fHeader);
336   fInputFJ->SetCalTrkEvent(GetCalTrkEvent());
337   fInputFJ->FillInput();
338   
339   // Find Jets
340   FindJets(); 
341
342   // Background
343   if(fAODEvBkg){
344     AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
345
346     fJetBkg->SetHeader(fHeader);
347     fJetBkg->SetFastJetInput(fInputFJ);
348     
349      Int_t count = 0;  
350      if(header->GetBkgFastJetb()){
351         Double_t sigma1 = 0 , meanarea1= 0, bkg1 = 0;
352         fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
353         fAODEvBkg->SetBackground(count,bkg1,sigma1,meanarea1);
354         count++;
355      }
356
357      if(header->GetBkgFastJetWoHardest()){
358         Double_t sigma2 = 0 , meanarea2 = 0, bkg2 = 0;
359         fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
360         fAODEvBkg->SetBackground(count,bkg2,sigma2,meanarea2);
361         count++;
362      }
363
364   }  
365
366   Reset();  
367   return kTRUE;
368
369 }