Keep track of missing DCS points in DDL maps (flagged by 'x')
[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   Int_t debug  = header->GetDebug();     // debug option
89   Int_t fOpt    = fReader->GetReaderHeader()->GetDetector();
90
91   // Check if we are reading AOD jets
92   TRefArray *refs = 0;
93   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
94   if (fromAod) { refs = fReader->GetReferences(); }
95   
96
97 //******************************** SISCONE PLUGIN CONFIGURATION
98 // Here we look for SISCone parameters in the header and we define our plugin.  
99
100   Double_t coneRadius                = header->GetConeRadius();                 // cone radius
101   Double_t overlapThreshold          = header->GetOverlapThreshold();           // overlap parameter
102   Int_t    nPassMax                  = header->GetNPassMax();                   // maximum number of passes
103   Double_t ptProtoJetMin             = header->GetPtProtojetMin();              // pT min of protojets
104   Double_t caching                   = header->GetCaching();                    // do we record found cones for this set of data?
105
106 //  if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
107 //  Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
108
109   fastjet::JetDefinition::Plugin * plugin;
110   plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
111
112   vector<fastjet::PseudoJet> inputParticles; 
113 //********************************  READING OF INPUT PARTICLES
114 // 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. 
115
116     if(fOpt==0) 
117       {
118         TClonesArray *lvArray = fReader->GetMomentumArray();
119         Int_t nIn = lvArray->GetEntries();
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         if(nIn == 0)// nIn = Number of particles in the event
130           {
131             if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
132             delete plugin;
133             return;
134           }
135         
136         Float_t px,py,pz,en;
137         
138         // Load input vectors
139         for(Int_t i = 0; i < nIn; i++)
140           { 
141             TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
142             px = lv->Px();
143             py = lv->Py();
144             pz = lv->Pz();
145             en = lv->Energy();
146             
147             fastjet::PseudoJet inputPart(px,py,pz,en); 
148             inputPart.set_user_index(i);
149             inputParticles.push_back(inputPart); 
150             
151           }
152       }
153     else {
154       TClonesArray* fUnit = fReader->GetUnitArray();
155       if(fUnit == 0) { cout << "Could not get the momentum array" << endl;          delete plugin; return; }
156       Int_t         nIn = fUnit->GetEntries();
157       if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;           delete plugin; return; }
158         // Information extracted from fUnitArray
159       // load input vectors and calculate total energy in array
160       Float_t pt,eta,phi,theta,px,py,pz,en;
161       Int_t ipart = 0;
162       for(Int_t i=0; i<nIn; i++) 
163         {
164           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
165           
166           if(uArray->GetUnitEnergy()>0.){
167             
168             // It is not necessary anymore to cut on particle pt
169             pt    = uArray->GetUnitEnergy();
170             eta   = uArray->GetUnitEta();
171             phi   = uArray->GetUnitPhi();
172             theta = EtaToTheta(eta);
173             en    = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
174             px    = TMath::Cos(phi)*pt;
175             py    = TMath::Sin(phi)*pt;
176             pz    = en*TMath::TanH(eta);
177             if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
178             
179             fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
180             inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
181             inputParticles.push_back(inputPart);  // back of the input_particles vector 
182             ipart++;
183           }
184         } // End loop on UnitArray 
185     }
186   
187 //******************************** CHOICE OF JET AREA  
188 // Here we determine jets area for subtracting background later
189 // For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
190
191   Double_t ghostEtamax        = header->GetGhostEtaMax();       // maximum eta in which a ghost can be generated
192   Double_t ghostArea          = header->GetGhostArea();         // area of a ghost
193   Int_t    activeAreaRepeats  = header->GetActiveAreaRepeats(); // do we repeat area calculation?
194   Double_t gridScatter        = header->GetGridScatter();       // fractional random fluctuations of the position of the ghosts on the y-phi grid
195   Double_t ktScatter          = header->GetKtScatter();         // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
196   Double_t meanGhostKt        = header->GetMeanGhostKt();       // average transverse momentum of the ghosts.
197
198   Double_t areaTypeNumber  = header->GetAreaTypeNumber();       // the number determines jet area type 
199   fastjet::AreaType areaType = fastjet::active_area;
200   if (areaTypeNumber == 1)  areaType = fastjet::active_area;                 
201   if (areaTypeNumber == 2)  areaType = fastjet::active_area_explicit_ghosts; 
202   if (areaTypeNumber == 3)  areaType = fastjet::one_ghost_passive_area;      
203   if (areaTypeNumber == 4)  areaType = fastjet::passive_area;                
204   if (areaTypeNumber == 5)  areaType = fastjet::voronoi_area;                
205
206   fastjet::AreaDefinition areaDef;
207
208   if (areaTypeNumber < 5) 
209   {
210     fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
211     areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
212   }
213
214   if (areaTypeNumber == 5)
215   {
216     Double_t effectiveRFact = header->GetEffectiveRFact();
217     fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
218     areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
219   }
220
221 //********************************
222 //********************************
223
224   Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
225
226 //********************************
227 //********************************
228
229 //********************************
230 //******************************** BG SUBTRACTION
231   if (bgMode == 1)// BG subtraction
232   {
233     //******************************** JETS FINDING AND EXTRACTION
234     fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
235     // Here we extract inclusive jets with pt > ptmin, sorted by pt 
236     Double_t ptMin = header->GetMinJetPt(); 
237     vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
238     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
239   
240     //***************************** BACKGROUND SUBTRACTION
241
242     // Set the rapidity-azimuth range within which to study background
243     Double_t rapMin = header->GetRapMin();
244     Double_t rapMax = header->GetRapMax();
245     Double_t phiMin = header->GetPhiMin();
246     Double_t phiMax = header->GetPhiMax();
247     fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
248
249     // rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas)
250     Int_t algo = header->GetBGAlgorithm();
251     fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
252     if (algo == 0) algorithm = fastjet::kt_algorithm;
253     if (algo == 1) algorithm = fastjet::cambridge_algorithm;
254     fastjet::JetDefinition jetDefForRho(algorithm, 0.5);
255     fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef);
256     Double_t rho = csForRho.median_pt_per_unit_area_4vector(range);
257     cout<<"rho = "<<rho<<endl;
258
259     // Vector of corrected jets
260     vector<fastjet::PseudoJet> corrJets = clustSeq.subtracted_jets(rho, ptMin);
261   
262     //***************************** JETS DISPLAY
263
264     for (size_t j = 0; j < jets.size(); j++)
265     {
266       // If the jet is only made of ghosts, continue.
267       if (clustSeq.is_pure_ghost(jets[j]) == 1) continue;
268
269       // If the correction is > jet energy px = py = pz = e = 0
270       if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;
271
272       if(debug){
273         cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
274         cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
275         //     cout<<"px = "<<jets[j].px()<<endl;
276         //     cout<<"py = "<<jets[j].py()<<endl;
277         //     cout<<"pz = "<<jets[j].pz()<<endl;
278         cout<<"e = "<<jets[j].E()<<endl;
279         
280         cout<<"********************************** Corrected jet(s)"<<endl;
281         cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
282         //     cout<<"px = "<<corrJets[j].px()<<endl;
283         //     cout<<"py = "<<corrJets[j].py()<<endl;
284         //     cout<<"pz = "<<corrJets[j].pz()<<endl;
285         cout<<"e = "<<corrJets[j].E()<<endl;
286       }
287       // Go to write AOD info
288       AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
289       if(debug) aodjet.Print("");
290       AddJet(aodjet);
291     }
292   }
293
294 //********************************
295 //******************************** NO BG SUBTRACTION
296
297   if (bgMode == 0)// No BG subtraction
298   {
299     //******************************** JETS FINDING AND EXTRACTION
300     fastjet::ClusterSequence clustSeq(inputParticles, plugin);
301     // Here we extract inclusive jets with pt > ptmin, sorted by pt 
302     Double_t ptMin = header->GetMinJetPt(); 
303     vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
304     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
305   
306     //***************************** JETS DISPLAY
307
308     for (size_t k = 0; k < jets.size(); k++)
309     {
310       if(debug){
311         cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
312         cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
313         //     cout<<"px = "<<jets[k].px()<<endl;
314         //     cout<<"py = "<<jets[k].py()<<endl;
315         //     cout<<"pz = "<<jets[k].pz()<<endl;
316         cout<<"e = "<<jets[k].E()<<endl;
317       }
318       // Go to write AOD info
319       AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
320       if(debug) aodjet.Print("");
321       AddJet(aodjet);
322     }
323   }
324
325   delete plugin;
326
327 }
328  
329 //____________________________________________________________________________
330
331 Float_t  AliSISConeJetFinder::EtaToTheta(Float_t arg)
332 {
333   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
334   return 2.*atan(exp(-arg));
335
336
337 }
338
339 //____________________________________________________________________________
340
341 void AliSISConeJetFinder::InitTask(TChain *tree)
342 {
343
344   printf("SISCone jet finder initialization ******************");
345   fReader->CreateTasks(tree);
346
347 }
348
349
350 //____________________________________________________________________________
351
352 void AliSISConeJetFinder::WriteJHeaderToFile() const
353 {
354   fHeader->Write();
355 }
356
357 //____________________________________________________________________________
358