Verbous print-out in debug mode only.
[u/mrichter/AliRoot.git] / JETAN / AliJetBkg.cxx
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():
60   TObject(),
61   fReader(0),
62   fHeader(0),
63   fInputFJ(0)
64 {
65   // Default constructor
66
67 }
68 //______________________________________________________________________
69 AliJetBkg::AliJetBkg(const AliJetBkg& input):
70   TObject(input),
71     fReader(input.fReader),
72     fHeader(input.fHeader),
73     fInputFJ(input.fInputFJ)
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 }
86 //___________________________________________________________________
87 Float_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 //___________________________________________________________________
107 Float_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
130 Float_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 /////////////////////////////////
145 Float_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 ////////////////////////////////////////////////////////////////////////
297 Float_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
435 Double_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); 
480
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   
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;
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
517 Float_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
524 Double_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
530 Bool_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 }