return *this;
}
+
+
//___________________________________________________________________
-Float_t AliJetBkg::BkgFastJet(){
+ void AliJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma,
+Double_t& meanarea){
+
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+ if(debug)cout<<"=============== AliJetBkg::BkgFastJetb() =========== "<<endl;
+ vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+
+ //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
+
+
+
+
+ double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
+
+ Double_t medianb,sigmab,meanareab;
+ CalcRhob(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
+ rho=medianb;
+ sigma=sigmab;
+ meanarea=meanareab;
+
+}
+//_________________________________________________________________
+
+ void AliJetBkg::BkgFastJetWoHardest(Double_t& rho,Double_t& sigma,
+Double_t& meanarea){
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+ if(debug)cout<<"=============== AliJetBkg::BkgWoHardest() =========== "<<endl;
+ vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+
+ //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
+ double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
+ Double_t medianb,sigmab,meanareab;
+ CalcRhoWoHardest(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
+ rho=medianb;
+ sigma=sigmab;
+ meanarea=meanareab;
+
+ }
+
+//____________________________________________________________________
+
+ void AliJetBkg::CalcRhob(Double_t& median,Double_t&
+ sigma,Double_t&
+ meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
+ rParamBkg,TString method){
+ //calculate rho using the fastjet method
+
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+
+ fastjet::Strategy strategy = header->GetStrategy();
+ fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
+ fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
+ fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
+
+ // create an object that specifies how we to define the area
+ fastjet::AreaDefinition areaDef;
+ double ghostEtamax = header->GetGhostEtaMax();
+ double ghostArea = header->GetGhostArea();
+ int activeAreaRepeats = header->GetActiveAreaRepeats();
+
+ // now create the object that holds info about ghosts
+
+
+
+ fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
+ // and from that get an area definition
+ fastjet::AreaType areaType = header->GetAreaType();
+ areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
+
+ //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. ";
+ comment+= "Jet definition: ";
+ comment+= TString(jetDef.description());
+ // comment+= ". Area definition: ";
+ //comment+= TString(areaDef.description());
+ comment+= ". Strategy adopted by FastJet: ";
+ comment+= TString(clust_seq.strategy_string());
+ comment+= Form("Method: %s",method.Data());
+ header->SetComment(comment);
+ if(debug){
+ cout << "--------------------------------------------------------" << endl;
+ cout << comment << endl;
+ cout << "--------------------------------------------------------" << endl;
+ }
+
+ vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
+ 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 phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+
+
+ phiMin = 0;
+ phiMax = 2*TMath::Pi();
+
+ rapMax = ghostEtamax - rParamBkg;
+ rapMin = - ghostEtamax + rParamBkg;
+
+
+ fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
+
+
+ double medianb, sigmab, meanareab;
+ clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, medianb, sigmab, meanareab, false);
+ median=medianb;
+ sigma=sigmab;
+ meanarea=meanareab;
+
+
+ }
+
+//____________________________________________________________________
+
+
+ void AliJetBkg::CalcRhoWoHardest(Double_t& median,Double_t&
+sigma,Double_t&
+meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t
+rParamBkg,TString method){
+ //calculate rho using the fastjet method
+
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
+
+ fastjet::Strategy strategy = header->GetStrategy();
+ fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
+ fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
+ fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
+
+ // create an object that specifies how we to define the area
+ fastjet::AreaDefinition areaDef;
+ double ghostEtamax = header->GetGhostEtaMax();
+ double ghostArea = header->GetGhostArea();
+ int activeAreaRepeats = header->GetActiveAreaRepeats();
+
+ // now create the object that holds info about ghosts
+
+
+
+ fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
+ // 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;
+ //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. ";
+ comment+= "Jet definition: ";
+ comment+= TString(jetDef.description());
+ // comment+= ". Area definition: ";
+ //comment+= TString(areaDef.description());
+ comment+= ". Strategy adopted by FastJet: ";
+ comment+= TString(clust_seq.strategy_string());
+ comment+= Form("Method: %s",method.Data());
+ header->SetComment(comment);
+ if(debug){
+ cout << "--------------------------------------------------------" << endl;
+ cout << comment << endl;
+ cout << "--------------------------------------------------------" << endl;
+ }
+
+ vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
+ vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
+ vector<fastjet::PseudoJet> jets2=sorted_by_pt(inclusiveJets);
+ if(jets2.size()>=2) jets2.erase(jets2.begin(),jets2.begin()+1);
+
+
- cout<<"=============== AliJetBkg::BkgFastJet() =========== "<<endl;
+ double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+
+
+ phiMin = 0;
+ phiMax = 2*TMath::Pi();
+
+ rapMax = ghostEtamax - rParamBkg;
+ rapMin = - ghostEtamax + rParamBkg;
+
+
+ fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
+
+
+ double medianb, sigmab, meanareab;
+ clust_seq.get_median_rho_and_sigma(jets2, range, false, medianb, sigmab,
+meanareab, false);
+ median=medianb;
+ sigma=sigmab;
+ meanarea=meanareab;
+
+
+ }
+
+
+//___________________________________________________________________
+Float_t AliJetBkg::BkgFastJet(){
+// Return background
+ 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(){
+// Background for charged jets
+ AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+ Bool_t debug = header->GetDebug(); // debug option
- cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "<<endl;
+ 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");
- cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
+ if(debug)cout<<"-------- rho (from CHARGED part)="<<rho<<endl;
return rho;
}
{
//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
Float_t rCone=0.4;
if(nJ==1) {
- AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
+ AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
jeteta=jettmp->Eta();
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(nJ>=2) {
- AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
- AliAODJet *jettmp1 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
+ AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
+ AliAODJet *jettmp1 = (AliAODJet*)(fAODJets->At(1));
jeteta=jettmp->Eta();
jetphi=jettmp->Phi();
jeteta1=jettmp1->Eta();
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. ";
// 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;
}
}
-Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
-
+Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius) const{
+// Apply emcal acceptance cuts
Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
Float_t deltaphi=110./180.*TMath::Pi();
Float_t phicut=deltaphi/2.-radius;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff