X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliFastJetFinder.cxx;h=ddce4776bc806f0de163ca771109b53f887af6cd;hb=edf8b2d403af3ebab01765b0b41d62d0b8e470a0;hp=c61d90ad1d0354ed0df32ded1fae50fddfecbde9;hpb=a38beabda48835b1fda13db8c0b1cab0b0408820;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliFastJetFinder.cxx b/JETAN/AliFastJetFinder.cxx index c61d90ad1d0..ddce4776bc8 100644 --- a/JETAN/AliFastJetFinder.cxx +++ b/JETAN/AliFastJetFinder.cxx @@ -13,6 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ //--------------------------------------------------------------------- // FastJet v2.3.4 finder algorithm interface @@ -21,93 +22,74 @@ // Authors: Rafael.Diaz.Valdes@cern.ch // Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option) // +// ** 2011 magali.estienne@subatech.in2p3.fr & alexandre.shabetai@cern.ch +// new implementation of background subtraction +// allowing to subtract bkg using a different algo than the one used for signal jets //--------------------------------------------------------------------- - #include -#include -#include -#include -#include -#include -#include #include "AliFastJetFinder.h" #include "AliFastJetHeaderV1.h" -#include "AliJetReaderHeader.h" -#include "AliJetReader.h" -#include "AliJetUnitArray.h" #include "AliFastJetInput.h" -#include "AliJetBkg.h" +#include "AliFastJetBkg.h" #include "AliAODJetEventBackground.h" -#include "AliAODTrack.h" +#include "AliAODJet.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 using namespace std; - ClassImp(AliFastJetFinder) - -//____________________________________________________________________________ +//////////////////////////////////////////////////////////////////////// AliFastJetFinder::AliFastJetFinder(): AliJetFinder(), fInputFJ(new AliFastJetInput()), - fJetBkg(new AliJetBkg()) + fJetBkg(new AliFastJetBkg()) { // Constructor } //____________________________________________________________________________ - AliFastJetFinder::~AliFastJetFinder() { // destructor - delete fInputFJ; fInputFJ = 0; - delete fJetBkg; fJetBkg = 0; + delete fInputFJ; + delete fJetBkg; + } //______________________________________________________________________________ void AliFastJetFinder::FindJets() { - cout<<"----------in AliFastJetFinder::FindJets() ------------------"<GetDebug(); // debug option + Int_t debug = header->GetDebug(); // debug option Bool_t bgMode = header->GetBGMode(); // choose to subtract BG or not + if(debug>0) cout<<"----------in AliFastJetFinder::FindJets() ------------------"<ClassName(),"AliJetAODReader"); - if (fromAod) { refs = fReader->GetReferences(); } - // RUN ALGORITHM // read input particles ----------------------------- vector inputParticles=fInputFJ->GetInputParticles(); if(inputParticles.size()==0){ - Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); + if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); return; } // create an object that represents your choice of jet algorithm, and // the associated parameters double rParam = header->GetRparam(); + double rBkgParam = header->GetRparamBkg(); fastjet::Strategy strategy = header->GetStrategy(); fastjet::RecombinationScheme recombScheme = header->GetRecombScheme(); fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); @@ -115,8 +97,8 @@ void AliFastJetFinder::FindJets() // create an object that specifies how we to define the area fastjet::AreaDefinition areaDef; - double ghostEtamax = header->GetGhostEtaMax(); - double ghostArea = header->GetGhostArea(); + double ghostEtamax = header->GetGhostEtaMax(); + double ghostArea = header->GetGhostArea(); int activeAreaRepeats = header->GetActiveAreaRepeats(); // now create the object that holds info about ghosts @@ -125,40 +107,40 @@ void AliFastJetFinder::FindJets() fastjet::AreaType areaType = header->GetAreaType(); areaDef = fastjet::AreaDefinition(areaType,ghostSpec); + //***************************** JETS FINDING + // run the jet clustering with the above jet definition + fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); + vector jets; - - - if(bgMode) // BG subtraction!!!!!!!!!! Not entering here + if(bgMode) // Do BG subtraction directly with the same algorithm (cambridge or kt) for jet signal and background { - //***************************** JETS FINDING AND EXTRACTION + //***************************** JETS FINDING FOR RHO ESTIMATION // run the jet clustering with the above jet definition - fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); + fastjet::JetAlgorithm algorithmBkg = header->GetBGAlgorithm(); + fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy); + fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef); // save a comment in the header - TString comment = "Running FastJet algorithm with the following setup. "; comment+= "Jet definition: "; comment+= TString(jetDef.description()); + comment+= "Jet bckg definition: "; + comment+= TString(jetDefBkg.description()); comment+= ". Area definition: "; comment+= TString(areaDef.description()); comment+= ". Strategy adopted by FastJet: "; comment+= TString(clust_seq.strategy_string()); header->SetComment(comment); - if(debug){ + if(debug>0){ cout << "--------------------------------------------------------" << endl; cout << comment << endl; cout << "--------------------------------------------------------" << endl; } - //header->PrintParameters(); - - // extract the inclusive jets with pt > ptmin, sorted by pt + // extract the inclusive jets sorted by pt double ptmin = header->GetPtMin(); - vector inclusiveJets = clust_seq.inclusive_jets(ptmin); - - //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; - + vector inclusiveJets = clust_seq.inclusive_jets(); //subtract background // =========================================== // set the rapididty , phi range within which to study the background @@ -167,132 +149,98 @@ void AliFastJetFinder::FindJets() double phiMax = header->GetPhiMax(); double phiMin = header->GetPhiMin(); fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); - - // subtract background - vector subJets = clust_seq.subtracted_jets(range,ptmin); - - // print out - //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n"; - //cout << "---------------------------------------\n"; - //cout << endl; - //printf(" ijet rap phi Pt area +- err\n"); + + // Extract rho and sigma + Double_t rho = 0.; + Double_t sigma = 0.; + Double_t meanarea = 0.; + Bool_t kUse4VectorArea = header->Use4VectorArea(); + vector bkgJets = clust_seq_bkg.inclusive_jets(); + clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, kUse4VectorArea, rho, sigma, meanarea, false); + + // subtract background and extract jets bkg subtracted + vector subJets = clust_seq.subtracted_jets(rho,ptmin); // sort jets into increasing pt - vector jets = sorted_by_pt(subJets); - for (size_t j = 0; j < jets.size(); j++) { // loop for jets - - double area = clust_seq.area(jets[j]); - double areaError = clust_seq.area_error(jets[j]); - - if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError); - - // go to write AOD info - AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); - //cout << "Printing jet " << endl; - if(debug) aodjet.Print(""); - //cout << "Adding jet ... " ; - AddJet(aodjet); - //cout << "added \n" << endl; - - } - } - - - - + jets = sorted_by_pt(subJets); + } else { // No BG subtraction!!!!!!!! Default header is bgmode=0. - - TClonesArray* fUnit = fReader->GetUnitArray(); //Big mmentum array - if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; } - Int_t nIn = fUnit->GetEntries(); - - - //fastjet::ClusterSequence clust_seq(inputParticles, jetDef); - fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); - - // save a comment in the header + // save a comment in the header TString comment = "Running FastJet algorithm with the following setup. "; comment+= "Jet definition: "; comment+= TString(jetDef.description()); comment+= ". Strategy adopted by FastJet: "; comment+= TString(clust_seq.strategy_string()); header->SetComment(comment); - if(debug){ + if(debug>0){ cout << "--------------------------------------------------------" << endl; cout << comment << endl; cout << "--------------------------------------------------------" << endl; } - //header->PrintParameters(); - // extract the inclusive jets with pt > ptmin, sorted by pt + // extract the inclusive jets with pt > ptmin, sorted by pt double ptmin = header->GetPtMin(); vector inclusiveJets = clust_seq.inclusive_jets(ptmin); - //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; - - vector jets = sorted_by_pt(inclusiveJets); // Added by me - for (size_t j = 0; j < jets.size(); j++) { // loop for jets // Added by me + jets = sorted_by_pt(inclusiveJets); + + } + + for (size_t j = 0; j < jets.size(); j++) { // loop for jets + + double area = clust_seq.area(jets[j]); + double areaError = clust_seq.area_error(jets[j]); + + if(debug>0) printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp()); + + vector constituents = clust_seq.constituents(jets[j]); + int nCon= constituents.size(); + TArrayI ind(nCon); - if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp()); - - vector constituents = clust_seq.constituents(jets[j]); - int nCon= constituents.size(); - TArrayI ind(nCon); - Double_t area=clust_seq.area(jets[j]); - // go to write AOD info - AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); - aodjet.SetEffArea(area,0); - //cout << "Printing jet " << endl; - if(debug) aodjet.Print(""); - Int_t count=0; - for (int i=0; i < nCon; i++) - { + if ((jets[j].eta() > (header->GetJetEtaMax())) || + (jets[j].eta() < (header->GetJetEtaMin())) || + (jets[j].phi() > (header->GetJetPhiMax())) || + (jets[j].phi() < (header->GetJetPhiMin())) || + (jets[j].perp() < header->GetPtMin())) continue; // acceptance eta range and etmin + + // go to write AOD info + AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); + aodjet.SetEffArea(area,areaError); + //cout << "Printing jet " << endl; + if(debug>0) aodjet.Print(""); + + for (int i=0; i < nCon; i++) + { fastjet::PseudoJet mPart=constituents[i]; ind[i]=mPart.user_index(); - // cout<> px >> py >> pz >> en; + if (!in.good()) break; + //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz); + nlines++; + inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en)); + } + //printf(" found %d pointsn",nlines); + in.close(); + ////////////////////////////////////////////////// // create an object that represents your choice of jet algorithm, and // the associated parameters @@ -318,32 +266,26 @@ void AliFastJetFinder::RunTest(const char* datafile) fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme; fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy); - - // create an object that specifies how we to define the area fastjet::AreaDefinition areaDef; double ghostEtamax = 7.0; double ghostArea = 0.05; int activeAreaRepeats = 1; - // now create the object that holds info about ghosts fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea); // and from that get an area definition areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec); - // run the jet clustering with the above jet definition fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); - // tell the user what was done cout << "--------------------------------------------------------" << endl; cout << "Jet definition was: " << jetDef.description() << endl; cout << "Area definition was: " << areaDef.description() << endl; cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()< 5 GeV, sorted by pt double ptmin = 5.0; @@ -351,7 +293,6 @@ void AliFastJetFinder::RunTest(const char* datafile) cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; - //subtract background // =========================================== // set the rapididty range within which to study the background double rapMax = ghostEtamax - rParam; @@ -367,143 +308,62 @@ void AliFastJetFinder::RunTest(const char* datafile) printf(" ijet rap phi Pt area +- err\n"); for (size_t j = 0; j < jets.size(); j++) { - - double area = clust_seq.area(jets[j]); + double area = clust_seq.area(jets[j]); double areaError = clust_seq.area_error(jets[j]); - printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(), jets[j].phi(),jets[j].perp(), area, areaError); } cout << endl; // ================================================================ - - } //____________________________________________________________________________ - -void AliFastJetFinder::WriteJHeaderToFile() +void AliFastJetFinder::WriteJHeaderToFile() const { + // Write Jet Header to file fHeader->Write(); -} - -//____________________________________________________________________________ - -Float_t AliFastJetFinder::EtaToTheta(Float_t arg) -{ - // return (180./TMath::Pi())*2.*atan(exp(-arg)); - return 2.*atan(exp(-arg)); - } //____________________________________________________________________________ - -void AliFastJetFinder::InitTask(TChain *tree) -{ - - printf("Fast jet finder initialization ******************"); - fReader->CreateTasks(tree); - -} - - Bool_t AliFastJetFinder::ProcessEvent() { - // - // Process one event - // from meomntum array - - Bool_t ok = fReader->FillMomentumArray(); - - if (!ok) return kFALSE; - fInputFJ->SetHeader(fHeader); - fInputFJ->SetReader(fReader); - fInputFJ->FillInput(); - // Jets - FindJets(); - - fJetBkg->SetHeader(fHeader); - fJetBkg->SetReader(fReader); - /* - fJetBkg->SetFastJetInput(fInputFJ); - Double_t bkg1=fJetBkg->BkgFastJet(); - Double_t bkg2=fJetBkg->BkgChargedFastJet(); - Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets); - Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets); - - fAODEvBkg->SetBackground(0,bkg1); - fAODEvBkg->SetBackground(1,bkg2); - fAODEvBkg->SetBackground(2,bkg3); - fAODEvBkg->SetBackground(3,bkg4); - */ - Reset(); - return kTRUE; - -} - -Bool_t AliFastJetFinder::ProcessEvent2() -{ - // // Process one event // Charged only or charged+neutral jets - // - - TRefArray* ref = new TRefArray(); - Bool_t procid = kFALSE; - Bool_t ok = fReader->ExecTasks(procid,ref); - // Delete reference pointer - if (!ok) {delete ref; return kFALSE;} - - // Leading particles fInputFJ->SetHeader(fHeader); - fInputFJ->SetReader(fReader); + fInputFJ->SetCalTrkEvent(GetCalTrkEvent()); fInputFJ->FillInput(); - // Jets - FindJets(); - - fJetBkg->SetHeader(fHeader); - fJetBkg->SetReader(fReader); - fJetBkg->SetFastJetInput(fInputFJ); - Double_t bkg1=fJetBkg->BkgFastJet(); - Double_t bkg2=fJetBkg->BkgChargedFastJet(); - Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets); - Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets); - - fAODEvBkg->SetBackground(0,bkg1); - fAODEvBkg->SetBackground(1,bkg2); - fAODEvBkg->SetBackground(2,bkg3); - fAODEvBkg->SetBackground(3,bkg4); - - Int_t nEntRef = ref->GetEntries(); - - for(Int_t i=0; iAt(i))->SetUnitTrackID(0); - ((AliJetUnitArray*)ref->At(i))->SetUnitEnergy(0.); - ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller); - ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller); - ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad); - ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlagC(kTRUE,kBad); - ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc); - ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet); - ((AliJetUnitArray*)ref->At(i))->ClearUnitTrackRef(); - - // Reset process ID - AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i); - uA->ResetBit(kIsReferenced); - uA->SetUniqueID(0); - } + // Find Jets + FindJets(); - // Delete the reference pointer - ref->Delete(); - delete ref; + // Background + if(fAODEvBkg){ + AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; - Reset(); + fJetBkg->SetHeader(fHeader); + fJetBkg->SetFastJetInput(fInputFJ); + + Int_t count = 0; + if(header->GetBkgFastJetb()){ + Double_t sigma1 = 0 , meanarea1= 0, bkg1 = 0; + fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1); + fAODEvBkg->SetBackground(count,bkg1,sigma1,meanarea1); + count++; + } + + if(header->GetBkgFastJetWoHardest()){ + Double_t sigma2 = 0 , meanarea2 = 0, bkg2 = 0; + fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2); + fAODEvBkg->SetBackground(count,bkg2,sigma2,meanarea2); + count++; + } + + } + Reset(); return kTRUE; + }