]>
Commit | Line | Data |
---|---|---|
d89b8229 | 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 | // Class to calculate the background per unit area | |
20 | // manages the search for jets | |
21 | // Authors: Elena Bruna elena.bruna@yale.edu | |
22 | // Sevil Salur ssalur@lbl.gov | |
23 | // | |
24 | // 2011 : | |
25 | // renamed from AliJetBkg to AliFastJetBkg as this class uses only FASTJET based algos. | |
26 | //--------------------------------------------------------------------- | |
27 | ||
28 | #include <Riostream.h> | |
29 | #include <TClonesArray.h> | |
30 | #include <TF1.h> | |
31 | #include <TString.h> | |
32 | ||
33 | #include "AliJetHeader.h" | |
34 | #include "AliFastJetHeaderV1.h" | |
35 | #include "AliAODJet.h" | |
36 | #include "AliFastJetInput.h" | |
37 | #include "AliFastJetBkg.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(AliFastJetBkg) | |
49 | ||
50 | //////////////////////////////////////////////////////////////////////// | |
51 | ||
52 | AliFastJetBkg::AliFastJetBkg(): | |
53 | TObject(), | |
54 | fHeader(0), | |
55 | fInputFJ(0) | |
56 | { | |
57 | // Default constructor | |
58 | } | |
59 | ||
60 | //______________________________________________________________________ | |
61 | AliFastJetBkg::AliFastJetBkg(const AliFastJetBkg& input): | |
62 | TObject(input), | |
63 | fHeader(input.fHeader), | |
64 | fInputFJ(input.fInputFJ) | |
65 | { | |
66 | // copy constructor | |
67 | } | |
68 | ||
69 | //______________________________________________________________________ | |
70 | AliFastJetBkg& AliFastJetBkg::operator=(const AliFastJetBkg& source) | |
71 | { | |
72 | // Assignment operator. | |
73 | if(this!=&source){ | |
74 | TObject::operator=(source); | |
75 | fHeader = source.fHeader; | |
76 | fInputFJ = source.fInputFJ; | |
77 | } | |
78 | ||
79 | return *this; | |
80 | ||
81 | } | |
82 | ||
83 | //___________________________________________________________________ | |
84 | void AliFastJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma, | |
85 | Double_t& meanarea) | |
86 | { | |
87 | // Bkg estimation | |
88 | ||
89 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 90 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 91 | if(debug>0) cout<<"=============== AliFastJetBkg::BkgFastJetb() =========== "<<endl; |
92 | vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles(); | |
93 | ||
94 | double rParamBkg = header->GetRparamBkg(); //Radius for background calculation | |
95 | ||
96 | Double_t medianb,sigmab,meanareab; | |
97 | CalcRhob(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All"); | |
98 | rho=medianb; | |
99 | sigma=sigmab; | |
100 | meanarea=meanareab; | |
101 | ||
102 | } | |
103 | ||
104 | //_________________________________________________________________ | |
105 | void AliFastJetBkg::BkgFastJetWoHardest(Double_t& rho,Double_t& sigma, | |
106 | Double_t& meanarea) | |
107 | { | |
108 | ||
109 | // Bkg estimation without hardest jet | |
110 | ||
111 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 112 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 113 | if(debug) cout<<"=============== AliFastJetBkg::BkgWoHardest() =========== "<<endl; |
114 | vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles(); | |
115 | ||
116 | double rParamBkg = header->GetRparamBkg(); //Radius for background calculation | |
117 | Double_t medianb,sigmab,meanareab; | |
118 | CalcRhoWoHardest(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All"); | |
119 | rho=medianb; | |
120 | sigma=sigmab; | |
121 | meanarea=meanareab; | |
122 | ||
123 | } | |
124 | ||
125 | //____________________________________________________________________ | |
126 | void AliFastJetBkg::CalcRhob(Double_t& median,Double_t& | |
127 | sigma,Double_t& | |
128 | meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t | |
129 | rParamBkg,TString method) | |
130 | { | |
131 | // calculate rho using the fastjet method | |
132 | ||
133 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 134 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 135 | |
136 | fastjet::Strategy strategy = header->GetStrategy(); | |
137 | fastjet::RecombinationScheme recombScheme = header->GetRecombScheme(); | |
138 | fastjet::JetAlgorithm algorithm = header->GetBGAlgorithm(); | |
139 | fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy); | |
140 | ||
141 | // create an object that specifies how we to define the area | |
142 | fastjet::AreaDefinition areaDef; | |
143 | double ghostEtamax = header->GetGhostEtaMax(); | |
144 | double ghostArea = header->GetGhostArea(); | |
145 | int activeAreaRepeats = header->GetActiveAreaRepeats(); | |
146 | ||
147 | // now create the object that holds info about ghosts | |
148 | ||
149 | fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea); | |
150 | // and from that get an area definition | |
151 | fastjet::AreaType areaType = header->GetAreaType(); | |
152 | areaDef = fastjet::AreaDefinition(areaType,ghost_spec); | |
153 | ||
154 | //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef); | |
155 | fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef); | |
156 | TString comment = "Running FastJet algorithm for BKG calculation with the following setup. "; | |
157 | comment+= "Jet definition: "; | |
158 | comment+= TString(jetDef.description()); | |
159 | // comment+= ". Area definition: "; | |
160 | // comment+= TString(areaDef.description()); | |
161 | comment+= ". Strategy adopted by FastJet: "; | |
162 | comment+= TString(clust_seq.strategy_string()); | |
163 | comment+= Form("Method: %s",method.Data()); | |
164 | header->SetComment(comment); | |
165 | if(debug>0){ | |
166 | cout << "--------------------------------------------------------" << endl; | |
167 | cout << comment << endl; | |
168 | cout << "--------------------------------------------------------" << endl; | |
169 | } | |
170 | ||
171 | vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.); | |
172 | vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); | |
173 | ||
174 | double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0; | |
175 | ||
176 | phiMin = 0; | |
177 | phiMax = 2*TMath::Pi(); | |
178 | ||
179 | rapMax = ghostEtamax - rParamBkg; | |
180 | rapMin = - ghostEtamax + rParamBkg; | |
181 | ||
182 | fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); | |
183 | ||
184 | double medianb, sigmab, meanareab; | |
185 | clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, medianb, sigmab, meanareab, false); | |
186 | median=medianb; | |
187 | sigma=sigmab; | |
188 | meanarea=meanareab; | |
189 | ||
190 | } | |
191 | ||
192 | //____________________________________________________________________ | |
193 | void AliFastJetBkg::CalcRhoWoHardest(Double_t& median,Double_t& | |
194 | sigma,Double_t& meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t | |
195 | rParamBkg,TString method) | |
196 | { | |
197 | // calculate rho (without the hardest jet) using the fastjet method | |
198 | ||
199 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 200 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 201 | |
202 | fastjet::Strategy strategy = header->GetStrategy(); | |
203 | fastjet::RecombinationScheme recombScheme = header->GetRecombScheme(); | |
204 | fastjet::JetAlgorithm algorithm = header->GetBGAlgorithm(); | |
205 | fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy); | |
206 | ||
207 | // create an object that specifies how we to define the area | |
208 | fastjet::AreaDefinition areaDef; | |
209 | double ghostEtamax = header->GetGhostEtaMax(); | |
210 | double ghostArea = header->GetGhostArea(); | |
211 | int activeAreaRepeats = header->GetActiveAreaRepeats(); | |
212 | ||
213 | // now create the object that holds info about ghosts | |
214 | ||
215 | fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea); | |
216 | // and from that get an area definition | |
217 | fastjet::AreaType areaType = header->GetAreaType(); | |
218 | areaDef = fastjet::AreaDefinition(areaType,ghost_spec); | |
219 | //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef); | |
220 | fastjet::ClusterSequenceArea clust_seq(inputParticles,jetDef,areaDef); | |
221 | TString comment = "Running FastJet algorithm for BKG calculation with the following setup. "; | |
222 | comment+= "Jet definition: "; | |
223 | comment+= TString(jetDef.description()); | |
224 | // comment+= ". Area definition: "; | |
225 | // comment+= TString(areaDef.description()); | |
226 | comment+= ". Strategy adopted by FastJet: "; | |
227 | comment+= TString(clust_seq.strategy_string()); | |
228 | comment+= Form("Method: %s",method.Data()); | |
229 | header->SetComment(comment); | |
230 | if(debug>0){ | |
231 | cout << "--------------------------------------------------------" << endl; | |
232 | cout << comment << endl; | |
233 | cout << "--------------------------------------------------------" << endl; | |
234 | } | |
235 | ||
236 | vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.); | |
237 | vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); | |
238 | vector<fastjet::PseudoJet> jets2=sorted_by_pt(inclusiveJets); | |
239 | if(jets2.size()>=2) jets2.erase(jets2.begin(),jets2.begin()+1); | |
240 | ||
241 | double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0; | |
242 | ||
243 | phiMin = 0; | |
244 | phiMax = 2*TMath::Pi(); | |
245 | ||
246 | rapMax = ghostEtamax - rParamBkg; | |
247 | rapMin = - ghostEtamax + rParamBkg; | |
248 | ||
249 | fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); | |
250 | ||
251 | double medianb, sigmab, meanareab; | |
252 | clust_seq.get_median_rho_and_sigma(jets2, range, false, medianb, sigmab, | |
253 | meanareab, false); | |
254 | median=medianb; | |
255 | sigma=sigmab; | |
256 | meanarea=meanareab; | |
257 | ||
258 | } | |
259 | ||
260 | //___________________________________________________________________ | |
261 | Float_t AliFastJetBkg::BkgFastJet() | |
262 | { | |
263 | // Return background | |
264 | ||
265 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 266 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 267 | |
268 | if(debug>0) cout<<"=============== AliFastJetBkg::BkgFastJet() =========== "<<endl; | |
269 | vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles(); | |
270 | ||
271 | if(debug>0) cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl; | |
272 | ||
273 | double rParamBkg = header->GetRparamBkg(); //Radius for background calculation | |
274 | Double_t rho=CalcRho(inputParticles,rParamBkg,"All"); | |
275 | if(debug) cout<<"-------- rho (from all part)="<<rho<<endl; | |
276 | return rho; | |
277 | ||
278 | } | |
279 | ||
280 | //___________________________________________________________________ | |
281 | Float_t AliFastJetBkg::BkgChargedFastJet() | |
282 | { | |
283 | // Background for charged jets | |
284 | ||
285 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 286 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 287 | |
288 | if(debug>0) cout<<"=============== AliFastJetBkg::BkgChargedFastJet() =========== "<<endl; | |
289 | ||
290 | vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh(); | |
291 | ||
292 | if(debug>0) cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl; | |
293 | ||
294 | double rParam = header->GetRparam(); | |
295 | ||
296 | Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg"); | |
297 | ||
298 | if(debug>0) cout<<"-------- rho (from CHARGED part)="<<rho<<endl; | |
299 | return rho; | |
300 | ||
301 | } | |
302 | ||
303 | //___________________________________________________________________ | |
304 | Float_t AliFastJetBkg::BkgStat() | |
305 | { | |
306 | // background subtraction using statistical method | |
307 | ||
308 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 309 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 310 | |
311 | if(debug>0) cout<<"==============AliFastJetBkg::BkgStat()============="<<endl; | |
312 | //TO BE IMPLEMENTED | |
313 | Int_t nTracks= 0; | |
314 | TF1 fun("fun",BkgFunction,0,800,1); | |
315 | Double_t enTot=fun.Eval(nTracks); | |
316 | Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal | |
317 | return enTot/accEMCal; | |
318 | ||
319 | } | |
320 | ||
321 | //___________________________________________________________________ | |
322 | Float_t AliFastJetBkg::BkgFastJetCone(TClonesArray* fAODJets) | |
323 | { | |
324 | // Cone background subtraction method applied on the fastjet: REmove the particles of the | |
325 | // two largest jets with the given R from the estimation of new rho. | |
326 | ||
327 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 328 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 329 | |
330 | if(debug>0) cout<<"==============AliFastJetBkg::SubtractFastJetBackgCone()============="<<endl; | |
331 | ||
332 | Float_t rc= header->GetRparam(); | |
333 | ||
334 | //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h) | |
335 | Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal | |
336 | ||
337 | Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets... | |
338 | if(debug>0) cout<<"nJets: "<<nJ<<endl; | |
339 | ||
340 | // Information extracted from fInputParticle | |
341 | // load input vectors and calculate total energy in array | |
342 | Float_t pt,eta,phi; | |
343 | Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0; | |
344 | Float_t rhoback=0.0; | |
345 | ||
346 | Float_t ptallback=0.0; //particles without the jet | |
347 | Float_t restarea=accEMCal; //initial area set | |
348 | Bool_t acc=0; | |
349 | Bool_t acc1=0; | |
350 | Float_t rCone=0.4; | |
351 | ||
352 | if(nJ==1) { | |
353 | AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0)); | |
354 | jeteta=jettmp->Eta(); | |
355 | jetphi=jettmp->Phi(); | |
356 | acc=EmcalAcceptance(jeteta,jetphi,rCone); | |
357 | if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc; | |
358 | if(debug) cout<<" acc "<<acc<<endl; | |
359 | } | |
360 | ||
361 | if(nJ>=2) { | |
362 | AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0)); | |
363 | AliAODJet *jettmp1 = (AliAODJet*)(fAODJets->At(1)); | |
364 | jeteta=jettmp->Eta(); | |
365 | jetphi=jettmp->Phi(); | |
366 | jeteta1=jettmp1->Eta(); | |
367 | jetphi1=jettmp1->Phi(); | |
368 | acc=EmcalAcceptance(jeteta,jetphi,rCone); | |
369 | acc1=EmcalAcceptance(jeteta1,jetphi1,rCone); | |
370 | if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc; | |
371 | if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc; | |
372 | if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc; | |
373 | ||
374 | if(debug) cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl; | |
375 | ||
376 | } | |
377 | ||
378 | // cout<<" nIn = "<<nIn<<endl; | |
379 | Float_t sumpt=0; | |
380 | vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles(); | |
381 | for(UInt_t i=0; i<inputParticles.size(); i++) | |
382 | { // Loop over input list of particles | |
383 | pt = inputParticles[i].perp(); | |
384 | eta = inputParticles[i].eta(); | |
385 | phi = inputParticles[i].phi(); | |
386 | ||
387 | // To be updated | |
388 | //cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl; | |
389 | ||
390 | Float_t deta=0.0, dphi=0.0, dr=100.0; | |
391 | Float_t deta1=0.0, dphi1=0.0, dr1=100.0; | |
392 | ||
393 | //cout<<i<<" pt="<<pt<<" eta="<<eta<<" phi="<<phi<<endl; | |
394 | if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){ | |
395 | sumpt+=pt; | |
396 | //if(i<30)cout<<i<<" pt = "<<pt<<endl; | |
397 | ||
398 | if(nJ==1 && acc==1) { | |
399 | deta = eta - jeteta; | |
400 | dphi = phi - jetphi; | |
401 | if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); | |
402 | if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; | |
403 | dr = TMath::Sqrt(deta * deta + dphi * dphi); | |
404 | if(dr<=rc)sumpt-=pt; | |
405 | } | |
406 | ||
407 | if(nJ>=2) { | |
408 | if(acc==1){ | |
409 | deta = eta - jeteta; | |
410 | dphi = phi - jetphi; | |
411 | if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); | |
412 | if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; | |
413 | dr = TMath::Sqrt(deta * deta + dphi * dphi); | |
414 | if(dr<=rc)sumpt-=pt; | |
415 | } | |
416 | if(acc1==1){ | |
417 | deta1 = eta - jeteta1; | |
418 | dphi1 = phi - jetphi1; | |
419 | if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi(); | |
420 | if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1; | |
421 | dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1); | |
422 | if(dr1<=rc)sumpt-=pt; | |
423 | } | |
424 | } | |
425 | ||
426 | if(dr >= rc && dr1 >=rc) { | |
427 | // particles outside both cones | |
428 | if(debug>1) cout<<" out of the cone "<<dr<<" "<<deta<<" deltaeta "<<dphi<<" dphi "<<i<<" particle "<<endl; | |
429 | if(debug>1) cout<<" out of the cone "<<dr1<<" "<<deta1<<" deltaeta1 "<<dphi1<<" dphi1 "<<i<<" particle "<<endl; | |
430 | ptallback+=pt; | |
431 | } | |
432 | } | |
433 | } // End loop on input list of particles | |
434 | ||
435 | if(debug>0) cout<<"total area left "<<restarea<<endl; | |
436 | if(debug>0) cout<<"sumpt="<<sumpt<<endl; | |
437 | // if(acc==1 || acc1==1) rhoback= ptallback/restarea; | |
438 | // else rhoback=ptallback; | |
439 | ||
440 | rhoback= ptallback/restarea; | |
441 | if(debug)cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl; | |
442 | ||
443 | return rhoback; | |
444 | ||
445 | } | |
446 | ||
447 | //___________________________________________________________________ | |
448 | Double_t AliFastJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method) | |
449 | { | |
450 | // calculate rho using the fastjet method | |
451 | ||
452 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
37f45f89 | 453 | Int_t debug = header->GetDebug(); // debug option |
d89b8229 | 454 | |
455 | fastjet::Strategy strategy = header->GetStrategy(); | |
456 | fastjet::RecombinationScheme recombScheme = header->GetRecombScheme(); | |
457 | fastjet::JetAlgorithm algorithm = header->GetBGAlgorithm(); | |
458 | fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy); | |
459 | ||
460 | // create an object that specifies how we to define the area | |
461 | fastjet::AreaDefinition areaDef; | |
462 | double ghostEtamax = header->GetGhostEtaMax(); | |
463 | double ghostArea = header->GetGhostArea(); | |
464 | int activeAreaRepeats = header->GetActiveAreaRepeats(); | |
465 | ||
466 | // now create the object that holds info about ghosts | |
467 | ||
468 | if (method.Contains("Charg"))ghostEtamax=0.9; | |
469 | ||
470 | fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea); | |
471 | // and from that get an area definition | |
472 | fastjet::AreaType areaType = header->GetAreaType(); | |
473 | areaDef = fastjet::AreaDefinition(areaType,ghost_spec); | |
474 | if(debug>0) cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<< endl; | |
475 | //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef); | |
476 | fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef); | |
477 | TString comment = "Running FastJet algorithm for BKG calculation with the following setup. "; | |
478 | comment+= "Jet definition: "; | |
479 | comment+= TString(jetDef.description()); | |
480 | // comment+= ". Area definition: "; | |
481 | // comment+= TString(areaDef.description()); | |
482 | comment+= ". Strategy adopted by FastJet: "; | |
483 | comment+= TString(clust_seq.strategy_string()); | |
484 | header->SetComment(comment); | |
485 | if(debug>0){ | |
486 | cout << "--------------------------------------------------------" << endl; | |
487 | cout << comment << endl; | |
488 | cout << "--------------------------------------------------------" << endl; | |
489 | } | |
490 | ||
491 | double ptmin = header->GetPtMin(); | |
492 | vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin); | |
493 | vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); | |
494 | ||
495 | if (debug>0) { | |
496 | cout<<"# of BKG jets = "<<jets.size()<<endl; | |
497 | for (size_t j = 0; j < jets.size(); j++) { // loop for jets | |
498 | 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])); | |
499 | } | |
500 | } | |
501 | ||
502 | double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0; | |
503 | ||
504 | if (method.Contains("All")){ | |
505 | phiMin = 80.*TMath::Pi()/180+rParamBkg; | |
506 | phiMax = 190.*TMath::Pi()/180-rParamBkg; | |
507 | } | |
508 | if (method.Contains("Charg")){ | |
509 | phiMin = 0; | |
510 | phiMax = 2*TMath::Pi(); | |
511 | } | |
512 | rapMax = ghostEtamax - rParamBkg; | |
513 | rapMin = - ghostEtamax + rParamBkg; | |
514 | ||
515 | fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); | |
516 | ||
517 | Double_t rho=clust_seq.median_pt_per_unit_area(range); | |
518 | // double median, sigma, meanArea; | |
519 | // clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true); | |
520 | // fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea); | |
521 | ||
522 | // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec); | |
523 | ||
524 | if(debug>0) cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl; | |
525 | ||
526 | return rho; | |
527 | ||
528 | } | |
529 | ||
530 | //___________________________________________________________________ | |
531 | Double_t AliFastJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/) | |
532 | { | |
533 | // to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity) | |
534 | return 1; | |
535 | ||
536 | } | |
537 | ||
538 | //___________________________________________________________________ | |
7aa72c24 | 539 | Bool_t AliFastJetBkg::EmcalAcceptance(Float_t eta, Float_t phi, Float_t radius) const |
d89b8229 | 540 | { |
541 | // Apply emcal acceptance cuts | |
542 | // To be updated | |
543 | ||
544 | Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2; | |
545 | Float_t deltaphi=110./180.*TMath::Pi(); | |
546 | Float_t phicut=deltaphi/2.-radius; | |
547 | Float_t etacut=0.7-radius; | |
548 | //cout<<" eta "<<eta<<" phi "<<phi<<endl; | |
549 | //cout<<etacut<<" "<<phicut<<" "<<meanPhi<<" "<<endl; | |
550 | if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1; | |
551 | else return 0; | |
552 | ||
553 | } |