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