update to adapt to QA framework modifications
[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
89 cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
90 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
91
92 cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
93
94 for(UInt_t i=0;i<inputParticles.size();i++){
95 // cout<<" "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
96
97 }
98
99 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
100 double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
101 Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
102 cout<<"-------- rho (from all part)="<<rho<<endl;
103 return rho;
104
105}
106//___________________________________________________________________
107Float_t AliJetBkg::BkgChargedFastJet(){
108
109 cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
110
111 vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
112
113 cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
114
115 for(UInt_t i=0;i<inputParticlesCharged.size();i++){
116 // cout<<" "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;
117
118 }
119 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
120 double rParam = header->GetRparam();
121
122 Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
123
124 cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
125 return rho;
126}
127
128
129
130Float_t AliJetBkg::BkgStat()
131{
132 //background subtraction using statistical method
133
134 cout<<"==============AliJetBkg::BkgStat()============="<<endl;
135 //TO BE IMPLEMENTED
136 // Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
137 Int_t nTracks= 0;
138 TF1 fun("fun",BkgFunction,0,800,1);
139 Double_t enTot=fun.Eval(nTracks);
140 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
141 return enTot/accEMCal;
142
143}
144/////////////////////////////////
145Float_t AliJetBkg::BkgRemoveJetLeading(TClonesArray* fAODJets)
146{
147 // Remove the particles of the
148 // two largest jets using the track references stored in the AODJet from the estimation of new rho.
149
150 cout<<"==============AliJetBkg::BkgRemoveJetLeading()============="<<endl;
151
152
153 // check if we are reading AOD jets
154 TRefArray *refs = 0;
155 Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetESDReader");
156 if (fromAod) { refs = fReader->GetReferences(); }
157
158 //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
159 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
160
161 Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
162 cout<<"nJets: "<<nJ<<endl;
163
164
165 //begin unit array
166 TClonesArray* fUnit = fReader->GetUnitArray();
167 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
168
169 Int_t nIn = fUnit->GetEntries();
170 if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
171
172 Float_t rhoback=0.0;
173 Float_t jetarea1=0.0,jetarea2=0.0;
174
175 Int_t particlejet1=-99;
176 Int_t particlejet2=-99;
177 TRefArray *refarray1 = 0;
178 TRefArray *refarray2 = 0;
179 Int_t nJettracks1 = 0, nJettracks2 = 0;
180 Int_t acc=0,acc1=0;
181 AliAODJet *jet1;
182 AliAODJet *jet2;
183
184 if(nJ==1){
185 jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
186 jetarea1=jet1->EffectiveAreaCharged();
187 Float_t jetPhi=jet1->Phi();
188 Float_t jetEta=jet1->Eta();
189 if(jetPhi>1.396 && jetPhi<3.316 && jetEta>-0.7 && jetEta<0.7)acc=1;
190 refarray1=jet1->GetRefTracks();
191 nJettracks1=refarray1->GetEntries();
192 cout<<"nJ = 1, acc="<<acc<<" jetarea1="<<jetarea1<<endl;
193
194 }
195
196 if(nJ>=2){
197 jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
198 jetarea1=jet1->EffectiveAreaCharged();
199 Float_t jetPhi1=jet1->Phi();
200 Float_t jetEta1=jet1->Eta();
201 if(jetPhi1>1.396 && jetPhi1<3.316 && jetEta1>-0.7 && jetEta1<0.7)acc=1;
202 refarray1=jet1->GetRefTracks();
203 nJettracks1=refarray1->GetEntries();
204 cout<<"npart = "<<nJettracks1<<endl;
205
206 jet2 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
207 jetarea2=jet2->EffectiveAreaCharged();
208 Float_t jetPhi2=jet2->Phi();
209 Float_t jetEta2=jet2->Eta();
210 if(jetPhi2>1.396 && jetPhi2<3.316 && jetEta2>-0.7 && jetEta2<0.7)acc1=1;
211 refarray2=jet2->GetRefTracks();
212 nJettracks2=refarray2->GetEntries();
213 cout<<"nJ = "<<nJ<<", acc="<<acc<<" acc1="<<acc1<<" jetarea1="<<jetarea1<<" jetarea2="<<jetarea2<<endl;
214 }
215
216
217
218 // cout<<" nIn = "<<nIn<<endl;
219 Float_t sumPt=0;
220 Float_t eta,phi,pt;
221 Int_t ipart=0;
222
223 for(Int_t i=0; i<nIn; i++)
224 { //Unit Array Loop
225 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
226
227 if(uArray->GetUnitEnergy()>0.){
228 eta = uArray->GetUnitEta();
229 phi = uArray->GetUnitPhi();
230 pt = uArray->GetUnitEnergy();
231 // cout<<"ipart = "<<ipart<<" eta="<<eta<<" phi="<<phi<<endl;
232 if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
233 //cout<<sumPt<<endl;
234 sumPt+=pt;
235
236
237 if(nJ==1 && acc==1){
238 for(Int_t ii=0; ii<nJettracks1;ii++){
239
240 particlejet1 = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
241
242 if(ipart==particlejet1) {
243 sumPt-=pt;
244 }
245 }
246 }
247
248
249 if(nJ>=2){
250
251 //first jet
252 if(acc==1){
253 for(Int_t ii=0; ii<nJettracks1;ii++){
254 particlejet1 = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
255
256 //cout<<"uArr loop = "<<i<<" ipart in uArr (1/2)="<<ipart<<" part in jet="<<ii<<" partID="<<particlejet1<<" sumPt="<<sumPt<<endl;
257 if(ipart==particlejet1) {
258 sumPt-=pt;
259 }
260 }
261 }
262 if(acc1==1){
263 //second jet
264 for(Int_t ii=0; ii<nJettracks2;ii++){
265 particlejet2 = ((AliJetUnitArray*)refarray2->At(ii))->GetUnitTrackID();
266 //cout<<"uArr loop = "<<i<<" ipart in uArr (2/2)="<<ipart<<" part in jet="<<ii<<" partID="<<particlejet2<<" sumPt="<<sumPt<<endl;
267 if(ipart==particlejet2) {
268 sumPt-=pt;
269 }
270 }
271 }
272
273
274 }
275
276 }//if phi,eta
277 ipart++;
278 }//end if energy
279 }// end unit array loop
280
281
282 Float_t areasum=areasum=accEMCal-acc*jetarea1-acc1*jetarea2;
283 cout<<"pt sum "<<sumPt<<" area "<<areasum<<endl;
284
285 if(nJ>0) rhoback=sumPt/areasum;
286 else rhoback=0.;
287 cout<<" rho from leading jet paricle array removed "<<rhoback<<endl;
288
289 return rhoback;
290
291}
292
293
294
295
296////////////////////////////////////////////////////////////////////////
297Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
298{
299
300 // Cone background subtraction method applied on the fastjet: REmove the particles of the
301 // two largest jets with the given R from the estimation of new rho.
302 cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
303
304 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
305 Float_t rc= header->GetRparam();
306
307 //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
308 Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
309
310 Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
311 cout<<"nJets: "<<nJ<<endl;
312
313
314 //begin unit array
315 TClonesArray* fUnit = fReader->GetUnitArray();
316 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
317
318 Int_t nIn = fUnit->GetEntries();
319 if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
320
321 // Information extracted from fUnitArray
322 // load input vectors and calculate total energy in array
323 Float_t pt,eta,phi;
324 Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0;
325 Float_t rhoback=0.0;
326
327 Float_t ptallback=0.0; //particles without the jet
328 Float_t restarea=accEMCal; //initial area set
329 Bool_t acc=0;
330 Bool_t acc1=0;
331 Float_t rCone=0.4;
332
333 if(nJ==1) {
334 AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
335 jeteta=jettmp->Eta();
336 jetphi=jettmp->Phi();
337 acc=EmcalAcceptance(jeteta,jetphi,rCone);
338 if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
339 cout<<" acc "<<acc<<endl;
340 }
341
342
343 if(nJ>=2) {
344 AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
345 AliAODJet *jettmp1 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
346 jeteta=jettmp->Eta();
347 jetphi=jettmp->Phi();
348 jeteta1=jettmp1->Eta();
349 jetphi1=jettmp1->Phi();
350 acc=EmcalAcceptance(jeteta,jetphi,rCone);
351 acc1=EmcalAcceptance(jeteta1,jetphi1,rCone);
352 if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc;
353 if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc;
354 if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
355
356 cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl;
357
358 }
359
360 // cout<<" nIn = "<<nIn<<endl;
361 Float_t sumpt=0;
362 for(Int_t i=0; i<nIn; i++)
363 { //Unit Array Loop
364 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
365 if(uArray->GetUnitEnergy()>0.){
366
367 pt = uArray->GetUnitEnergy();
368 eta = uArray->GetUnitEta();
369 phi = uArray->GetUnitPhi();
370
371 //cout<<"test emcal acceptance for particles "<<EmcalAcceptance(eta,phi,0.)<<endl;
372
373 Float_t deta=0.0, dphi=0.0, dr=100.0;
374 Float_t deta1=0.0, dphi1=0.0, dr1=100.0;
375
376 //cout<<i<<" pt="<<pt<<" eta="<<eta<<" phi="<<phi<<endl;
377 if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
378 sumpt+=pt;
379 //if(i<30)cout<<i<<" pt unit = "<<pt<<endl;
380
381 if(nJ==1 && acc==1) {
382 deta = eta - jeteta;
383 dphi = phi - jetphi;
384 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
385 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
386 dr = TMath::Sqrt(deta * deta + dphi * dphi);
387 if(dr<=rc)sumpt-=pt;
388 }
389
390 if(nJ>=2) {
391 if(acc==1){
392 deta = eta - jeteta;
393 dphi = phi - jetphi;
394 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
395 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
396 dr = TMath::Sqrt(deta * deta + dphi * dphi);
397 if(dr<=rc)sumpt-=pt;
398 }
399 if(acc1==1){
400 deta1 = eta - jeteta1;
401 dphi1 = phi - jetphi1;
402 if (dphi1 < -TMath::Pi()) dphi1= -dphi1 - 2.0 * TMath::Pi();
403 if (dphi1 > TMath::Pi()) dphi1 = 2.0 * TMath::Pi() - dphi1;
404 dr1 = TMath::Sqrt(deta1 * deta1 + dphi1 * dphi1);
405 if(dr1<=rc)sumpt-=pt;
406 }
407
408 }
409
410 if(dr >= rc && dr1 >=rc) {
411 // particles outside both cones
412
413 //cout<<" out of the cone "<<dr<<" "<<deta<<" deltaeta "<<dphi<<" dphi "<<i<<" particle "<<endl;
414 //cout<<" out of the cone "<<dr1<<" "<<deta1<<" deltaeta1 "<<dphi1<<" dphi1 "<<i<<" particle "<<endl;
415 ptallback+=pt;
416 }
417 }
418 //cout<<" ipart "<<ipart<<" rhointegral "<<rhoback <<endl;
419 }
420 } // End loop on UnitArray
421
422 cout<<"total area left "<<restarea<<endl;
423 cout<<"sumpt="<<sumpt<<endl;
424 // if(acc==1 || acc1==1) rhoback= ptallback/restarea;
425 //else rhoback=ptallback;
426
427 rhoback= ptallback/restarea;
428 cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl;
429
430 return rhoback;
431
432}
433
434
435Double_t AliJetBkg::CalcRho(vector<fastjet::PseudoJet> inputParticles,Double_t rParamBkg,TString method){
436 //calculate rho using the fastjet method
437
438 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
439 Bool_t debug = header->GetDebug(); // debug option
440
441 fastjet::Strategy strategy = header->GetStrategy();
442 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
443 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
444 fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
445
446 // create an object that specifies how we to define the area
447 fastjet::AreaDefinition areaDef;
448 double ghostEtamax = header->GetGhostEtaMax();
449 double ghostArea = header->GetGhostArea();
450 int activeAreaRepeats = header->GetActiveAreaRepeats();
451
452 // now create the object that holds info about ghosts
453
454 if (method.Contains("Charg"))ghostEtamax=0.9;
455
456 fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
457 // and from that get an area definition
458 fastjet::AreaType areaType = header->GetAreaType();
459 areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
460 cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
461 //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
462 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
463 TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
464 comment+= "Jet definition: ";
465 comment+= TString(jetDef.description());
466 // comment+= ". Area definition: ";
467 //comment+= TString(areaDef.description());
468 comment+= ". Strategy adopted by FastJet: ";
469 comment+= TString(clust_seq.strategy_string());
470 header->SetComment(comment);
471 if(debug){
472 cout << "--------------------------------------------------------" << endl;
473 cout << comment << endl;
474 cout << "--------------------------------------------------------" << endl;
475 }
476
477 double ptmin = header->GetPtMin();
478 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
479 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
50c7d9f7 480
db2272d5 481 if (debug) {
482 cout<<"# of BKG jets = "<<jets.size()<<endl;
483 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
484 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]));
485 }
486 }
487
50c7d9f7 488 double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
489
490 if (method.Contains("All")){
491 phiMin = 80.*TMath::Pi()/180+rParamBkg;
492 phiMax = 190.*TMath::Pi()/180-rParamBkg;
50c7d9f7 493 }
494 if (method.Contains("Charg")){
495 phiMin = 0;
496 phiMax = 2*TMath::Pi();
497 }
498 rapMax = ghostEtamax - rParamBkg;
499 rapMin = - ghostEtamax + rParamBkg;
500
501
502 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
503
504
505 Double_t rho=clust_seq.median_pt_per_unit_area(range);
506 // double median, sigma, meanArea;
507 //clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true);
508 //fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea);
509
510 // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);
511
512
513 cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl;
514 return rho;
515}
516
517Float_t AliJetBkg::EtaToTheta(Float_t arg)
518{
519 // return (180./TMath::Pi())*2.*atan(exp(-arg));
520 return 2.*atan(exp(-arg));
521}
522
523
524Double_t AliJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/){
525 //to be implemented--- (pT + Energy in EMCal Acceptance vs Multiplicity)
526 return 1;
527}
528
529
530Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
531
532 Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
533 Float_t deltaphi=110./180.*TMath::Pi();
534 Float_t phicut=deltaphi/2.-radius;
535 Float_t etacut=0.7-radius;
536 //cout<<" eta "<<eta<<" phi "<<phi<<endl;
537 //cout<<etacut<<" "<<phicut<<" "<<meanPhi<<" "<<endl;
538 if(TMath::Abs(eta)<etacut && TMath::Abs(phi-meanPhi)<phicut) return 1;
539 else return 0;
540}