X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliUA1JetFinderV1.cxx;h=d8f74f9a49cf3dd1d38cc47fbc7da1bf63f60cca;hb=1afa951f28e5a1d66cb0b78bb90e895caebe1579;hp=378c8c6d6d7a0502cb980bf38befad725732b139;hpb=952b368cd77fe683068c6d2bdcc0f4bec7827115;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliUA1JetFinderV1.cxx b/JETAN/AliUA1JetFinderV1.cxx index 378c8c6d6d7..d8f74f9a49c 100644 --- a/JETAN/AliUA1JetFinderV1.cxx +++ b/JETAN/AliUA1JetFinderV1.cxx @@ -46,9 +46,14 @@ ClassImp(AliUA1JetFinderV1) AliUA1JetFinderV1::AliUA1JetFinderV1() : AliJetFinder(), - fLego(0) + fLego(0), + fhEtBackg(0), + fhAreaBackg(0) { // Constructor + for(int i = 0;i < kMaxJets;i++){ + fhAreaJet[i] = fhEtJet[i] = 0; + } } //////////////////////////////////////////////////////////////////////// @@ -59,6 +64,16 @@ 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; + } + } //////////////////////////////////////////////////////////////////////// @@ -82,19 +97,29 @@ 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 + // 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(); @@ -102,37 +127,44 @@ void AliUA1JetFinderV1::FindJets() 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]); + npart += 1; etbgTotal+= ptT[i]; + etbg2 += ptT[i]*ptT[i]; } // 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 @@ -181,13 +213,13 @@ void AliUA1JetFinderV1::FindJets() } // add jets to list - Int_t* idxjets = new Int_t[nj]; + Int_t idxjets[kMaxJets]; Int_t nselectj = 0; -// printf("Found %d jets \n", nj); - + 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())) || @@ -208,6 +240,22 @@ void AliUA1JetFinderV1::FindJets() //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; @@ -215,9 +263,9 @@ void AliUA1JetFinderV1::FindJets() } //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]]; @@ -242,20 +290,6 @@ void AliUA1JetFinderV1::FindJets() 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; - - } //////////////////////////////////////////////////////////////////////// @@ -308,15 +342,15 @@ 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; @@ -438,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]]; @@ -459,11 +498,6 @@ void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& ncellsJet[p] = ncellsAlgoJet[idx[p]]; } - - //delete - delete[] index; - delete[] idx; - } //////////////////////////////////////////////////////////////////////// @@ -477,7 +511,7 @@ void AliUA1JetFinderV1::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t& //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 @@ -502,7 +536,7 @@ void AliUA1JetFinderV1::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t& } //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 @@ -569,7 +603,7 @@ void AliUA1JetFinderV1::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float } //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 @@ -636,14 +678,14 @@ void AliUA1JetFinderV1::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float 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 @@ -673,10 +715,10 @@ void AliUA1JetFinderV1::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float 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 @@ -685,29 +727,20 @@ void AliUA1JetFinderV1::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float 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 @@ -756,13 +798,13 @@ void AliUA1JetFinderV1::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ, Flo 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 @@ -792,10 +834,10 @@ void AliUA1JetFinderV1::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ, Flo 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 @@ -804,29 +846,20 @@ void AliUA1JetFinderV1::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ, Flo 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; ljetWrite();