fix of some rule violations
[u/mrichter/AliRoot.git] / JETAN / AliJetBkg.cxx
CommitLineData
50c7d9f7 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#include <Riostream.h>
17#include <TChain.h>
18#include <TFile.h>
19#include <TList.h>
20#include <TArrayF.h>
21#include <TClonesArray.h>
22#include <TF1.h>
23#include <TString.h>
24
25#include "AliJetHeader.h"
26#include "AliJetReader.h"
27#include "AliJetReaderHeader.h"
28#include "AliFastJetFinder.h"
29#include "AliFastJetHeaderV1.h"
30#include "AliJetReaderHeader.h"
31#include "AliJetReader.h"
32#include "AliJetUnitArray.h"
33#include "AliFastJetInput.h"
34#include "AliESDEvent.h"
35
36
37#include "fastjet/PseudoJet.hh"
38#include "fastjet/ClusterSequenceArea.hh"
39#include "fastjet/AreaDefinition.hh"
40#include "fastjet/JetDefinition.hh"
41// get info on how fastjet was configured
42#include "fastjet/config.h"
43
44#ifdef ENABLE_PLUGIN_SISCONE
45#include "fastjet/SISConePlugin.hh"
46#endif
47
48#include<sstream> // needed for internal io
49#include<vector>
50#include <cmath>
51
52using namespace std;
53#include "AliJetBkg.h"
54
55ClassImp(AliJetBkg)
56
57////////////////////////////////////////////////////////////////////////
58
59AliJetBkg::AliJetBkg():
856618e7 60 TObject(),
61 fReader(0),
62 fHeader(0),
63 fInputFJ(0)
50c7d9f7 64{
65 // Default constructor
66
67}
68//______________________________________________________________________
69AliJetBkg::AliJetBkg(const AliJetBkg& input):
856618e7 70 TObject(input),
50c7d9f7 71 fReader(input.fReader),
72 fHeader(input.fHeader),
856618e7 73 fInputFJ(input.fInputFJ)
50c7d9f7 74{
75 // copy constructor
76
77}
78//______________________________________________________________________
79AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){
80 // Assignment operator.
81 this->~AliJetBkg();
82 new(this) AliJetBkg(source);
83 return *this;
84
85}
86//___________________________________________________________________
87Float_t AliJetBkg::BkgFastJet(){
88
967c266c 89 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
90 Bool_t debug = header->GetDebug(); // debug option
91
92 if(debug)cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
50c7d9f7 93 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
94
967c266c 95 if(debug)cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
50c7d9f7 96
97 for(UInt_t i=0;i<inputParticles.size();i++){
98 // cout<<" "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
99
100 }
101
50c7d9f7 102 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
103 Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
967c266c 104 if(debug)cout<<"-------- rho (from all part)="<<rho<<endl;
50c7d9f7 105 return rho;
106
107}
108//___________________________________________________________________
109Float_t AliJetBkg::BkgChargedFastJet(){
110
967c266c 111 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
112 Bool_t debug = header->GetDebug(); // debug option
113
114 if(debug)cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
50c7d9f7 115
116 vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
117
967c266c 118 if(debug)cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
50c7d9f7 119
120 for(UInt_t i=0;i<inputParticlesCharged.size();i++){
121 // cout<<" "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;
122
123 }
967c266c 124
50c7d9f7 125 double rParam = header->GetRparam();
126
127 Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
128
129 cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
130 return rho;
131}
132
133
134
135Float_t AliJetBkg::BkgStat()
136{
137 //background subtraction using statistical method
138
967c266c 139 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
140 Bool_t debug = header->GetDebug(); // debug option
141
142 if(debug)cout<<"==============AliJetBkg::BkgStat()============="<<endl;
50c7d9f7 143 //TO BE IMPLEMENTED
144 // Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
145 Int_t nTracks= 0;
146 TF1 fun("fun",BkgFunction,0,800,1);
147 Double_t enTot=fun.Eval(nTracks);
148 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
149 return enTot/accEMCal;
150
151}
50c7d9f7 152
153////////////////////////////////////////////////////////////////////////
154Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
155{
156
157 // Cone background subtraction method applied on the fastjet: REmove the particles of the
158 // two largest jets with the given R from the estimation of new rho.
50c7d9f7 159
160 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
967c266c 161 Bool_t debug = header->GetDebug(); // debug option
162
163 if(debug)cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
164
50c7d9f7 165 Float_t rc= header->GetRparam();
166
167 //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
168 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
169
170 Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
967c266c 171 if(debug)cout<<"nJets: "<<nJ<<endl;
50c7d9f7 172
173
174 //begin unit array
175 TClonesArray* fUnit = fReader->GetUnitArray();
176 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
177
178 Int_t nIn = fUnit->GetEntries();
179 if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
180
181 // Information extracted from fUnitArray
182 // load input vectors and calculate total energy in array
183 Float_t pt,eta,phi;
184 Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0;
185 Float_t rhoback=0.0;
186
187 Float_t ptallback=0.0; //particles without the jet
188 Float_t restarea=accEMCal; //initial area set
189 Bool_t acc=0;
190 Bool_t acc1=0;
191 Float_t rCone=0.4;
192
193 if(nJ==1) {
194 AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
195 jeteta=jettmp->Eta();
196 jetphi=jettmp->Phi();
197 acc=EmcalAcceptance(jeteta,jetphi,rCone);
198 if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
967c266c 199 if(debug)cout<<" acc "<<acc<<endl;
50c7d9f7 200 }
201
202
203 if(nJ>=2) {
204 AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
205 AliAODJet *jettmp1 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
206 jeteta=jettmp->Eta();
207 jetphi=jettmp->Phi();
208 jeteta1=jettmp1->Eta();
209 jetphi1=jettmp1->Phi();
210 acc=EmcalAcceptance(jeteta,jetphi,rCone);
211 acc1=EmcalAcceptance(jeteta1,jetphi1,rCone);
212 if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc;
213 if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc;
214 if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
215
967c266c 216 if(debug)cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl;
50c7d9f7 217
218 }
219
220 // cout<<" nIn = "<<nIn<<endl;
221 Float_t sumpt=0;
222 for(Int_t i=0; i<nIn; i++)
223 { //Unit Array Loop
224 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
225 if(uArray->GetUnitEnergy()>0.){
226
227 pt = uArray->GetUnitEnergy();
228 eta = uArray->GetUnitEta();
229 phi = uArray->GetUnitPhi();
230
231 //cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl;
232
233 Float_t deta=0.0, dphi=0.0, dr=100.0;
234 Float_t deta1=0.0, dphi1=0.0, dr1=100.0;
235
236 //cout<<i<<" pt="<<pt<<" eta="<<eta<<" phi="<<phi<<endl;
237 if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
238 sumpt+=pt;
239 //if(i<30)cout<<i<<" pt unit = "<<pt<<endl;
240
241 if(nJ==1 && acc==1) {
242 deta = eta - jeteta;
243 dphi = phi - jetphi;
244 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
245 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
246 dr = TMath::Sqrt(deta * deta + dphi * dphi);
247 if(dr<=rc)sumpt-=pt;
248 }
249
250 if(nJ>=2) {
251 if(acc==1){
252 deta = eta - jeteta;
253 dphi = phi - jetphi;
254 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
255 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
256 dr = TMath::Sqrt(deta * deta + dphi * dphi);
257 if(dr<=rc)sumpt-=pt;
258 }
259 if(acc1==1){
260 deta1 = eta - jeteta1;
261 dphi1 = phi - jetphi1;
262 if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi();
263 if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1;
264 dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1);
265 if(dr1<=rc)sumpt-=pt;
266 }
267
268 }
269
270 if(dr >= rc && dr1 >=rc) {
271 // particles outside both cones
272
273 //cout<<" out of the cone "<<dr<<" "<<deta<<" deltaeta "<<dphi<<" dphi "<<i<<" particle "<<endl;
274 //cout<<" out of the cone "<<dr1<<" "<<deta1<<" deltaeta1 "<<dphi1<<" dphi1 "<<i<<" particle "<<endl;
275 ptallback+=pt;
276 }
277 }
278 //cout<<" ipart "<<ipart<<" rhointegral "<<rhoback <<endl;
279 }
280 } // End loop on UnitArray
281
967c266c 282 if(debug)cout<<"total area left "<<restarea<<endl;
283 if(debug)cout<<"sumpt="<<sumpt<<endl;
50c7d9f7 284 // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
285 //else rhoback=ptallback;
286
287 rhoback= ptallback/restarea;
967c266c 288 if(debug)cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl;
50c7d9f7 289
290 return rhoback;
291
292}
293
294
295Double_t AliJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method){
296 //calculate rho using the fastjet method
297
298 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
299 Bool_t debug = header->GetDebug(); // debug option
300
301 fastjet::Strategy strategy = header->GetStrategy();
302 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
303 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
304 fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
305
306 // create an object that specifies how we to define the area
307 fastjet::AreaDefinition areaDef;
308 double ghostEtamax = header->GetGhostEtaMax();
309 double ghostArea = header->GetGhostArea();
310 int activeAreaRepeats = header->GetActiveAreaRepeats();
311
312 // now create the object that holds info about ghosts
313
314 if (method.Contains("Charg"))ghostEtamax=0.9;
315
316 fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
317 // and from that get an area definition
318 fastjet::AreaType areaType = header->GetAreaType();
319 areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
967c266c 320 if(debug)cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
50c7d9f7 321 //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
322 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
323 TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
324 comment+= "Jet definition: ";
325 comment+= TString(jetDef.description());
326 // comment+= ". Area definition: ";
327 //comment+= TString(areaDef.description());
328 comment+= ". Strategy adopted by FastJet: ";
329 comment+= TString(clust_seq.strategy_string());
330 header->SetComment(comment);
331 if(debug){
332 cout << "--------------------------------------------------------" << endl;
333 cout << comment << endl;
334 cout << "--------------------------------------------------------" << endl;
335 }
336
337 double ptmin = header->GetPtMin();
338 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
339 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
50c7d9f7 340
db2272d5 341 if (debug) {
342 cout<<"# of BKG jets = "<<jets.size()<<endl;
343 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
344 printf("BKG Jet found %5d %9.5f %8.5f %10.3f %4.4f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(),clust_seq.area(jets[j]));
345 }
346 }
347
50c7d9f7 348 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
349
350 if (method.Contains("All")){
351 phiMin = 80.*TMath::Pi()/180+rParamBkg;
352 phiMax = 190.*TMath::Pi()/180-rParamBkg;
50c7d9f7 353 }
354 if (method.Contains("Charg")){
355 phiMin = 0;
356 phiMax = 2*TMath::Pi();
357 }
358 rapMax = ghostEtamax - rParamBkg;
359 rapMin = - ghostEtamax + rParamBkg;
360
361
362 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
363
364
365 Double_t rho=clust_seq.median_pt_per_unit_area(range);
366 // double median, sigma, meanArea;
367 //clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true);
368 //fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea);
369
370 // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);
371
372
967c266c 373 if(debug)cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl;
50c7d9f7 374 return rho;
375}
376
377Float_t AliJetBkg::EtaToTheta(Float_t arg)
378{
379 // return (180./TMath::Pi())*2.*atan(exp(-arg));
380 return 2.*atan(exp(-arg));
381}
382
383
384Double_t AliJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/){
385 //to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity)
386 return 1;
387}
388
389
390Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
391
392 Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
393 Float_t deltaphi=110./180.*TMath::Pi();
394 Float_t phicut=deltaphi/2.-radius;
395 Float_t etacut=0.7-radius;
396 //cout<<" eta "<<eta<<" phi "<<phi<<endl;
397 //cout<<etacut<<" "<<phicut<<" "<<meanPhi<<" "<<endl;
398 if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1;
399 else return 0;
400}