//___________________________________________________________________
Float_t AliJetBkg::BkgFastJet(){
- cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+
+ if(debug)cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
- cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
+ if(debug)cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
for(UInt_t i=0;i<inputParticles.size();i++){
// cout<<" "<<inputParticles[i].px()<<" "<<inputParticles[i].py()<<" "<<inputParticles[i].pz()<<endl;
}
- AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
Double_t rho=CalcRho(inputParticles,rParamBkg,"All");
- cout<<"-------- rho (from all part)="<<rho<<endl;
+ if(debug)cout<<"-------- rho (from all part)="<<rho<<endl;
return rho;
}
//___________________________________________________________________
Float_t AliJetBkg::BkgChargedFastJet(){
- cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+
+ if(debug)cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
vector<fastjet::PseudoJet> inputParticlesCharged=fInputFJ->GetInputParticlesCh();
- cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
+ if(debug)cout<<"printing CHARGED inputParticles for BKG "<<inputParticlesCharged.size()<<endl;
for(UInt_t i=0;i<inputParticlesCharged.size();i++){
// cout<<" "<<inputParticlesCharged[i].px()<<" "<<inputParticlesCharged[i].py()<<" "<<inputParticlesCharged[i].pz()<<endl;
}
- AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+
double rParam = header->GetRparam();
Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
{
//background subtraction using statistical method
- cout<<"==============AliJetBkg::BkgStat()============="<<endl;
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+
+ if(debug)cout<<"==============AliJetBkg::BkgStat()============="<<endl;
//TO BE IMPLEMENTED
// Int_t nTracks= fReader->GetESD()->GetNumberOfTracks();
Int_t nTracks= 0;
return enTot/accEMCal;
}
-/////////////////////////////////
-Float_t AliJetBkg::BkgRemoveJetLeading(TClonesArray* fAODJets)
-{
- // Remove the particles of the
- // two largest jets using the track references stored in the AODJet from the estimation of new rho.
-
- cout<<"==============AliJetBkg::BkgRemoveJetLeading()============="<<endl;
-
-
- // check if we are reading AOD jets
- TRefArray *refs = 0;
- Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetESDReader");
- if (fromAod) { refs = fReader->GetReferences(); }
-
- //Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
- Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
-
- Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
- cout<<"nJets: "<<nJ<<endl;
-
-
- //begin unit array
- TClonesArray* fUnit = fReader->GetUnitArray();
- if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return -99; }
-
- Int_t nIn = fUnit->GetEntries();
- if(nIn == 0) { cout << "entries = 0 ; Event empty !!!" << endl ; return -99; }
-
- Float_t rhoback=0.0;
- Float_t jetarea1=0.0,jetarea2=0.0;
-
- Int_t particlejet1=-99;
- Int_t particlejet2=-99;
- TRefArray *refarray1 = 0;
- TRefArray *refarray2 = 0;
- Int_t nJettracks1 = 0, nJettracks2 = 0;
- Int_t acc=0,acc1=0;
- AliAODJet *jet1;
- AliAODJet *jet2;
-
- if(nJ==1){
- jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
- jetarea1=jet1->EffectiveAreaCharged();
- Float_t jetPhi=jet1->Phi();
- Float_t jetEta=jet1->Eta();
- if(jetPhi>1.396 && jetPhi<3.316 && jetEta>-0.7 && jetEta<0.7)acc=1;
- refarray1=jet1->GetRefTracks();
- nJettracks1=refarray1->GetEntries();
- cout<<"nJ = 1, acc="<<acc<<" jetarea1="<<jetarea1<<endl;
-
- }
-
- if(nJ>=2){
- jet1 = dynamic_cast<AliAODJet*>(fAODJets->At(0));
- jetarea1=jet1->EffectiveAreaCharged();
- Float_t jetPhi1=jet1->Phi();
- Float_t jetEta1=jet1->Eta();
- if(jetPhi1>1.396 && jetPhi1<3.316 && jetEta1>-0.7 && jetEta1<0.7)acc=1;
- refarray1=jet1->GetRefTracks();
- nJettracks1=refarray1->GetEntries();
- cout<<"npart = "<<nJettracks1<<endl;
-
- jet2 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
- jetarea2=jet2->EffectiveAreaCharged();
- Float_t jetPhi2=jet2->Phi();
- Float_t jetEta2=jet2->Eta();
- if(jetPhi2>1.396 && jetPhi2<3.316 && jetEta2>-0.7 && jetEta2<0.7)acc1=1;
- refarray2=jet2->GetRefTracks();
- nJettracks2=refarray2->GetEntries();
- cout<<"nJ = "<<nJ<<", acc="<<acc<<" acc1="<<acc1<<" jetarea1="<<jetarea1<<" jetarea2="<<jetarea2<<endl;
- }
-
-
-
- // cout<<" nIn = "<<nIn<<endl;
- Float_t sumPt=0;
- Float_t eta,phi,pt;
- Int_t ipart=0;
-
- for(Int_t i=0; i<nIn; i++)
- { //Unit Array Loop
- AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
-
- if(uArray->GetUnitEnergy()>0.){
- eta = uArray->GetUnitEta();
- phi = uArray->GetUnitPhi();
- pt = uArray->GetUnitEnergy();
- // cout<<"ipart = "<<ipart<<" eta="<<eta<<" phi="<<phi<<endl;
- if(phi>1.396 && phi<3.316 && eta>-0.7 && eta<0.7){
- //cout<<sumPt<<endl;
- sumPt+=pt;
-
-
- if(nJ==1 && acc==1){
- for(Int_t ii=0; ii<nJettracks1;ii++){
-
- particlejet1 = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
-
- if(ipart==particlejet1) {
- sumPt-=pt;
- }
- }
- }
-
-
- if(nJ>=2){
-
- //first jet
- if(acc==1){
- for(Int_t ii=0; ii<nJettracks1;ii++){
- particlejet1 = ((AliJetUnitArray*)refarray1->At(ii))->GetUnitTrackID();
-
- //cout<<"uArr loop = "<<i<<" ipart in uArr (1/2)="<<ipart<<" part in jet="<<ii<<" partID="<<particlejet1<<" sumPt="<<sumPt<<endl;
- if(ipart==particlejet1) {
- sumPt-=pt;
- }
- }
- }
- if(acc1==1){
- //second jet
- for(Int_t ii=0; ii<nJettracks2;ii++){
- particlejet2 = ((AliJetUnitArray*)refarray2->At(ii))->GetUnitTrackID();
- //cout<<"uArr loop = "<<i<<" ipart in uArr (2/2)="<<ipart<<" part in jet="<<ii<<" partID="<<particlejet2<<" sumPt="<<sumPt<<endl;
- if(ipart==particlejet2) {
- sumPt-=pt;
- }
- }
- }
-
-
- }
-
- }//if phi,eta
- ipart++;
- }//end if energy
- }// end unit array loop
-
-
- Float_t areasum=areasum=accEMCal-acc*jetarea1-acc1*jetarea2;
- cout<<"pt sum "<<sumPt<<" area "<<areasum<<endl;
-
- if(nJ>0) rhoback=sumPt/areasum;
- else rhoback=0.;
- cout<<" rho from leading jet paricle array removed "<<rhoback<<endl;
-
- return rhoback;
-
-}
-
-
-
////////////////////////////////////////////////////////////////////////
Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
// Cone background subtraction method applied on the fastjet: REmove the particles of the
// two largest jets with the given R from the estimation of new rho.
- cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+
+ if(debug)cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<<endl;
+
Float_t rc= header->GetRparam();
//Hard wired Calorimeter area (get it later from the AliJetReaderHeader.h)
Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal
Int_t nJ=fAODJets->GetEntries(); //this must be the # of jets...
- cout<<"nJets: "<<nJ<<endl;
+ if(debug)cout<<"nJets: "<<nJ<<endl;
//begin unit array
jetphi=jettmp->Phi();
acc=EmcalAcceptance(jeteta,jetphi,rCone);
if(acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
- cout<<" acc "<<acc<<endl;
+ if(debug)cout<<" acc "<<acc<<endl;
}
if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc;
if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc;
- cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl;
+ if(debug)cout<<" acc1="<<acc<<" acc2="<<acc1<<" restarea="<<restarea<<endl;
}
}
} // End loop on UnitArray
- cout<<"total area left "<<restarea<<endl;
- cout<<"sumpt="<<sumpt<<endl;
+ if(debug)cout<<"total area left "<<restarea<<endl;
+ if(debug)cout<<"sumpt="<<sumpt<<endl;
// if(acc==1 || acc1==1) rhoback= ptallback/restarea;
//else rhoback=ptallback;
rhoback= ptallback/restarea;
- cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl;
+ if(debug)cout<<"rhoback "<<rhoback<<" "<<nJ<<" "<<endl;
return rhoback;
// and from that get an area definition
fastjet::AreaType areaType = header->GetAreaType();
areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
- cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
+ if(debug)cout<<"rParamBkg="<<rParamBkg<<" ghostEtamax="<<ghostEtamax<<" ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
//fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
double ptmin = header->GetPtMin();
vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
- cout<<"# of BKG jets = "<<jets.size()<<endl;
- for (size_t j = 0; j < jets.size(); j++) { // loop for jets
-
- 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]));
- }
- // double phiMax = header->GetPhiMax();
- //double phiMin = header->GetPhiMin();
+ if (debug) {
+ cout<<"# of BKG jets = "<<jets.size()<<endl;
+ for (size_t j = 0; j < jets.size(); j++) { // loop for jets
+ 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]));
+ }
+ }
+
double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
if (method.Contains("All")){
phiMin = 80.*TMath::Pi()/180+rParamBkg;
phiMax = 190.*TMath::Pi()/180-rParamBkg;
- //phiMin = 0;
- //phiMax = 2*TMath::Pi();
}
if (method.Contains("Charg")){
phiMin = 0;
// fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec);
- cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl;
+ if(debug)cout<<"bkg in R="<<rParamBkg<<" : "<<rho<<" range: Rap="<<rapMin<<","<<rapMax<<" -- phi="<<phiMin<<","<<phiMax<<endl;
return rho;
}