X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliFastJetFinder.cxx;h=f15989bffd59ea8ea7adb715a2276608bcf8ac6b;hb=f754f2c452d0ee08a9a0aaef1933655e6f407604;hp=68ef13aefbdbccdf9e65acf1123ac7ac81800a09;hpb=6cbd6b57710e627498828d98eca7a726febd135a;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliFastJetFinder.cxx b/JETAN/AliFastJetFinder.cxx index 68ef13aefbd..f15989bffd5 100644 --- a/JETAN/AliFastJetFinder.cxx +++ b/JETAN/AliFastJetFinder.cxx @@ -96,7 +96,7 @@ void AliFastJetFinder::FindJets() TRefArray *refs = 0; Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader"); if (fromAod) { refs = fReader->GetReferences(); } - + // RUN ALGORITHM // read input particles ----------------------------- @@ -109,6 +109,7 @@ void AliFastJetFinder::FindJets() // 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(); @@ -126,41 +127,42 @@ 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 + //***************************** CLUSTER 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+= ". Strategy adopted by FastJet and bkg: "; comment+= TString(clust_seq.strategy_string()); header->SetComment(comment); if(debug){ - cout << "--------------------------------------------------------" << endl; - cout << comment << endl; - cout << "--------------------------------------------------------" << endl; + cout << "--------------------------------------------------------" << endl; + cout << comment << endl; + cout << "--------------------------------------------------------" << endl; } //header->PrintParameters(); - // 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 inclusiveJets = clust_seq.inclusive_jets(); + //subtract background // =========================================== // set the rapididty , phi range within which to study the background double rapMax = header->GetRapMax(); @@ -169,121 +171,128 @@ void AliFastJetFinder::FindJets() double phiMin = header->GetPhiMin(); fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); - // subtract background - vector subJets = clust_seq.subtracted_jets(range,ptmin); - - // print out + // Extract rho and sigma + Double_t rho = 0.; + Double_t sigma = 0.; + Double_t meanarea = 0.; + Bool_t Use4VectorArea = header->Use4VectorArea(); + vector bkgJets = clust_seq_bkg.inclusive_jets(); + clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, Use4VectorArea, rho, sigma, meanarea, false); + + // subtract background and extract jets bkg subtracted + vector subJets = clust_seq.subtracted_jets(rho,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"); // 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. + // 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){ + cout << "--------------------------------------------------------" << endl; + cout << comment << endl; + cout << "--------------------------------------------------------" << endl; + } + //header->PrintParameters(); + // extract the inclusive jets with pt > ptmin, sorted by pt + double ptmin = header->GetPtMin(); + vector inclusiveJets = clust_seq.inclusive_jets(ptmin); + jets = sorted_by_pt(inclusiveJets); // Added by me + } + for (size_t j = 0; j < jets.size(); j++) { // loop for jets - 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(); + double area = clust_seq.area(jets[j]); + double areaError = clust_seq.area_error(jets[j]); - - //fastjet::ClusterSequence clust_seq(inputParticles, jetDef); - fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); - - // 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){ - cout << "--------------------------------------------------------" << endl; - cout << comment << endl; - cout << "--------------------------------------------------------" << endl; - } - //header->PrintParameters(); - - // 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 - - 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()); + 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); 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++) - { - fastjet::PseudoJet mPart=constituents[i]; - ind[i]=mPart.user_index(); - // cout<At(ii); - if(uArray->GetUnitEnergy()>0.){ - uArray->SetUnitTrackID(ipart);//used to have the index available in AliJEtBkg - if(ipart==ind[i]){ - TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef(); - Int_t tracksInCell = trackArray->GetEntries(); - for(int ji = 0; ji < tracksInCell; ji++){ - AliAODTrack * track = (AliAODTrack*)trackArray->At(ji); - aodjet.AddTrack(track); - } - - count++; - } - ipart++; - } - } - } - - - - AddJet(aodjet); - - } // end loop for jets - } + fastjet::PseudoJet mPart=constituents[i]; + ind[i]=mPart.user_index(); + // cout<GetEntries(); + + //internal loop over all the unit cells + Int_t ipart = 0; + + for(Int_t ii=0; iiAt(ii); + if(uArray->GetUnitEnergy()>0.){ + uArray->SetUnitTrackID(ipart);//used to have the index available in AliJEtBkg + if(ipart==ind[i]){ + TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef(); + Int_t tracksInCell = trackArray->GetEntries(); + for(int ji = 0; ji < tracksInCell; ji++){ + AliAODTrack * track = (AliAODTrack*)trackArray->At(ji); + aodjet.AddTrack(track); + } + + count++; + } + ipart++; + } + } + } + } + } + + AddJet(aodjet); + + + } // End loop on jets } @@ -425,17 +434,19 @@ Bool_t AliFastJetFinder::ProcessEvent() // Jets FindJets(); - fJetBkg->SetHeader(fHeader); - fJetBkg->SetReader(fReader); - 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); + if( fAODEvBkg){ + fJetBkg->SetHeader(fHeader); + fJetBkg->SetReader(fReader); + 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); + } @@ -479,16 +490,17 @@ Bool_t AliFastJetFinder::ProcessEvent2() // Jets FindJets(); - fJetBkg->SetHeader(fHeader); - fJetBkg->SetReader(fReader); - fJetBkg->SetFastJetInput(fInputFJ); - Double_t sigma1,meanarea1,sigma2,meanarea2; - Double_t bkg1,bkg2; - fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1); - fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2); - fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1); - fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2); - + if( fAODEvBkg){ + fJetBkg->SetHeader(fHeader); + fJetBkg->SetReader(fReader); + fJetBkg->SetFastJetInput(fInputFJ); + Double_t sigma1,meanarea1,sigma2,meanarea2; + Double_t bkg1,bkg2; + fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1); + fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2); + fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1); + fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2); + } // Double_t bkg1=fJetBkg->BkgFastJet();