X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliSISConeJetFinder.cxx;h=f9ef55f8b91a461755e8f2f5b9edb23519d3ff94;hb=cb0a52e080e53665e8d4500dee3ce41be919a8a6;hp=307a2b77300666c23d86fe9f8ce81dd862b86e6e;hpb=2551af9d9203c2cd9471980cdbc1205d7bad018d;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliSISConeJetFinder.cxx b/JETAN/AliSISConeJetFinder.cxx index 307a2b77300..f9ef55f8b91 100644 --- a/JETAN/AliSISConeJetFinder.cxx +++ b/JETAN/AliSISConeJetFinder.cxx @@ -15,26 +15,30 @@ //--------------------------------------------------------------------- -// SISCone (FastJet v2.3.4) finder algorithm interface +// 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 //--------------------------------------------------------------------- + #include #include -#include #include #include #include #include #include +#include #include "AliHeader.h" -#include "AliJet.h" #include "AliJetKineReader.h" #include "AliJetReader.h" #include "AliJetReaderHeader.h" +#include "AliJetUnitArray.h" #include "AliSISConeJetFinder.h" #include "AliSISConeJetHeader.h" @@ -43,21 +47,23 @@ #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" -// Get info on how fastjet was configured +// 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 // needed for internal io #include #include using namespace std; + ClassImp(AliSISConeJetFinder) + //____________________________________________________________________________ AliSISConeJetFinder::AliSISConeJetFinder(): @@ -77,16 +83,11 @@ AliSISConeJetFinder::~AliSISConeJetFinder() void AliSISConeJetFinder::FindJets() { - Bool_t debug = kFALSE; - // Pick up siscone header AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader; + Int_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(); } - //******************************** SISCONE PLUGIN CONFIGURATION // Here we look for SISCone parameters in the header and we define our plugin. @@ -97,83 +98,254 @@ void AliSISConeJetFinder::FindJets() Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets Double_t caching = header->GetCaching(); // do we record found cones for this set of data? - if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale +// 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. - TClonesArray *lvArray = fReader->GetMomentumArray(); - Int_t nIn = lvArray->GetEntries(); - - // We check if lvArray is ok - if(lvArray == 0) + if(fOpt==0) + { + TClonesArray *lvArray = fReader->GetMomentumArray(); + + + // We check if lvArray is ok + if(lvArray == 0) + { + cout << "Could not get the momentum array" << endl; + delete plugin; + return; + } + + Int_t nIn = lvArray->GetEntries(); + + if(nIn == 0)// nIn = Number of particles in the event + { + if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; + delete plugin; + 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; delete plugin; return; } + Int_t nIn = fUnit->GetEntries(); + if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; delete plugin; 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 + + Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated + Double_t ghostArea = header->GetGhostArea(); // area of a ghost + Int_t activeAreaRepeats = header->GetActiveAreaRepeats(); // do we repeat area calculation? + Double_t gridScatter = header->GetGridScatter(); // fractional random fluctuations of the position of the ghosts on the y-phi grid + Double_t ktScatter = header->GetKtScatter(); // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid + Double_t meanGhostKt = header->GetMeanGhostKt(); // average transverse momentum of the ghosts. + + Double_t areaTypeNumber = header->GetAreaTypeNumber(); // the number determines jet area type + fastjet::AreaType areaType = fastjet::active_area; + if (areaTypeNumber == 1) areaType = fastjet::active_area; + if (areaTypeNumber == 2) areaType = fastjet::active_area_explicit_ghosts; + if (areaTypeNumber == 3) areaType = fastjet::one_ghost_passive_area; + if (areaTypeNumber == 4) areaType = fastjet::passive_area; + if (areaTypeNumber == 5) areaType = fastjet::voronoi_area; + + fastjet::AreaDefinition areaDef; + + if (areaTypeNumber < 5) { - cout << "Could not get the momentum array" << endl; - return; + fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt); + areaDef = fastjet::AreaDefinition(areaType,ghostSpec); } - if(nIn == 0)// nIn = Number of particles in the event + if (areaTypeNumber == 5) { - if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; - return; + Double_t effectiveRFact = header->GetEffectiveRFact(); + fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact); + areaDef = fastjet::AreaDefinition(areaType,ghostSpec); } - Int_t nJets = 0; // Number of jets in this event - fJets->SetNinput(nIn) ; // fJets = AliJet number of input objects - Float_t px,py,pz,en; - vector inputParticles; +//******************************** +//******************************** - // 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); - - } - -//******************************** JETS FINDING + Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not - fastjet::ClusterSequence clustSeq(inputParticles, plugin); +//******************************** +//******************************** -//***************************** JETS EXTRACTION AND CORRECTION +//******************************** +//******************************** 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; + Double_t RRho = header->GetRForRho(); + fastjet::JetDefinition jetDefForRho(algorithm, RRho); + 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 + + for (size_t j = 0; j < jets.size(); j++) + { + // 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)"< ptmin, sorted by pt - Double_t ptMin = header->GetMinJetPt(); - vector inclusiveJets = clustSeq.inclusive_jets(ptMin); - vector jets = sorted_by_pt(inclusiveJets); +//******************************** +//******************************** NO BG SUBTRACTION - for (Int_t j = 0; j < jets.size(); j++) + if (bgMode == 0)// No BG subtraction { - cout<<"********************************** Reconstructed jet(s) (non corrected)"< ptmin, sorted by pt + Double_t ptMin = header->GetMinJetPt(); + vector inclusiveJets = clustSeq.inclusive_jets(ptMin); + vector jets = sorted_by_pt(inclusiveJets); + + //***************************** JETS DISPLAY + + for (size_t k = 0; k < jets.size(); k++) + { + if(debug) + { + cout<<"********************************** Reconstructed jet(s) (non corrected)"<CreateTasks(tree); + +} + + +//____________________________________________________________________________ + +void AliSISConeJetFinder::WriteJHeaderToFile() const { fHeader->Write(); } + +//____________________________________________________________________________ +