X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliJetBkg.cxx;h=30e7a9b8f2ac46a7e3851262dd3bc88c0428216f;hb=839f7ade62f5ce48da901123aabee41f56a795a3;hp=69185cb7158006735b168537778e23f21bb8cd23;hpb=225f0094b9880be274df542c89d74937d9a5dcd9;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliJetBkg.cxx b/JETAN/AliJetBkg.cxx index 69185cb7158..30e7a9b8f2a 100644 --- a/JETAN/AliJetBkg.cxx +++ b/JETAN/AliJetBkg.cxx @@ -13,594 +13,553 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ + +//-------------------------------------------------- +// Method implementation for background studies and background subtraction with UA1 algorithms +// +// Author: magali.estienne@subatech.in2p3.fr +//------------------------------------------------- + #include -#include -#include #include -#include #include -#include -#include - -#include "AliJetHeader.h" -#include "AliJetReader.h" -#include "AliJetReaderHeader.h" -#include "AliFastJetFinder.h" -#include "AliFastJetHeaderV1.h" -#include "AliJetReaderHeader.h" -#include "AliJetReader.h" -#include "AliJetUnitArray.h" -#include "AliFastJetInput.h" -#include "AliESDEvent.h" - - -#include "fastjet/PseudoJet.hh" -#include "fastjet/ClusterSequenceArea.hh" -#include "fastjet/AreaDefinition.hh" -#include "fastjet/JetDefinition.hh" -// get info on how fastjet was configured -#include "fastjet/config.h" - -#ifdef ENABLE_PLUGIN_SISCONE -#include "fastjet/SISConePlugin.hh" -#endif - -#include // needed for internal io -#include -#include +#include +#include -using namespace std; +#include "AliAODJetEventBackground.h" +#include "AliUA1JetHeaderV1.h" +#include "AliJetCalTrk.h" #include "AliJetBkg.h" + +using namespace std; + ClassImp(AliJetBkg) //////////////////////////////////////////////////////////////////////// AliJetBkg::AliJetBkg(): TObject(), - fReader(0), - fHeader(0), - fInputFJ(0) + fEvent(0x0), + fHeader(0x0), + fDebug(0), + fhEtBackg(0x0), + fhAreaBackg(0x0) { // Default constructor - + for(int i = 0;i < kMaxJets;i++){ + fhAreaJet[i] = fhEtJet[i] = 0; + } } -//______________________________________________________________________ + +//---------------------------------------------------------------- AliJetBkg::AliJetBkg(const AliJetBkg& input): TObject(input), - fReader(input.fReader), - fHeader(input.fHeader), - fInputFJ(input.fInputFJ) + fEvent(input.fEvent), + fHeader(input.fHeader), + fDebug(input.fDebug), + fhEtBackg(input.fhEtBackg), + fhAreaBackg(input.fhAreaBackg) { // copy constructor - -} -//______________________________________________________________________ -AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){ - // Assignment operator. - this->~AliJetBkg(); - new(this) AliJetBkg(source); - return *this; + for(int i = 0;i < kMaxJets;i++){ + fhAreaJet[i] = input.fhAreaJet[i]; + fhEtJet[i] = input.fhEtJet[i]; + } } +//---------------------------------------------------------------- +AliJetBkg::~AliJetBkg() +{ + // Destructor + if(fhEtBackg) delete fhEtBackg; + if(fhAreaBackg) delete fhAreaBackg; + for(int i = 0;i < kMaxJets;i++){ + if(fhAreaJet[i]) delete fhAreaJet[i]; + if(fhEtJet[i]) delete fhEtJet[i]; + } -//___________________________________________________________________ - 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() =========== "< inputParticles=fInputFJ->GetInputParticles(); - - //cout<<"printing inputParticles for BKG "<GetCalTrkTrack(id)->GetCutFlag() == 1) + return kTRUE; + else return kFALSE; - +} - double rParamBkg = header->GetRparamBkg(); //Radius for background calculation +//---------------------------------------------------------------- +Bool_t AliJetBkg::SignalCutPass(Int_t id, Int_t nTracks) +{ + // Check if track or cell passes the cut flag + if(id < nTracks && fEvent->GetCalTrkTrack(id)->GetSignalFlag() == 1) + return kTRUE; + else return kFALSE; - 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() =========== "< inputParticles=fInputFJ->GetInputParticles(); - - //cout<<"printing inputParticles for BKG "<GetRparamBkg(); //Radius for background calculation - Double_t medianb,sigmab,meanareab; - CalcRhoWoHardest(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All"); - rho=medianb; - sigma=sigmab; - meanarea=meanareab; +//---------------------------------------------------------------- +Float_t AliJetBkg::CalcJetAreaEtaCut(Float_t radius, Float_t etaJet) +{ + // Calculate jet area taking into account an acceptance cut in eta + AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; + Float_t detamax = etaJet + radius; + Float_t detamin = etaJet - radius; + Float_t accmax = 0.0; Float_t accmin = 0.0; + if(detamax > header->GetLegoEtaMax()){ // sector outside etamax + Float_t h = header->GetLegoEtaMax() - etaJet; + accmax = radius*radius*TMath::ACos(h/radius) - h*TMath::Sqrt(radius*radius - h*h); } - -//____________________________________________________________________ - - void AliJetBkg::CalcRhob(Double_t& median,Double_t& - sigma,Double_t& - meanarea,vector 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; + if(detamin < header->GetLegoEtaMin()){ // sector outside etamin + Float_t h = header->GetLegoEtaMax() + etaJet; + accmin = radius*radius*TMath::ACos(h/radius) - h*TMath::Sqrt(radius*radius - h*h); } - - vector inclusiveJets = clust_seq.inclusive_jets(0.); - vector jets = sorted_by_pt(inclusiveJets); - - - //cout<<"# of BKG jets = "< 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 +} - +//---------------------------------------------------------------- +void AliJetBkg::CalcJetAndBckgAreaEtaCut(Bool_t calcOutsideArea, Float_t radius, Int_t nJets, const Float_t* etaJet, Float_t* &areaJet, Float_t &areaOut) +{ + // Calculate jet and bacground areas taking into account an acceptance cut in eta - 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="<SetComment(comment); - if(debug){ - cout << "--------------------------------------------------------" << endl; - cout << comment << endl; - cout << "--------------------------------------------------------" << endl; + AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; + areaOut = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax() - header->GetLegoPhiMin()); + for(Int_t k=0; k inclusiveJets = clust_seq.inclusive_jets(0.); - vector jets = sorted_by_pt(inclusiveJets); - vector jets2=sorted_by_pt(inclusiveJets); - if(jets2.size()>=2) jets2.erase(jets2.begin(),jets2.begin()+1); - - - - 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(){ - - AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; - Bool_t debug = header->GetDebug(); // debug option - - if(debug)cout<<"=============== AliJetBkg::BkgFastJet() =========== "< inputParticles=fInputFJ->GetInputParticles(); - - if(debug)cout<<"printing inputParticles for BKG "<GetRparamBkg(); //Radius for background calculation - Double_t rho=CalcRho(inputParticles,rParamBkg,"All"); - if(debug)cout<<"-------- rho (from all part)="<GetDebug(); // debug option - - if(debug)cout<<"=============== AliJetBkg::BkgChargedFastJet() =========== "< inputParticlesCharged=fInputFJ->GetInputParticlesCh(); +//---------------------------------------------------------------- +void AliJetBkg::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, Float_t&sigmaN, + const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, + Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, + Float_t* etsigJet, Int_t* multJetT, Int_t* multJetC, Int_t* multJet, + Int_t* injet, Float_t* &areaJet) +{ + // + // Background subtraction using cone method but without correction in dE/deta distribution + // Cases to take into account the EMCal geometry are included + // - if(debug)cout<<"printing CHARGED inputParticles for BKG "<GetRparam(); + AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; + //calculate energy inside and outside cones + fDebug = header->GetDebug(); + Float_t rc = header->GetRadius(); + Float_t etOut = 0; + // Get number of tracks from EventCalTrk + Int_t nTracks = fEvent->GetNCalTrkTracks(); + + Float_t etIn[kMaxJets] = {0}; + Float_t areaOut = 0.; + + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + if(jpart < nTracks) multJetT[ijet]++; + else multJetC[ijet]++; + multJet[ijet]++; + injet[jpart] = ijet; + if(PtCutPass(jpart,nTracks)){ // pt cut + etIn[ijet] += ptT[jpart]; + if(SignalCutPass(jpart,nTracks)) + etsigJet[ijet]+= ptT[jpart]; + } + break; + } + }// end jets loop + + if((injet[jpart] == -1) && + (PtCutPass(jpart,nTracks))){ + etOut += ptT[jpart]; // particle outside cones and pt cut + } + } //end particle loop + + // Calculate jet and background areas + Bool_t calcAreaOut = kTRUE; + CalcJetAndBckgAreaEtaCut(calcAreaOut,rc, nJ, etaJet, areaJet, areaOut); + + //subtract background using area method + for(Int_t ljet=0; ljetGetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()); + etbgTotalN = etOut*areaT/areaOut; + + // estimate standard deviation of background + Int_t count = 0; + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + if((injet[jpart] == -1) && + (PtCutPass(jpart,nTracks))){ + sigmaN += etbgTotalN/areaT - ptT[jpart]; + // To be checked (Division by jet area to obtain standard deviation of rho ?) + + count=count+1; + } + } + if (count>0) + sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count); - if(debug)cout<<"-------- rho (from CHARGED part)="<GetDebug(); // debug option - - if(debug)cout<<"==============AliJetBkg::BkgStat()============="<GetESD()->GetNumberOfTracks(); - Int_t nTracks= 0; - TF1 fun("fun",BkgFunction,0,800,1); - Double_t enTot=fun.Eval(nTracks); - Double_t accEMCal=2*0.7*110./180*TMath::Pi();//2.68 area of EMCal - return enTot/accEMCal; - -} - -//////////////////////////////////////////////////////////////////////// -Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets) -{ + // Cases to take into account the EMCal geometry are included + // - // 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. - - AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; - Bool_t debug = header->GetDebug(); // debug option - - if(debug)cout<<"==============AliJetBkg::SubtractFastJetBackgCone()============="<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... - if(debug)cout<<"nJets: "<GetBackgStat(); // pre-calculated background - - //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; } + //calculate energy inside + Float_t rc= header->GetRadius(); + Float_t etIn[kMaxJets] = {0.0}; + // Get number of tracks from EventCalTrk + Int_t nTracks = fEvent->GetNCalTrkTracks(); + Float_t areaOut = 0.; + + for(Int_t jpart = 0; jpart < nIn; jpart++) + { // loop for all particles in array + + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + if(jpart < nTracks) multJetT[ijet]++; + else multJetC[ijet]++; + multJet[ijet]++; + injet[jpart] = ijet; + + if(PtCutPass(jpart,nTracks)){ // pt cut + etIn[ijet] += ptT[jpart]; + if(SignalCutPass(jpart,nTracks)) + etsigJet[ijet]+= ptT[jpart]; + } + break; + } + }// end jets loop + } //end particle loop - // Information extracted from fUnitArray - // load input vectors and calculate total energy in array - Float_t pt,eta,phi; - Float_t jeteta = 0,jetphi = 0,jeteta1 = 0, jetphi1 = 0; - Float_t rhoback=0.0; - - Float_t ptallback=0.0; //particles without the jet - Float_t restarea=accEMCal; //initial area set - Bool_t acc=0; - Bool_t acc1=0; - Float_t rCone=0.4; - - if(nJ==1) { - 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; - if(debug)cout<<" acc "<=2) { - AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0)); - AliAODJet *jettmp1 = (AliAODJet*)(fAODJets->At(1)); - jeteta=jettmp->Eta(); - jetphi=jettmp->Phi(); - jeteta1=jettmp1->Eta(); - jetphi1=jettmp1->Phi(); - acc=EmcalAcceptance(jeteta,jetphi,rCone); - acc1=EmcalAcceptance(jeteta1,jetphi1,rCone); - if(acc1==1 && acc==1)restarea= accEMCal-2*TMath::Pi()*rc*rc; - if(acc1==1 && acc==0)restarea= accEMCal-TMath::Pi()*rc*rc; - if(acc1==0 && acc==1)restarea= accEMCal-TMath::Pi()*rc*rc; - - if(debug)cout<<" acc1="<0)sigmaN=TMath::Sqrt(TMath::Abs(sigmaN)/count); - return rhoback; - } - -Double_t AliJetBkg::CalcRho(vector 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 - - if (method.Contains("Charg"))ghostEtamax=0.9; - - 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); - if(debug)cout<<"rParamBkg="<SetComment(comment); - if(debug){ - cout << "--------------------------------------------------------" << endl; - cout << comment << endl; - cout << "--------------------------------------------------------" << endl; +//---------------------------------------------------------------- +void AliJetBkg::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t&sigmaN, + const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, Float_t* etJet, + const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet, + Int_t* multJetT, Int_t* multJetC, Int_t* multJet, Int_t* injet, Float_t* &/*areaJet*/) +{ + // + // Cone background subtraction method taking into acount dEt/deta distribution + // Cases to take into account the EMCal geometry are not included + // + + AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader; + //general + Float_t rc= header->GetRadius(); + Float_t etamax = header->GetLegoEtaMax(); + Float_t etamin = header->GetLegoEtaMin(); + Int_t ndiv = 100; + // Get number of tracks from EventCalTrk + Int_t nTracks = fEvent->GetNCalTrkTracks(); + + // jet energy and area arrays + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + for(Int_t mjet=0; mjetReset(); + fhAreaJet[mjet]->Reset(); } - - double ptmin = header->GetPtMin(); - vector inclusiveJets = clust_seq.inclusive_jets(ptmin); - vector jets = sorted_by_pt(inclusiveJets); - - if (debug) { - cout<<"# of BKG jets = "<Reset(); + if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); + fhAreaBackg->Reset(); + TH1::AddDirectory(oldStatus); + + //fill energies + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + if(jpart < nTracks) multJetT[ijet]++; + else multJetC[ijet]++; + multJet[ijet]++; + injet[jpart] = ijet; + + if(PtCutPass(jpart,nTracks)){ // pt cut + fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone + if(SignalCutPass(jpart,nTracks)) + etsigJet[ijet]+= ptT[jpart]; + } + break; + } + }// end jets loop + + if((injet[jpart] == -1) && + (PtCutPass(jpart,nTracks) == 1)) + fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + } //end particle loop + + //calc areas + Float_t eta0 = etamin; + Float_t etaw = (etamax - etamin)/((Float_t)ndiv); + Float_t eta1 = eta0 + etaw; + for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins + Float_t etac = eta0 + etaw/2.0; + Float_t areabg = etaw*2.0*TMath::Pi(); + for(Int_t ijet=0; ijet rc && deta1 < rc){ + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + areaj = acc1; + } + if(deta0 < rc && deta1 > rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + areaj = acc0; + } + if(deta0 < rc && deta1 < rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + if(eta1Fill(etac,areaj); + areabg = areabg - areaj; + } // end jets loop + fhAreaBackg->Fill(etac,areabg); + eta0 = eta1; + eta1 = eta1 + etaw; + } // end loop for all eta bins + + //subtract background + for(Int_t kjet=0; kjetGetBinContent(bin)){ + Float_t areab = fhAreaBackg->GetBinContent(bin); + Float_t etb = fhEtBackg->GetBinContent(bin); + Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab; + etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction } + } } - - 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; - } - if (method.Contains("Charg")){ - phiMin = 0; - phiMax = 2*TMath::Pi(); + // calc background total + Double_t etOut = fhEtBackg->Integral(); + Double_t areaOut = fhAreaBackg->Integral(); + Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()); + etbgTotalN = etOut*areaT/areaOut; + + Int_t count=0; + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + if((injet[jpart] == -1) && + (PtCutPass(jpart,nTracks))){ + sigmaN += etbgTotalN/areaT - ptT[jpart]; + count=count+1; + } } - rapMax = ghostEtamax - rParamBkg; - rapMin = - ghostEtamax + rParamBkg; - - - fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); - - - Double_t rho=clust_seq.median_pt_per_unit_area(range); - // double median, sigma, meanArea; - //clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, median, sigma, meanArea, true); - //fastjet::ActiveAreaSpec area_spec(ghostEtamax,activeAreaRepeats,ghostArea); - - // fastjet::ClusterSequenceActiveArea clust_seq_bkg(inputParticles, jetDef,area_spec); - - - if(debug)cout<<"bkg in R="<GetBackgCutRatio(); + + //general + Float_t rc= header->GetRadius(); + Float_t etamax = header->GetLegoEtaMax(); + Float_t etamin = header->GetLegoEtaMin(); + Int_t ndiv = 100; + // Get number of tracks from EventCalTrk + Int_t nTracks = fEvent->GetNCalTrkTracks(); + + // jet energy and area arrays + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + for(Int_t mjet=0; mjetReset(); + fhAreaJet[mjet]->Reset(); + } + // background energy and area + if(!fhEtBackg)fhEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); + fhEtBackg->Reset(); + if(!fhAreaBackg) fhAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); + fhAreaBackg->Reset(); + TH1::AddDirectory(oldStatus); + + //fill energies + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + if(jpart < nTracks) multJetT[ijet]++; + else multJetC[ijet]++; + multJet[ijet]++; + injet[jpart] = ijet; + + if(PtCutPass(jpart,nTracks)){ // pt cut + fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone + if(SignalCutPass(jpart,nTracks)) + etsigJet[ijet]+= ptT[jpart]; + } + break; + } + }// end jets loop + if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + } //end particle loop + + //calc areas + Float_t eta0 = etamin; + Float_t etaw = (etamax - etamin)/((Float_t)ndiv); + Float_t eta1 = eta0 + etaw; + for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins + Float_t etac = eta0 + etaw/2.0; + Float_t areabg = etaw*2.0*TMath::Pi(); + for(Int_t ijet=0; ijet rc && deta1 < rc){ + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + areaj = acc1; + } + if(deta0 < rc && deta1 > rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + areaj = acc0; + } + if(deta0 < rc && deta1 < rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + if(eta1Fill(etac,areaj); + areabg = areabg - areaj; + } // end jets loop + fhAreaBackg->Fill(etac,areabg); + eta0 = eta1; + eta1 = eta1 + etaw; + } // end loop for all eta bins + + //subtract background + for(Int_t kjet=0; kjetGetBinContent(bin)){ + Float_t areab = fhAreaBackg->GetBinContent(bin); + Float_t etb = fhEtBackg->GetBinContent(bin); + Float_t areaR = (fhAreaJet[kjet]->GetBinContent(bin))/areab; + etJet[kjet] = etJet[kjet] + ((fhEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction + } + } + } -Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){ + // calc background total + Double_t etOut = fhEtBackg->Integral(); + Double_t areaOut = fhAreaBackg->Integral(); + Float_t areaT = (header->GetLegoEtaMax()-header->GetLegoEtaMin())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()); + etbgTotalN = etOut*areaT/areaOut; - 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; - Float_t etacut=0.7-radius; - //cout<<" eta "<