1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 #include <Riostream.h>
21 #include <TClonesArray.h>
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"
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"
44 #ifdef ENABLE_PLUGIN_SISCONE
45 #include "fastjet/SISConePlugin.hh"
48 #include<sstream> // needed for internal io
53 #include "AliJetBkg.h"
57 ////////////////////////////////////////////////////////////////////////
59 AliJetBkg::AliJetBkg():
65 // Default constructor
68 //______________________________________________________________________
69 AliJetBkg::AliJetBkg(const AliJetBkg& input):
71 fReader(input.fReader),
72 fHeader(input.fHeader),
73 fInputFJ(input.fInputFJ)
78 //______________________________________________________________________
79 AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){
80 // Assignment operator.
82 new(this) AliJetBkg(source);
88 //___________________________________________________________________
89 void AliJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma,
92 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
93 Bool_t debug = header->GetDebug(); // debug option
94 if(debug)cout<<"=============== AliJetBkg::BkgFastJetb() =========== "<<endl;
95 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
97 //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
102 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
104 Double_t medianb,sigmab,meanareab;
105 CalcRhob(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
111 //_________________________________________________________________
113 void AliJetBkg::BkgFastJetWoHardest(Double_t& rho,Double_t& sigma,
115 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
116 Bool_t debug = header->GetDebug(); // debug option
117 if(debug)cout<<"=============== AliJetBkg::BkgWoHardest() =========== "<<endl;
118 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
120 //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
121 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
122 Double_t medianb,sigmab,meanareab;
123 CalcRhoWoHardest(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
130 //____________________________________________________________________
132 void AliJetBkg::CalcRhob(Double_t& median,Double_t&
134 meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
135 rParamBkg,TString method){
136 //calculate rho using the fastjet method
138 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
139 Bool_t debug = header->GetDebug(); // debug option
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);
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();
152 // now create the object that holds info about ghosts
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);
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());
170 header->SetComment(comment);
172 cout << "--------------------------------------------------------" << endl;
173 cout << comment << endl;
174 cout << "--------------------------------------------------------" << endl;
177 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
178 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
181 //cout<<"# of BKG jets = "<<jets.size()<<endl;
182 // for (size_t j = 0; j < jets.size(); j++) { // loop for jets
183 //printf("BKG Jet found %5d %9.5f %8.5f %10.3f %4.4f
184 //\n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(),
185 //clust_seq.area(jets[j]));
189 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
193 phiMax = 2*TMath::Pi();
195 rapMax = ghostEtamax - rParamBkg;
196 rapMin = - ghostEtamax + rParamBkg;
199 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
202 double medianb, sigmab, meanareab;
203 clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, medianb, sigmab, meanareab, false);
211 //____________________________________________________________________
214 void AliJetBkg::CalcRhoWoHardest(Double_t& median,Double_t&
216 meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
217 rParamBkg,TString method){
218 //calculate rho using the fastjet method
220 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
221 Bool_t debug = header->GetDebug(); // debug option
223 fastjet::Strategy strategy = header->GetStrategy();
224 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
225 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
226 fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
228 // create an object that specifies how we to define the area
229 fastjet::AreaDefinition areaDef;
230 double ghostEtamax = header->GetGhostEtaMax();
231 double ghostArea = header->GetGhostArea();
232 int activeAreaRepeats = header->GetActiveAreaRepeats();
234 // now create the object that holds info about ghosts
238 fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
239 // and from that get an area definition
240 fastjet::AreaType areaType = header->GetAreaType();
241 areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
242 //cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
243 //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
244 fastjet::ClusterSequenceArea clust_seq(inputParticles,jetDef,areaDef);
245 TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
246 comment+= "Jet definition: ";
247 comment+= TString(jetDef.description());
248 // comment+= ". Area definition: ";
249 //comment+= TString(areaDef.description());
250 comment+= ". Strategy adopted by FastJet: ";
251 comment+= TString(clust_seq.strategy_string());
252 header->SetComment(comment);
254 cout << "--------------------------------------------------------" << endl;
255 cout << comment << endl;
256 cout << "--------------------------------------------------------" << endl;
259 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
260 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
261 vector<fastjet::PseudoJet> jets2=sorted_by_pt(inclusiveJets);
262 if(jets2.size()>=2) jets2.erase(jets2.begin(),jets2.begin()+1);
266 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
270 phiMax = 2*TMath::Pi();
272 rapMax = ghostEtamax - rParamBkg;
273 rapMin = - ghostEtamax + rParamBkg;
276 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
279 double medianb, sigmab, meanareab;
280 clust_seq.get_median_rho_and_sigma(jets2, range, false, medianb, sigmab,
290 //___________________________________________________________________
291 Float_t AliJetBkg::BkgFastJet(){
293 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
294 Bool_t debug = header->GetDebug(); // debug option
296 if(debug)cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
297 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
299 if(debug)cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
301 for(UInt_t i=0;i<inputParticles.size();i++){
302 // cout<<" "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
306 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
307 Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
308 if(debug)cout<<"-------- rho (from all part)="<<rho<<endl;
312 //___________________________________________________________________
313 Float_t AliJetBkg::BkgChargedFastJet(){
315 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
316 Bool_t debug = header->GetDebug(); // debug option
318 if(debug)cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
320 vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
322 if(debug)cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
324 for(UInt_t i=0;i<inputParticlesCharged.size();i++){
325 // cout<<" "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;
329 double rParam = header->GetRparam();
331 Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
333 if(debug)cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
339 Float_t AliJetBkg::BkgStat()
341 //background subtraction using statistical method
343 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
344 Bool_t debug = header->GetDebug(); // debug option
346 if(debug)cout<<"==============AliJetBkg::BkgStat()============="<<endl;
348 // Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
350 TF1 fun("fun",BkgFunction,0,800,1);
351 Double_t enTot=fun.Eval(nTracks);
352 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
353 return enTot/accEMCal;
357 ////////////////////////////////////////////////////////////////////////
358 Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
361 // Cone background subtraction method applied on the fastjet: REmove the particles of the
362 // two largest jets with the given R from the estimation of new rho.
364 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
365 Bool_t debug = header->GetDebug(); // debug option
367 if(debug)cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
369 Float_t rc= header->GetRparam();
371 //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
372 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
374 Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
375 if(debug)cout<<"nJets: "<<nJ<<endl;
379 TClonesArray* fUnit = fReader->GetUnitArray();
380 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
382 Int_t nIn = fUnit->GetEntries();
383 if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
385 // Information extracted from fUnitArray
386 // load input vectors and calculate total energy in array
388 Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0;
391 Float_t ptallback=0.0; //particles without the jet
392 Float_t restarea=accEMCal; //initial area set
398 AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
399 jeteta=jettmp->Eta();
400 jetphi=jettmp->Phi();
401 acc=EmcalAcceptance(jeteta,jetphi,rCone);
402 if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
403 if(debug)cout<<" acc "<<acc<<endl;
408 AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
409 AliAODJet *jettmp1 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
410 jeteta=jettmp->Eta();
411 jetphi=jettmp->Phi();
412 jeteta1=jettmp1->Eta();
413 jetphi1=jettmp1->Phi();
414 acc=EmcalAcceptance(jeteta,jetphi,rCone);
415 acc1=EmcalAcceptance(jeteta1,jetphi1,rCone);
416 if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc;
417 if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc;
418 if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
420 if(debug)cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl;
424 // cout<<" nIn = "<<nIn<<endl;
426 for(Int_t i=0; i<nIn; i++)
428 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
429 if(uArray->GetUnitEnergy()>0.){
431 pt = uArray->GetUnitEnergy();
432 eta = uArray->GetUnitEta();
433 phi = uArray->GetUnitPhi();
435 //cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl;
437 Float_t deta=0.0, dphi=0.0, dr=100.0;
438 Float_t deta1=0.0, dphi1=0.0, dr1=100.0;
440 //cout<<i<<" pt="<<pt<<" eta="<<eta<<" phi="<<phi<<endl;
441 if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
443 //if(i<30)cout<<i<<" pt unit = "<<pt<<endl;
445 if(nJ==1 && acc==1) {
448 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
449 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
450 dr = TMath::Sqrt(deta * deta + dphi * dphi);
458 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
459 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
460 dr = TMath::Sqrt(deta * deta + dphi * dphi);
464 deta1 = eta - jeteta1;
465 dphi1 = phi - jetphi1;
466 if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi();
467 if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1;
468 dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1);
469 if(dr1<=rc)sumpt-=pt;
474 if(dr >= rc && dr1 >=rc) {
475 // particles outside both cones
477 //cout<<" out of the cone "<<dr<<" "<<deta<<" deltaeta "<<dphi<<" dphi "<<i<<" particle "<<endl;
478 //cout<<" out of the cone "<<dr1<<" "<<deta1<<" deltaeta1 "<<dphi1<<" dphi1 "<<i<<" particle "<<endl;
482 //cout<<" ipart "<<ipart<<" rhointegral "<<rhoback <<endl;
484 } // End loop on UnitArray
486 if(debug)cout<<"total area left "<<restarea<<endl;
487 if(debug)cout<<"sumpt="<<sumpt<<endl;
488 // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
489 //else rhoback=ptallback;
491 rhoback= ptallback/restarea;
492 if(debug)cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl;
499 Double_t AliJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method){
500 //calculate rho using the fastjet method
502 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
503 Bool_t debug = header->GetDebug(); // debug option
505 fastjet::Strategy strategy = header->GetStrategy();
506 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
507 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
508 fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
510 // create an object that specifies how we to define the area
511 fastjet::AreaDefinition areaDef;
512 double ghostEtamax = header->GetGhostEtaMax();
513 double ghostArea = header->GetGhostArea();
514 int activeAreaRepeats = header->GetActiveAreaRepeats();
516 // now create the object that holds info about ghosts
518 if (method.Contains("Charg"))ghostEtamax=0.9;
520 fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
521 // and from that get an area definition
522 fastjet::AreaType areaType = header->GetAreaType();
523 areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
524 if(debug)cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
525 //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
526 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
527 TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
528 comment+= "Jet definition: ";
529 comment+= TString(jetDef.description());
530 // comment+= ". Area definition: ";
531 //comment+= TString(areaDef.description());
532 comment+= ". Strategy adopted by FastJet: ";
533 comment+= TString(clust_seq.strategy_string());
534 header->SetComment(comment);
536 cout << "--------------------------------------------------------" << endl;
537 cout << comment << endl;
538 cout << "--------------------------------------------------------" << endl;
541 double ptmin = header->GetPtMin();
542 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
543 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
546 cout<<"# of BKG jets = "<<jets.size()<<endl;
547 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
548 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]));
552 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
554 if (method.Contains("All")){
555 phiMin = 80.*TMath::Pi()/180+rParamBkg;
556 phiMax = 190.*TMath::Pi()/180-rParamBkg;
558 if (method.Contains("Charg")){
560 phiMax = 2*TMath::Pi();
562 rapMax = ghostEtamax - rParamBkg;
563 rapMin = - ghostEtamax + rParamBkg;
566 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
569 Double_t rho=clust_seq.median_pt_per_unit_area(range);
570 // double median, sigma, meanArea;
571 //clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true);
572 //fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea);
574 // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);
577 if(debug)cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl;
581 Float_t AliJetBkg::EtaToTheta(Float_t arg)
583 // return (180./TMath::Pi())*2.*atan(exp(-arg));
584 return 2.*atan(exp(-arg));
588 Double_t AliJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/){
589 //to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity)
594 Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
596 Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
597 Float_t deltaphi=110./180.*TMath::Pi();
598 Float_t phicut=deltaphi/2.-radius;
599 Float_t etacut=0.7-radius;
600 //cout<<" eta "<<eta<<" phi "<<phi<<endl;
601 //cout<<etacut<<" "<<phicut<<" "<<meanPhi<<" "<<endl;
602 if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1;