X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliSISConeJetFinder.cxx;h=07e58db1ef3d4c8c6d5d628e13a882d59da2bd03;hb=e00c3b87295b86443585849c524cbfbffe242cb0;hp=771138a9983cf2cf65b1d45312b53eaa7e4aa155;hpb=9b72221bcd40f171d1eb49295e32750386854022;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliSISConeJetFinder.cxx b/JETAN/AliSISConeJetFinder.cxx index 771138a9983..07e58db1ef3 100644 --- a/JETAN/AliSISConeJetFinder.cxx +++ b/JETAN/AliSISConeJetFinder.cxx @@ -13,32 +13,23 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ //--------------------------------------------------------------------- // FastJet v2.3.4 finder algorithm interface // // Author: swensy.jangal@ires.in2p3.fr // -// Last modification: Neutral cell energy included in the jet reconstruction -// -// Author: Magali.estienne@subatech.in2p3.fr +// ** 2011 magali.estienne@subatech.in2p3.fr & alexandre.shabetai@cern.ch +// Modified accordingly to reader/finder splitting and new handling of neutral information (via FastJetInput) //--------------------------------------------------------------------- - #include -#include -#include -#include -#include -#include -#include -#include - -#include "AliHeader.h" -#include "AliJetKineReader.h" -#include "AliJetReader.h" -#include "AliJetReaderHeader.h" -#include "AliJetUnitArray.h" + +#include "AliFastJetInput.h" +#include "AliFastJetBkg.h" +#include "AliAODJetEventBackground.h" +#include "AliAODJet.h" #include "AliSISConeJetFinder.h" #include "AliSISConeJetHeader.h" @@ -46,145 +37,78 @@ #include "fastjet/ClusterSequenceArea.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" - // get info on how fastjet was configured #include "fastjet/config.h" -#ifdef ENABLE_PLUGIN_SISCONE +#if defined(ENABLE_PLUGIN_SISCONE) || defined(FASTJET_ENABLE_PLUGIN_SISCONE) #include "fastjet/SISConePlugin.hh" #endif -#include // needed for internal io #include -#include -using namespace std; +using namespace std; ClassImp(AliSISConeJetFinder) - -//____________________________________________________________________________ +//////////////////////////////////////////////////////////////////////// AliSISConeJetFinder::AliSISConeJetFinder(): -AliJetFinder() + AliJetFinder(), + fInputFJ(new AliFastJetInput()), + fJetBkg(new AliFastJetBkg()) { // Constructor } //____________________________________________________________________________ - AliSISConeJetFinder::~AliSISConeJetFinder() { // destructor + delete fInputFJ; + delete fJetBkg; + } //______________________________________________________________________________ void AliSISConeJetFinder::FindJets() { + // run the SISCone Jet finder - // Pick up siscone header + // Pick up siscone header AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader; - Bool_t debug = header->GetDebug(); // debug option - Int_t fOpt = fReader->GetReaderHeader()->GetDetector(); - - // Check if we are reading AOD jets - TRefArray *refs = 0; - Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader"); - if (fromAod) { refs = fReader->GetReferences(); } - + Int_t debug = header->GetDebug(); // debug option + Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not + + // Read input particles + vector inputParticles=fInputFJ->GetInputParticles(); + if(inputParticles.size()==0){ + if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); + return; + } -//******************************** SISCONE PLUGIN CONFIGURATION -// Here we look for SISCone parameters in the header and we define our plugin. + //------------------- SISCONE PLUGIN CONFIGURATION ---------------------- + // Look for SISCone parameters in the header and definition of the plugin. - Double_t coneRadius = header->GetConeRadius(); // cone radius - Double_t overlapThreshold = header->GetOverlapThreshold(); // overlap parameter - Int_t nPassMax = header->GetNPassMax(); // maximum number of passes - Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets - Double_t caching = header->GetCaching(); // do we record found cones for this set of data? + Double_t coneRadius = header->GetConeRadius(); // cone radius + Double_t overlapThreshold = header->GetOverlapThreshold(); // overlap parameter + Int_t nPassMax = header->GetNPassMax(); // maximum number of passes + Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets + Double_t caching = header->GetCaching(); // do we record found cones for this set of data? + // For bckg + double rBkgParam = header->GetRparamBkg(); + fastjet::Strategy strategy = header->GetStrategy(); + fastjet::RecombinationScheme recombScheme = header->GetRecombScheme(); -// if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale -// Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets + // if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale + // Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets fastjet::JetDefinition::Plugin * plugin; plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching); - vector inputParticles; -//******************************** READING OF INPUT PARTICLES -// Here we look for px, py pz and energy of each particle that we gather in a PseudoJet object, and we put all these PseudoJet in a vector of PseudoJets : input_particles. - - if(fOpt==0) - { - TClonesArray *lvArray = fReader->GetMomentumArray(); - Int_t nIn = lvArray->GetEntries(); - - // We check if lvArray is ok - if(lvArray == 0) - { - cout << "Could not get the momentum array" << endl; - return; - } - - if(nIn == 0)// nIn = Number of particles in the event - { - if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; - return; - } - - Float_t px,py,pz,en; - - // Load input vectors - for(Int_t i = 0; i < nIn; i++) - { - TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); - px = lv->Px(); - py = lv->Py(); - pz = lv->Pz(); - en = lv->Energy(); - - fastjet::PseudoJet inputPart(px,py,pz,en); - inputPart.set_user_index(i); - inputParticles.push_back(inputPart); - - } - } - else { - TClonesArray* fUnit = fReader->GetUnitArray(); - if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; } - Int_t nIn = fUnit->GetEntries(); - if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } - // Information extracted from fUnitArray - // load input vectors and calculate total energy in array - Float_t pt,eta,phi,theta,px,py,pz,en; - Int_t ipart = 0; - for(Int_t i=0; iAt(i); - - if(uArray->GetUnitEnergy()>0.){ - - // It is not necessary anymore to cut on particle pt - pt = uArray->GetUnitEnergy(); - eta = uArray->GetUnitEta(); - phi = uArray->GetUnitPhi(); - theta = EtaToTheta(eta); - en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta)); - px = TMath::Cos(phi)*pt; - py = TMath::Sin(phi)*pt; - pz = en*TMath::TanH(eta); - if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl; - - fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object - inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm - inputParticles.push_back(inputPart); // back of the input_particles vector - ipart++; - } - } // End loop on UnitArray - } - -//******************************** CHOICE OF JET AREA -// Here we determine jets area for subtracting background later -// For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez + //------------------- CHOICE OF JET AREA ---------------------- + // Definition of jet areas for background subtraction + // For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated Double_t ghostArea = header->GetGhostArea(); // area of a ghost @@ -204,150 +128,187 @@ void AliSISConeJetFinder::FindJets() fastjet::AreaDefinition areaDef; if (areaTypeNumber < 5) - { - fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt); - areaDef = fastjet::AreaDefinition(areaType,ghostSpec); - } + { + fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt); + areaDef = fastjet::AreaDefinition(areaType,ghostSpec); + } if (areaTypeNumber == 5) - { - Double_t effectiveRFact = header->GetEffectiveRFact(); - fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact); - areaDef = fastjet::AreaDefinition(areaType,ghostSpec); - } - -//******************************** -//******************************** - - Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not + { + Double_t effectiveRFact = header->GetEffectiveRFact(); + fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact); + areaDef = fastjet::AreaDefinition(areaType,ghostSpec); + } -//******************************** -//******************************** + //------------------- JETS FINDING AND EXTRACTION ---------------------- + fastjet::ClusterSequenceArea clust_seq(inputParticles, plugin, areaDef); -//******************************** -//******************************** BG SUBTRACTION - if (bgMode == 1)// BG subtraction - { - //******************************** JETS FINDING AND EXTRACTION - fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef); - // Here we extract inclusive jets with pt > ptmin, sorted by pt - Double_t ptMin = header->GetMinJetPt(); - vector inclusiveJets = clustSeq.inclusive_jets(ptMin); - vector jets = sorted_by_pt(inclusiveJets); - - //***************************** BACKGROUND SUBTRACTION - - // Set the rapidity-azimuth range within which to study background - Double_t rapMin = header->GetRapMin(); - Double_t rapMax = header->GetRapMax(); - Double_t phiMin = header->GetPhiMin(); - Double_t phiMax = header->GetPhiMax(); - fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); - - // rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas) - Int_t algo = header->GetBGAlgorithm(); - fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm; - if (algo == 0) algorithm = fastjet::kt_algorithm; - if (algo == 1) algorithm = fastjet::cambridge_algorithm; - fastjet::JetDefinition jetDefForRho(algorithm, 0.5); - fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef); - Double_t rho = csForRho.median_pt_per_unit_area_4vector(range); - cout<<"rho = "< corrJets = clustSeq.subtracted_jets(rho, ptMin); - - //***************************** JETS DISPLAY + vector jets; - for (size_t j = 0; j < jets.size(); j++) + if (bgMode == 1)// BG subtraction mode { - // If the jet is only made of ghosts, continue. - if (clustSeq.is_pure_ghost(jets[j]) == 1) continue; - - // If the correction is > jet energy px = py = pz = e = 0 - if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue; - - if(debug){ - cout<<"********************************** Reconstructed jet(s) (non corrected)"<SetComment(comment); + if(debug>0){ + cout << "--------------------------------------------------------" << endl; + cout << comment << endl; + cout << "--------------------------------------------------------" << endl; } - // Go to write AOD info - AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E()); - if(debug) aodjet.Print(""); - AddJet(aodjet); - } - } -//******************************** -//******************************** NO BG SUBTRACTION - - if (bgMode == 0)// No BG subtraction - { - //******************************** JETS FINDING AND EXTRACTION - fastjet::ClusterSequence clustSeq(inputParticles, plugin); - // Here we extract inclusive jets with pt > ptmin, sorted by pt - Double_t ptMin = header->GetMinJetPt(); - vector inclusiveJets = clustSeq.inclusive_jets(ptMin); - vector jets = sorted_by_pt(inclusiveJets); - - //***************************** JETS DISPLAY + // Here we extract inclusive jets with pt > ptmin, sorted by pt + Double_t ptMin = header->GetMinJetPt(); + vector inclusiveJets = clust_seq.inclusive_jets(); // ptMin removed + // vector jets = sorted_by_pt(inclusiveJets); - for (size_t k = 0; k < jets.size(); k++) - { - if(debug){ - cout<<"********************************** Reconstructed jet(s) (non corrected)"<GetRapMin(); + Double_t rapMax = header->GetRapMax(); + Double_t phiMin = header->GetPhiMin(); + Double_t phiMax = header->GetPhiMax(); + fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); + // 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 + jets = sorted_by_pt(subJets); -void AliSISConeJetFinder::InitTask(TChain *tree) -{ + } + else // No BG subtraction + { - printf("SISCone jet finder initialization ******************"); - fReader->CreateTasks(tree); + // save a comment in the header + TString comment = "Running Siscone algorithm with the following setup. "; + comment+= "Jet definition: "; + // comment+= TString(plugin.description()); + comment+= ". Strategy adopted by FastJet: "; + comment+= TString(clust_seq.strategy_string()); + header->SetComment(comment); + if(debug>0){ + cout << "--------------------------------------------------------" << endl; + cout << comment << endl; + cout << "--------------------------------------------------------" << endl; + } + //header->PrintParameters(); -} + // Here we extract inclusive jets with pt > ptmin, sorted by pt + Double_t ptMin = header->GetMinJetPt(); + vector inclusiveJets = clust_seq.inclusive_jets(ptMin); + jets = sorted_by_pt(inclusiveJets); + } -//____________________________________________________________________________ + //------------------- JET AND TRACK STORAGE ---------------------- + 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 %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError); + + vector constituents = clust_seq.constituents(jets[j]); + int nCon= constituents.size(); + TArrayI ind(nCon); + + 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->GetMinJetPt())) 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(); + + // Jet constituents (charged tracks) added to the AliAODJet + AliJetCalTrkEvent* calEvt = GetCalTrkEvent(); + for(Int_t itrack=0; itrackGetNCalTrkTracks(); itrack++) + { + if(itrack==ind[i]) + { + TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject(); + aodjet.AddTrack(track); + } + } + } // End loop on Constituents + + AddJet(aodjet); + + } // end loop for jets + + delete plugin; -void AliSISConeJetFinder::WriteJHeaderToFile() +} + +//____________________________________________________________________________ +void AliSISConeJetFinder::WriteJHeaderToFile() const { + // Write Jet Header To File( fHeader->Write(); } //____________________________________________________________________________ +Bool_t AliSISConeJetFinder::ProcessEvent() +{ + // Process one event + // Charged only or charged+neutral jets + + fInputFJ->SetHeader(fHeader); + fInputFJ->SetCalTrkEvent(GetCalTrkEvent()); + fInputFJ->FillInput(); + + // Jets + FindJets(); + + // Background + if( fAODEvBkg){ + fJetBkg->SetHeader(fHeader); + Double_t sigma1 = 0,meanarea1= 0,sigma2 = 0,meanarea2 = 0; + Double_t bkg1 = 0,bkg2 = 0; + + fJetBkg->SetFastJetInput(fInputFJ); + fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1); + fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2); + fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1); + fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2); + } + Reset(); + return kTRUE; + +}