X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliUA1JetFinderV1.cxx;h=d8f74f9a49cf3dd1d38cc47fbc7da1bf63f60cca;hb=1c13a3a0af41c8a8a0069dc356daa2f27e57d3cf;hp=dbcb9d5a70981ccaa598143dedafe35e5d5c61db;hpb=970a3bbcc28a0f316267735ad4d8b6a887b85d57;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliUA1JetFinderV1.cxx b/JETAN/AliUA1JetFinderV1.cxx index dbcb9d5a709..d8f74f9a49c 100644 --- a/JETAN/AliUA1JetFinderV1.cxx +++ b/JETAN/AliUA1JetFinderV1.cxx @@ -13,6 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ //--------------------------------------------------------------------- // UA1 Cone Algorithm Jet finder @@ -21,28 +22,38 @@ // (version in c++) //--------------------------------------------------------------------- -#include +#include +#include #include #include #include -#include +#include + #include "AliUA1JetFinderV1.h" #include "AliUA1JetHeaderV1.h" #include "AliJetReaderHeader.h" #include "AliJetReader.h" -#include "AliJet.h" +#include "AliJetHeader.h" -ClassImp(AliUA1JetFinderV1) +#include "AliAODJet.h" +#include "AliLog.h" -//////////////////////////////////////////////////////////////////////// -AliUA1JetFinderV1::AliUA1JetFinderV1() +ClassImp(AliUA1JetFinderV1) + +///////////////////////////////////////////////////////////////////// +AliUA1JetFinderV1::AliUA1JetFinderV1() : + AliJetFinder(), + fLego(0), + fhEtBackg(0), + fhAreaBackg(0) { // Constructor - fHeader = 0x0; - fLego = 0x0; + for(int i = 0;i < kMaxJets;i++){ + fhAreaJet[i] = fhEtJet[i] = 0; + } } //////////////////////////////////////////////////////////////////////// @@ -51,6 +62,18 @@ AliUA1JetFinderV1::~AliUA1JetFinderV1() { // destructor + delete fLego; + fLego = 0; + if(fhEtBackg)delete fhEtBackg; + fhEtBackg = 0; + if( fhAreaBackg) delete fhAreaBackg; + fhAreaBackg = 0; + for(int i = 0;i < kMaxJets;i++){ + if(fhAreaJet[i])delete fhAreaJet[i]; + if(fhEtJet[i]) delete fhEtJet[i]; + fhAreaJet[i] = fhEtJet[i] = 0; + } + } //////////////////////////////////////////////////////////////////////// @@ -74,58 +97,74 @@ void AliUA1JetFinderV1::FindJets() AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; TClonesArray *lvArray = fReader->GetMomentumArray(); Int_t nIn = lvArray->GetEntries(); - if (nIn == 0) return; + if (nIn <= 0) return; // local arrays for input - Float_t* ptT = new Float_t[nIn]; - Float_t* etaT = new Float_t[nIn]; - Float_t* phiT = new Float_t[nIn]; + // ToDo: check memory fragmentation, maybe better to + // define them globally and resize as needed + // Fragementation should be worse for low mult... + Float_t* ptT = new Float_t[nIn]; + Float_t* etaT = new Float_t[nIn]; + Float_t* phiT = new Float_t[nIn]; Int_t* injet = new Int_t[nIn]; + memset(ptT,0,sizeof(Float_t)*nIn); + memset(etaT,0,sizeof(Float_t)*nIn); + memset(phiT,0,sizeof(Float_t)*nIn); + + + // load input vectors and calculate total energy in array + //total energy in array Float_t etbgTotal = 0.0; - TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0); + Float_t npart = 0; + Float_t etbg2 = 0; - // load input vectors and calculate total energy in array for (Int_t i = 0; i < nIn; i++){ TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); ptT[i] = lv->Pt(); etaT[i] = lv->Eta(); phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi()); if (fReader->GetCutFlag(i) != 1) continue; - fLego->Fill(etaT[i], phiT[i], ptT[i]); - hPtTotal->Fill(ptT[i]); + fLego ->Fill(etaT[i], phiT[i], ptT[i]); + npart += 1; etbgTotal+= ptT[i]; + etbg2 += ptT[i]*ptT[i]; } - fJets->SetNinput(nIn); // calculate total energy and fluctuation in map - Double_t meanpt = hPtTotal->GetMean(); - Double_t ptRMS = hPtTotal->GetRMS(); - Double_t npart = hPtTotal->GetEntries(); + Double_t meanpt = 0; + Double_t ptRMS = 0; + if(npart>0){ + meanpt = etbgTotal/npart; + etbg2 = etbg2/npart; + if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities + ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt); + } + } Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS); // arrays to hold jets - Float_t* etaJet = new Float_t[30]; - Float_t* phiJet = new Float_t[30]; - Float_t* etJet = new Float_t[30]; - Float_t* etsigJet = new Float_t[30]; //signal et in jet - Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable) - Int_t* ncellsJet = new Int_t[30]; - Int_t* multJet = new Int_t[30]; + Float_t etaJet[kMaxJets]; + Float_t phiJet[kMaxJets]; + Float_t etJet[kMaxJets]; + Float_t etsigJet[kMaxJets]; //signal et in jet + Float_t etallJet[kMaxJets]; // total et in jet (tmp variable) + Int_t ncellsJet[kMaxJets]; + Int_t multJet[kMaxJets]; Int_t nJets; // to hold number of jets found by algorithm Int_t nj; // number of jets accepted Float_t prec = header->GetPrecBg(); Float_t bgprec = 1; while(bgprec > prec){ //reset jet arrays in memory - memset(etaJet,0,sizeof(Float_t)*30); - memset(phiJet,0,sizeof(Float_t)*30); - memset(etJet,0,sizeof(Float_t)*30); - memset(etallJet,0,sizeof(Float_t)*30); - memset(etsigJet,0,sizeof(Float_t)*30); - memset(ncellsJet,0,sizeof(Int_t)*30); - memset(multJet,0,sizeof(Int_t)*30); + memset(etaJet,0,sizeof(Float_t)*kMaxJets); + memset(phiJet,0,sizeof(Float_t)*kMaxJets); + memset(etJet,0,sizeof(Float_t)*kMaxJets); + memset(etallJet,0,sizeof(Float_t)*kMaxJets); + memset(etsigJet,0,sizeof(Float_t)*kMaxJets); + memset(ncellsJet,0,sizeof(Int_t)*kMaxJets); + memset(multJet,0,sizeof(Int_t)*kMaxJets); nJets = 0; nj = 0; // reset particles-jet array in memory @@ -155,9 +194,32 @@ void AliUA1JetFinderV1::FindJets() etbgTotal = etbgTotalN; // update with new background estimation } //end while + // add tracks to the jet if it wasn't yet done + if (header->GetBackgMode() == 0){ + Float_t rc= header->GetRadius(); + 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 + injet[jpart] = ijet; + break; + } + }// end jets loop + } //end particle loop + } + // add jets to list - Int_t* idxjets = new Int_t[nj]; + Int_t idxjets[kMaxJets]; Int_t nselectj = 0; + + TRefArray *refs = 0; + Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader"); + if (fromAod) refs = fReader->GetReferences(); + Float_t rc= header->GetRadius(); for(Int_t kj=0; kj (header->GetJetEtaMax())) || (etaJet[kj] < (header->GetJetEtaMin())) || @@ -167,14 +229,43 @@ void AliUA1JetFinderV1::FindJets() py = etJet[kj] * TMath::Sin(phiJet[kj]); pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj]))); en = TMath::Sqrt(px * px + py * py + pz * pz); - fJets->AddJet(px, py, pz, en); + + AliAODJet jet(px, py, pz, en); + + if (fromAod){ + for(Int_t jpart = 0; jpart < nIn; jpart++) // loop for all particles in array + if (injet[jpart] == kj && fReader->GetCutFlag(jpart) == 1) + jet.AddTrack(refs->At(jpart)); // check if the particle belongs to the jet and add the ref + } + + //jet.Print(""); + + // calculate the area of the jet + Float_t detamax = etaJet[kj] + rc; + Float_t detamin = etaJet[kj] - rc; + Float_t accmax = 0.0; Float_t accmin = 0.0; + if(detamax > header->GetLegoEtaMax()){ // sector outside etamax + Float_t h = header->GetLegoEtaMax() - etaJet[kj]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if(detamin < header->GetLegoEtaMin()){ // sector outside etamin + Float_t h = header->GetLegoEtaMax() + etaJet[kj]; + accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + Float_t areaJet = rc*rc*TMath::Pi() - accmax - accmin; + // set both areas + jet.SetEffArea(areaJet,areaJet); + + AddJet(jet); + idxjets[nselectj] = kj; nselectj++; - } + } //end particle loop + //add signal percentage and total signal in AliJets for analysis tool - Float_t* percentage = new Float_t[nselectj]; - Int_t* ncells = new Int_t[nselectj]; - Int_t* mult = new Int_t[nselectj]; + Float_t percentage[kMaxJets]; + Int_t ncells[kMaxJets]; + Int_t mult[kMaxJets]; for(Int_t i = 0; i< nselectj; i++){ percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]]; ncells[i] = ncellsJet[idxjets[i]]; @@ -193,35 +284,12 @@ void AliUA1JetFinderV1::FindJets() } if(bflag == 0) injet[bj] = -1; // set as background particle } - fJets->SetNCells(ncells); - fJets->SetPtFromSignal(percentage); - fJets->SetMultiplicities(mult); - fJets->SetInJet(injet); - fJets->SetEtaIn(etaT); - fJets->SetPhiIn(phiT); - fJets->SetPtIn(ptT); - fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi())); - //delete - delete ptT; - delete etaT; - delete phiT; - delete injet; - delete hPtTotal; - delete etaJet; - delete phiJet; - delete etJet; - delete etsigJet; - delete etallJet; - delete ncellsJet; - delete multJet; - delete idxjets; - delete percentage; - delete ncells; - delete mult; - - + delete [] ptT; + delete [] etaT; + delete [] phiT; + delete [] injet; } //////////////////////////////////////////////////////////////////////// @@ -232,19 +300,26 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& { //dump lego - // check enough space! *to be done* AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; - Float_t etCell[60000]; //! Cell Energy - Float_t etaCell[60000]; //! Cell eta - Float_t phiCell[60000]; //! Cell phi - Int_t flagCell[60000]; //! Cell flag + const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory + + const Int_t nBinEta = header->GetLegoNbinEta(); + const Int_t nBinPhi = header->GetLegoNbinPhi(); + if((nBinPhi*nBinEta)>nBinsMax){ + AliError("Too many bins of the ETA-PHI histogram"); + } + + Float_t etCell[nBinsMax]; //! Cell Energy + Float_t etaCell[nBinsMax]; //! Cell eta + Float_t phiCell[nBinsMax]; //! Cell phi + Short_t flagCell[nBinsMax]; //! Cell flag Int_t nCell = 0; TAxis* xaxis = fLego->GetXaxis(); TAxis* yaxis = fLego->GetYaxis(); Float_t e = 0.0; - for (Int_t i = 1; i <= header->GetLegoNbinEta(); i++) { - for (Int_t j = 1; j <= header->GetLegoNbinPhi(); j++) { + for (Int_t i = 1; i <= nBinEta; i++) { + for (Int_t j = 1; j <= nBinPhi; j++) { e = fLego->GetBinContent(i,j); if (e < 0.0) continue; // don't include this cells Float_t eta = xaxis->GetBinCenter(i); @@ -252,7 +327,7 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& etCell[nCell] = e; etaCell[nCell] = eta; phiCell[nCell] = phi; - flagCell[nCell] = 0; //default + flagCell[nCell] = 0; //default nCell++; } } @@ -267,33 +342,34 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& // tmp array of jets form algoritm - Float_t etaAlgoJet[30]; - Float_t phiAlgoJet[30]; - Float_t etAlgoJet[30]; - Int_t ncellsAlgoJet[30]; + Float_t etaAlgoJet[kMaxJets] = {0.0}; + Float_t phiAlgoJet[kMaxJets] = {0.0}; + Float_t etAlgoJet[kMaxJets] = {0.0}; + Int_t ncellsAlgoJet[kMaxJets] = {0}; //run algorithm// // sort cells by et - Int_t * index = new Int_t[nCell]; + Int_t index[nBinsMax]; TMath::Sort(nCell, etCell, index); // variable used in centroide loop - Float_t eta = 0.0; - Float_t phi = 0.0; - Float_t eta0 = 0.0; - Float_t phi0 = 0.0; - Float_t etab = 0.0; - Float_t phib = 0.0; - Float_t etas = 0.0; - Float_t phis = 0.0; - Float_t ets = 0.0; - Float_t deta = 0.0; - Float_t dphi = 0.0; - Float_t dr = 0.0; - Float_t etsb = 0.0; + Float_t eta = 0.0; + Float_t phi = 0.0; + Float_t eta0 = 0.0; + Float_t phi0 = 0.0; + Float_t etab = 0.0; + Float_t phib = 0.0; + Float_t etas = 0.0; + Float_t phis = 0.0; + Float_t ets = 0.0; + Float_t deta = 0.0; + Float_t dphi = 0.0; + Float_t dr = 0.0; + Float_t etsb = 0.0; Float_t etasb = 0.0; Float_t phisb = 0.0; - + Float_t dphib = 0.0; + for(Int_t icell = 0; icell < nCell; icell++){ Int_t jcell = index[icell]; @@ -313,21 +389,21 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& phisb = 0.0; for(Int_t kcell =0; kcell < nCell; kcell++){ Int_t lcell = index[kcell]; - if(lcell == jcell) continue; // cell itself - if(flagCell[lcell] != 0) continue; // cell used before - if(etCell[lcell] > etCell[jcell]) continue; + if(lcell == jcell) continue; // cell itself + if(flagCell[lcell] != 0) continue; // cell used before + if(etCell[lcell] > etCell[jcell]) continue; // can this happen //calculate dr deta = etaCell[lcell] - eta; - dphi = phiCell[lcell] - phi; - if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); - if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - dr = TMath::Sqrt(deta * deta + dphi * dphi); + dphi = TMath::Abs(phiCell[lcell] - phi); + if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi; + dr = TMath::Sqrt(deta * deta + dphi * dphi); if(dr <= rc){ // calculate offset from initiate cell deta = etaCell[lcell] - eta0; dphi = phiCell[lcell] - phi0; - if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); - if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + if (dphi < - TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi(); + etas = etas + etCell[lcell]*deta; phis = phis + etCell[lcell]*dphi; ets = ets + etCell[lcell]; @@ -335,24 +411,26 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& eta = eta0 + etas/ets; phi = phi0 + phis/ets; // if cone does not move much, just go to next step - dr = TMath::Sqrt((eta-etab)*(eta-etab) + (phi-phib)*(phi-phib)); + dphib = TMath::Abs(phi - phib); + if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib; + dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib); if(dr <= minmove) break; // cone should not move more than max_mov dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets)); if(dr > maxmove){ - eta = etab; - phi = phib; - ets = etsb; - etas = etasb; - phis = phisb; - }else{ // store this loop information - etab=eta; - phib=phi; - etsb = ets; - etasb = etas; - phisb = phis; + eta = etab; + phi = phib; + ets = etsb; + etas = etasb; + phis = phisb; + } else { // store this loop information + etab = eta; + phib = phi; + etsb = ets; + etasb = etas; + phisb = phis; } - } + } // inside cone }//end of cells loop looking centroide //avoid cones overloap (to be implemented in the future) @@ -394,19 +472,24 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& } //store tmp jet info !!! if(etbmax < etcmin) { - etaAlgoJet[nJets] = eta; - phiAlgoJet[nJets] = phi; - etAlgoJet[nJets] = etCone; - ncellsAlgoJet[nJets] = nCellIn; - nJets++; - } - + if(nJets %d) found by UA1JetFinder, adapt your cuts",kMaxJets)); + break; + } + } } // end of cells loop //reorder jets by et in cone //sort jets by energy - Int_t * idx = new Int_t[nJets]; - TMath::Sort(nJets, etAlgoJet, idx); + Int_t idx[kMaxJets]; + TMath::Sort(nJets, etAlgoJet, idx); // sort only the found jets for(Int_t p = 0; p < nJets; p++){ etaJet[p] = etaAlgoJet[idx[p]]; phiJet[p] = phiAlgoJet[idx[p]]; @@ -415,25 +498,20 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& ncellsJet[p] = ncellsAlgoJet[idx[p]]; } - - //delete - delete index; - delete idx; - } //////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV1::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, - Int_t* multJet, Int_t* injet) +void AliUA1JetFinderV1::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, + 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* multJet, Int_t* injet) { //background subtraction using cone method but without correction in dE/deta distribution //calculate energy inside and outside cones AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; Float_t rc= header->GetRadius(); - Float_t etIn[30]; + Float_t etIn[kMaxJets] = {0}; Float_t etOut = 0; for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut @@ -458,7 +536,7 @@ void AliUA1JetFinderV1::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, } //end particle loop //estimate jets and background areas - Float_t areaJet[30]; + Float_t areaJet[kMaxJets]; Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi(); for(Int_t k=0; kGetRadius(); - Float_t etIn[30]; + Float_t etIn[kMaxJets] = {0.0}; for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array //if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut @@ -525,7 +603,7 @@ void AliUA1JetFinderV1::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotal } //end particle loop //calc jets areas - Float_t areaJet[30]; + Float_t areaJet[kMaxJets]; Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi(); for(Int_t k=0; kReset(); + fhAreaJet[mjet]->Reset(); + } // background energy and area - TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); - TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); + 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 @@ -592,14 +678,14 @@ void AliUA1JetFinderV1::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTota injet[jpart] = ijet; multJet[ijet]++; if((fReader->GetCutFlag(jpart)) == 1){// pt cut - hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone + fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; } break; } }// end jets loop if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1) - hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones } //end particle loop //calc areas @@ -629,10 +715,10 @@ void AliUA1JetFinderV1::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTota if((eta0 < etaJet[ijet]) && (etaJet[ijet]Fill(etac,areaj); + fhAreaJet[ijet]->Fill(etac,areaj); areabg = areabg - areaj; } // end jets loop - hAreaBackg->Fill(etac,areabg); + fhAreaBackg->Fill(etac,areabg); eta0 = eta1; eta1 = eta1 + etaw; } // end loop for all eta bins @@ -641,38 +727,29 @@ void AliUA1JetFinderV1::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTota for(Int_t kjet=0; kjetGetBinContent(bin)){ - Float_t areab = hAreaBackg->GetBinContent(bin); - Float_t etb = hEtBackg->GetBinContent(bin); - Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; - etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction + if(fhAreaJet[kjet]->GetBinContent(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 } } } // calc background total - Double_t etOut = hEtBackg->Integral(); - Double_t areaOut = hAreaBackg->Integral(); + Double_t etOut = fhEtBackg->Integral(); + Double_t areaOut = fhAreaBackg->Integral(); Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); etbgTotalN = etOut*areaT/areaOut; - - //delete - for(Int_t ljet=0; ljetReset(); + fhAreaJet[mjet]->Reset(); + } // background energy and area - TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); // change range - TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); // change range + 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 @@ -712,13 +798,13 @@ void AliUA1JetFinderV1::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTot multJet[ijet]++; injet[jpart] = ijet; if((fReader->GetCutFlag(jpart)) == 1){ //pt cut - hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut + fhEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; } break; } }// end jets loop - if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + if(injet[jpart] == -1) fhEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones } //end particle loop //calc areas @@ -748,10 +834,10 @@ void AliUA1JetFinderV1::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTot if((eta0 < etaJet[ijet]) && (etaJet[ijet]Fill(etac,areaj); + fhAreaJet[ijet]->Fill(etac,areaj); areabg = areabg - areaj; } // end jets loop - hAreaBackg->Fill(etac,areabg); + fhAreaBackg->Fill(etac,areabg); eta0 = eta1; eta1 = eta1 + etaw; } // end loop for all eta bins @@ -760,29 +846,20 @@ void AliUA1JetFinderV1::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTot for(Int_t kjet=0; kjetGetBinContent(bin)){ - Float_t areab = hAreaBackg->GetBinContent(bin); - Float_t etb = hEtBackg->GetBinContent(bin); - Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; - etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction + if(fhAreaJet[kjet]->GetBinContent(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 } } } // calc background total - Double_t etOut = hEtBackg->Integral(); - Double_t areaOut = hAreaBackg->Integral(); + Double_t etOut = fhEtBackg->Integral(); + Double_t areaOut = fhAreaBackg->Integral(); Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); etbgTotalN = etOut*areaT/areaOut; - - //delete - for(Int_t ljet=0; ljetReset(); - fJets->ClearJets(); + AliJetFinder::Reset(); } //////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV1::WriteJHeaderToFile() +void AliUA1JetFinderV1::WriteJHeaderToFile() const { AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; - fOut->cd(); header->Write(); } @@ -809,11 +885,12 @@ void AliUA1JetFinderV1::Init() { // initializes some variables AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; - // book lego fLego = new TH2F("legoH","eta-phi", header->GetLegoNbinEta(), header->GetLegoEtaMin(), header->GetLegoEtaMax(), header->GetLegoNbinPhi(), header->GetLegoPhiMin(), header->GetLegoPhiMax()); + // Do not store in current dir + fLego->SetDirectory(0); }