Fixes to track selection for embedding in cluster task (Bastian), Do no select events...
[u/mrichter/AliRoot.git] / 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 "AliJetKineReader.h"
39 #include "AliJetReader.h"
40 #include "AliJetReaderHeader.h"
41 #include "AliJetUnitArray.h"
42 #include "AliSISConeJetFinder.h"
43 #include "AliSISConeJetHeader.h"
44
45 #include "fastjet/AreaDefinition.hh"
46 #include "fastjet/ClusterSequenceArea.hh"
47 #include "fastjet/JetDefinition.hh"
48 #include "fastjet/PseudoJet.hh"
49
50 // get info on how fastjet was configured
51 #include "fastjet/config.h"
52
53 #ifdef ENABLE_PLUGIN_SISCONE
54 #include "fastjet/SISConePlugin.hh"
55 #endif
56
57 #include<sstream>  // needed for internal io
58 #include<vector> 
59 #include <cmath> 
60
61 using namespace std;
62
63
64 ClassImp(AliSISConeJetFinder)
65
66
67 //____________________________________________________________________________
68
69 AliSISConeJetFinder::AliSISConeJetFinder():
70 AliJetFinder()
71 {
72   // Constructor
73 }
74
75 //____________________________________________________________________________
76
77 AliSISConeJetFinder::~AliSISConeJetFinder()
78 {
79   // destructor
80 }
81
82 //______________________________________________________________________________
83 void AliSISConeJetFinder::FindJets()
84 {
85
86   // Pick up siscone header
87   AliSISConeJetHeader *header      = (AliSISConeJetHeader*)fHeader;
88   AliJetReaderHeader *readerheader = (AliJetReaderHeader*)fReader->GetReaderHeader();
89
90   Int_t debug = header->GetDebug();                        // debug option
91   Int_t fOpt  = fReader->GetReaderHeader()->GetDetector();
92
93
94 //******************************** SISCONE PLUGIN CONFIGURATION
95 // Here we look for SISCone parameters in the header and we define our plugin.  
96
97   Double_t coneRadius                = header->GetConeRadius();       // cone radius
98   Double_t overlapThreshold          = header->GetOverlapThreshold(); // overlap parameter
99   Int_t    nPassMax                  = header->GetNPassMax();         // maximum number of passes
100   Double_t ptProtoJetMin             = header->GetPtProtojetMin();    // pT min of protojets
101   Double_t caching                   = header->GetCaching();          // do we record found cones for this set of data?
102
103 //  if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
104 //  Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
105
106   fastjet::JetDefinition::Plugin * plugin;
107   plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
108
109   vector<fastjet::PseudoJet> inputParticles; 
110
111 //********************************  READING OF INPUT PARTICLES
112 // 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. 
113
114   Double_t pttrackcut = readerheader->GetPtCut(); 
115   if (debug)cout<<"pT track cut for SISCone jet finder = "<<pttrackcut<<" GeV"<<endl;
116
117   if(fOpt==0) 
118   {
119     TClonesArray *lvArray    = fReader->GetMomentumArray();
120
121     // We check if lvArray is ok
122     if(lvArray == 0)
123     {
124       cout << "Could not get the momentum array" << endl;
125       delete plugin;
126       return;
127     }
128         
129     Int_t nIn = lvArray->GetEntries();
130         
131     if(nIn == 0)// nIn = Number of particles in the event
132     {
133       if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
134       delete plugin;
135       return;
136     }
137         
138     Float_t px,py,pz,en;
139         
140     // Load input vectors
141     for(Int_t i = 0; i < nIn; i++)
142     { 
143       TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
144       px = lv->Px();
145       py = lv->Py();
146       pz = lv->Pz();
147       en = lv->Energy();
148
149       if (pttrackcut > TMath::Sqrt(px*px+py*py)) continue; //added by syssy
150             
151       fastjet::PseudoJet inputPart(px,py,pz,en); 
152       inputPart.set_user_index(i);
153       inputParticles.push_back(inputPart); 
154     }
155   }
156   else 
157   {
158     TClonesArray* fUnit = fReader->GetUnitArray();
159     if(fUnit == 0) { cout << "Could not get the momentum array" << endl;            delete plugin; return; }
160     Int_t         nIn = fUnit->GetEntries();
161     if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;             delete plugin; return; }
162     // Information extracted from fUnitArray
163     // load input vectors and calculate total energy in array
164     Float_t pt,eta,phi,theta,px,py,pz,en;
165     Int_t ipart = 0;
166     for(Int_t i=0; i<nIn; i++) 
167     {
168       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
169           
170       if(uArray->GetUnitEnergy()>0.)
171       {
172         // It is not necessary anymore to cut on particle pt
173         pt    = uArray->GetUnitEnergy();
174         eta   = uArray->GetUnitEta();
175         phi   = uArray->GetUnitPhi();
176         theta = EtaToTheta(eta);
177         en    = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
178         px    = TMath::Cos(phi)*pt;
179         py    = TMath::Sin(phi)*pt;
180         pz    = en*TMath::TanH(eta);
181
182         if (pttrackcut > TMath::Sqrt(px*px+py*py)) continue; 
183
184         if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
185             
186         fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
187         inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
188         inputParticles.push_back(inputPart);  // back of the input_particles vector 
189         ipart++;
190       }
191     } // End loop on UnitArray 
192   }
193   
194 //******************************** CHOICE OF JET AREA  
195 // Here we determine jets area for subtracting background later
196 // For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
197
198   Double_t ghostEtamax        = header->GetGhostEtaMax();       // maximum eta in which a ghost can be generated
199   Double_t ghostArea          = header->GetGhostArea();         // area of a ghost
200   Int_t    activeAreaRepeats  = header->GetActiveAreaRepeats(); // do we repeat area calculation?
201   Double_t gridScatter        = header->GetGridScatter();       // fractional random fluctuations of the position of the ghosts on the y-phi grid
202   Double_t ktScatter          = header->GetKtScatter();         // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
203   Double_t meanGhostKt        = header->GetMeanGhostKt();       // average transverse momentum of the ghosts.
204
205   Double_t areaTypeNumber  = header->GetAreaTypeNumber();       // the number determines jet area type 
206   fastjet::AreaType areaType = fastjet::active_area;
207   if (areaTypeNumber == 1)  areaType = fastjet::active_area;                 
208   if (areaTypeNumber == 2)  areaType = fastjet::active_area_explicit_ghosts; 
209   if (areaTypeNumber == 3)  areaType = fastjet::one_ghost_passive_area;      
210   if (areaTypeNumber == 4)  areaType = fastjet::passive_area;                
211   if (areaTypeNumber == 5)  areaType = fastjet::voronoi_area;                
212
213   fastjet::AreaDefinition areaDef;
214
215   if (areaTypeNumber < 5) 
216   {
217     fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
218     areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
219   }
220
221   if (areaTypeNumber == 5)
222   {
223     Double_t effectiveRFact = header->GetEffectiveRFact();
224     fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
225     areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
226   }
227
228 //********************************
229 //********************************
230
231   Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
232
233 //********************************
234 //********************************
235
236 //********************************
237 //******************************** BG SUBTRACTION
238   if (bgMode == 1)// BG subtraction
239     {
240     //******************************** JETS FINDING AND EXTRACTION
241     fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
242
243     // Here we extract inclusive jets with pt > ptmin, sorted by pt 
244     Double_t ptMin = header->GetMinJetPt(); 
245     vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
246     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
247   
248     //***************************** BACKGROUND SUBTRACTION
249
250     // Set the rapidity-azimuth range within which to study background
251     Double_t rapMin = header->GetRapMin();
252     Double_t rapMax = header->GetRapMax();
253     Double_t phiMin = header->GetPhiMin();
254     Double_t phiMax = header->GetPhiMax();
255     fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
256
257     // As SISCone gives too small areas to estimate background,
258     // we have to choose between kT and Cambridge/Aachen algorithms
259     // to estimate mean background rho
260     fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
261     Int_t algo = header->GetBGAlgorithm();
262     if (algo == 0) algorithm = fastjet::kt_algorithm;
263     if (algo == 1) algorithm = fastjet::cambridge_algorithm;
264
265     // Rho estimation :
266     // Radius used to calculate rho, can be different from the one used to reconstruct jets
267     Double_t RRho = header->GetRForRho(); 
268     fastjet::JetDefinition jetDefinitionForRhoEstimation(algorithm, RRho);
269     fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefinitionForRhoEstimation, areaDef);
270     vector<fastjet::PseudoJet> inclusiveJetsForRhoEstimation = sorted_by_pt(csForRho.inclusive_jets());
271
272     // Number of hard jets not to be used to estimate rho
273     Int_t NHardJets = header->GetNumberOfJetsToErase();
274     if (debug) cout<<"Number of jets before subtraction : "<<inclusiveJetsForRhoEstimation.size()<<endl;
275     if (debug) cout<<"Number of jets not to count into rho estimation : "<<NHardJets<<endl;
276     Bool_t rhoOnJet = 0;
277     if (inclusiveJetsForRhoEstimation.size() == 1 && inclusiveJetsForRhoEstimation[0].perp() > 5) rhoOnJet = 1;//if there's only one jet (not bg kind), rho = 0
278     if (inclusiveJetsForRhoEstimation.size() <= (UInt_t)NHardJets && inclusiveJetsForRhoEstimation.size() != 1  && inclusiveJetsForRhoEstimation.size() != 0)
279     {
280       //      inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin(),inclusiveJetsForRhoEstimation.begin()+1);
281       inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin());
282       cout<<"Only the 1st (hardest) jet of the event hasn't been taken into account for rho estimation"<<endl;
283     }
284     if (inclusiveJetsForRhoEstimation.size() > (UInt_t)NHardJets) inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin(), inclusiveJetsForRhoEstimation.begin()+NHardJets);
285     
286     // Estimation of rho and fluctuations sigma
287
288     //1st method
289     Double_t rho      = 0;
290     Double_t sigma    = 0;
291     Double_t meanarea = 0;
292     // 3rd argument   : 1 = use area 4 vector rather than simple area
293     // last argument  : 0 = in case of explicit ghosts use
294     if (inclusiveJetsForRhoEstimation.size() != 0) csForRho.get_median_rho_and_sigma(inclusiveJetsForRhoEstimation, range, 1, rho, sigma, meanarea, 1); // this gives the fluctuations also
295
296     if (rhoOnJet) rho = 0;
297
298     if(debug) cout<<"rho = "<<rho<<endl;
299
300     // Vector of corrected jets
301     vector<fastjet::PseudoJet> corrJets = sorted_by_pt(clustSeq.subtracted_jets(rho, ptMin));
302   
303     //***************************** JETS DISPLAY
304
305     for (size_t j = 0; j < jets.size(); j++)
306     {
307 //       // If the jet is only made of ghosts, continue.
308 //       if (clustSeq.is_pure_ghost(corrJets[j]) == 1) continue;
309
310       // If the correction is > jet energy px = py = pz = e = 0
311       if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;
312
313       if(debug)
314       {
315         cout<<"********************************** Reconstructed jet(s) (not corrected)"<<endl;
316         cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
317         cout<<"e = "<<jets[j].E()<<endl;
318         
319         cout<<"********************************** Corrected jet(s)"<<endl;
320         cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
321         cout<<"e = "<<corrJets[j].E()<<endl;
322       }
323       // Go to write AOD info
324       AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
325       if(debug) aodjet.Print("");
326       AddJet(aodjet);
327     }
328   }
329
330 //********************************
331 //******************************** NO BG SUBTRACTION
332
333   if (bgMode == 0)// No BG subtraction
334   {
335     //******************************** JETS FINDING AND EXTRACTION
336     fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
337     // Here we extract inclusive jets with pt > ptmin, sorted by pt 
338     Double_t ptMin = header->GetMinJetPt(); 
339     vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
340     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
341
342     //***************************** JETS DISPLAY
343
344     for (size_t k = 0; k < jets.size(); k++)
345     {
346       if(debug)
347       {
348         cout<<"********************************** Reconstructed jet(s) (not corrected)"<<endl;
349         cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
350         cout<<"e = "<<jets[k].E()<<endl;
351       }
352
353       // Go to write AOD info
354       Double_t area = clustSeq.area(jets[k]);
355       AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
356       aodjet.SetEffArea(area,0);
357       if(debug) aodjet.Print("");
358       AddJet(aodjet);
359     }
360   }
361
362   delete plugin;
363
364 }
365  
366 //____________________________________________________________________________
367
368 Float_t  AliSISConeJetFinder::EtaToTheta(Float_t arg)
369 {
370   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
371   return 2.*atan(exp(-arg));
372
373
374 }
375
376 //____________________________________________________________________________
377
378 void AliSISConeJetFinder::InitTask(TChain *tree)
379 {
380
381   printf("SISCone jet finder initialization ******************");
382   fReader->CreateTasks(tree);
383
384 }
385
386
387 //____________________________________________________________________________
388
389 void AliSISConeJetFinder::WriteJHeaderToFile() const
390 {
391   fHeader->Write();
392 }
393
394 //____________________________________________________________________________
395