/////////////////////////////////////////////////////////////////////////////// // // Summary of fPrintLevel values: // ---------------------- // fPrintLevel = 1: print Electron block for all the events // fPrintLevel = 101: print line per event with good electron // fPrintLevel = 102: found matching non-trigger jet /////////////////////////////////////////////////////////////////////////////// #include "TF1.h" #include "TCanvas.h" #include "TPad.h" #include "TText.h" #include "Stntuple/obj/TStnHeaderBlock.hh" #include "Stntuple/obj/TStnTrackBlock.hh" #include "Stntuple/obj/TStnDBManager.hh" #include "Stntuple/alg/TStntuple.hh" #include "Stntuple/alg/TStnElectronID.hh" //#include "Stntuple/alg/TStnMuonID.hh" #include "Stntuple/obj/TStnTriggerTable.hh" #include "Stntuple/obj/TStnTrigger.hh" #include "Stntuple/loop/TStnAna.hh" #include "Stntuple/obj/TStnTrack.hh" #include "Stntuple/geom/TTrajectoryPoint.hh" #include "Stntuple/geom/TSimpleExtrapolator.hh" #include "TJpsiMonModule.hh" // #include "TTauUtils.hh" #include "Stntuple/fit/TCtvmft.hh" ClassImp(TJpsiMonModule) //_____________________________________________________________________________ TJpsiMonModule::TJpsiMonModule(const char* name, const char* title): TStnModule(name,title) { // fListOfGoodElectrons = new TObjArray(20); // fListOfLooseElectrons = new TObjArray(20); fListOfGoodMuons = new TObjArray(20); fListOfLooseMuons = new TObjArray(20); // dimuFit = new TCtvmft(); // fTightID = new TStnElectronID(); // fLooseID = new TStnElectronID(); //----------------------------------------------------------------------------- // modes: //----------------------------------------------------------------------------- fCutMode = 0; } //_____________________________________________________________________________ TJpsiMonModule::~TJpsiMonModule() { fListOfGoodMuons->Clear(); delete fListOfGoodMuons; fListOfLooseMuons->Clear(); delete fListOfLooseMuons; // delete dimuFit; } //_____________________________________________________________________________ int TJpsiMonModule::BeginJob() { // register the data blocks and book histograms RegisterDataBlock("MuonBlock" , "TStnMuonBlock" , &fMuonBlock); RegisterDataBlock("TriggerBlock" , "TStnTriggerBlock", &fTriggerBlock); RegisterDataBlock("TrackBlock" , "TStnTrackBlock" , &fTrackBlock); RegisterDataBlock("MetBlock" , "TStnMetBlock" , &fMetBlock); BookHistograms(); return 0; } //_____________________________________________________________________________ int TJpsiMonModule::BeginRun() { int rn = GetHeaderBlock()->RunNumber(); TStntuple::Init(rn); TStntuple::InitListOfL3Triggers(this); return 0; } //_____________________________________________________________________________ void TJpsiMonModule::FillEventHistograms(EventHist_t& Hist) { int nmuons = fMuonBlock->NMuons(); Hist.fNMu->Fill(nmuons); } //_____________________________________________________________________________ void TJpsiMonModule::FillMuMuHistograms(MuMuHist_t& Hist, TStnMuon* Mu1 , TStnMuon* Mu2 ) { // Fill simple histograms TLorentzVector dm; dm = *Mu1->Momentum(); dm += *Mu2->Momentum(); double charge = Mu1->Charge()+Mu2->Charge(); double mass = dm.M(); double pt = dm.Pt(); Hist.fMass->Fill(mass); Hist.fCharge->Fill(charge); Hist.fPt->Fill(pt); // Now deal with COT-only (using beam constrained tracks) // and SVX-only int trkmu1_index = Mu1->TrackNumber(); int trkmu2_index = Mu2->TrackNumber(); const TStnTrack* trk1 = fTrackBlock->Track(trkmu1_index); const TStnTrack* trk2 = fTrackBlock->Track(trkmu2_index); // SVX only int N_svxhits1 = trk1->NSvxHits(); int N_svxhits2 = trk2->NSvxHits(); if ((N_svxhits1 > 0) && (N_svxhits2 > 0)){ double svxmass = dm.M(); Hist.fMassSvx->Fill(svxmass); } // mass based on COT beam-constrained tracks static const double crvpt = 0.002116; static const float MU_MASS = 0.105658357; static const float TINY = 3.0e-09; bool goodbc = true; TLorentzVector bcdm; double pt1 = abs(crvpt/trk1->BcC0()); if (pt1 < TINY){ goodbc = false;} else{ double px = pt1*cos(trk1->BcPhi0()); double py = pt1*sin(trk1->BcPhi0()); double pz = pt1*(trk1->BcLam0()); double en = sqrt(pt1*pt1 + pz*pz + MU_MASS*MU_MASS); TLorentzVector p4mu(px,py,pz,en); bcdm = p4mu; } double pt2 = abs(crvpt/trk2->BcC0()); if (pt2 < TINY){ goodbc = false;} else{ double px = pt2*cos(trk2->BcPhi0()); double py = pt2*sin(trk2->BcPhi0()); double pz = pt2*(trk2->BcLam0()); double en = sqrt(pt2*pt2 + pz*pz + MU_MASS*MU_MASS); TLorentzVector p4mu(px,py,pz,en); bcdm += p4mu; } // int tralg1 = trk1->Algorithm(); if (goodbc){ double cotmass = bcdm.M(); Hist.fMassCot->Fill(cotmass); } // if (charge == 0) // printf("good BC: mass,cotmass,trkalg1,trkalg2 %g, %g, %i, %i\n",mass,cotmass,alg1,alg2); // } // else // printf("BAD BC: mass,cotmass,ptbc1,ptbc2,trkalg1,trkalg2:%g, %g, %g, %g, %i, %i,\n", mass,charge,pt1,pt2,alg1,alg2); // To use CTVMFT, uncomment below and also the line in the beginning // of the file created dimuFit and deleted dimuFit // Now do dimuon CTVMFT fit // if (charge == 0){ // dimuFit->Init(); // dimuFit->AddTrack(trk1,MU_MASS,1); // dimuFit->AddTrack(trk2,MU_MASS,1); // bool status = dimuFit->Fit(); // float dmass = -999; // float massfit = -999; // if (status){ // massfit = dimuFit->GetMass(trk1,trk2,dmass); // Hist.fFitMass->Fill(massfit); // } // } // if (mass < 0.5) GetHeaderBlock()->Print(Form(" dm mass = %10.3f\n",mass)); } //_____________________________________________________________________________ void TJpsiMonModule::FillMuonHistograms(MuoHist_t& Hist, TStnMuon* Muo) { int nmuons = fMuonBlock->NMuons(); const TLorentzVector* mom; for (int i=0; iMuon(i); mom = muon->Momentum(); Hist.fDetCode->Fill(muon->DetCode()); Hist.fPt->Fill(mom->Pt()); Hist.fEta->Fill(mom->Eta()); Hist.fPhi->Fill(TVector2::Phi_0_2pi(mom->Phi())); Hist.fEmEnergy->Fill(muon->EmEnergy()); Hist.fHadEnergy->Fill(muon->HadEnergy()); int mu_index = muon->TrackNumber(); TStnTrack* trk = fTrackBlock->Track(mu_index); int tralg = trk->Algorithm(); Hist.fTrkAlg->Fill(trk->Algorithm()); Hist.fNSVXHits->Fill(trk->NSvxHits()); } } //_____________________________________________________________________________ void TJpsiMonModule::FillEScaleHistograms(EScaleHist_t& Hist, TStnMuon* Mu1 , TStnMuon* Mu2 ) { // Select good muon pairs TLorentzVector dm; dm = *Mu1->Momentum(); dm += *Mu2->Momentum(); double mass = dm.M(); if (mass > 3.0 && mass < 3.2 && Mu1->Charge() != Mu2->Charge() && IsGoodMuon(Mu1) && IsGoodMuon(Mu2)) { Hist.fEmEnergyCMU->Fill(Mu1->EmEnergy()); Hist.fEmEnergyCMU->Fill(Mu2->EmEnergy()); Hist.fHadEnergyCMU->Fill(Mu1->HadEnergy()); Hist.fHadEnergyCMU->Fill(Mu2->HadEnergy()); } } //_____________________________________________________________________________ bool TJpsiMonModule::IsGoodMuon(TStnMuon* Mu) { bool good = false; if (Mu->HasCmuStub() && Mu->EmEnergy() < 2. && Mu->HadEnergy() < 6. && Mu->TrackPt() > 1.5) good = true; return good; } //_____________________________________________________________________________ int TJpsiMonModule::Event(int ientry) { Int_t nmu, passed; TStnHeaderBlock *h; h = GetHeaderBlock(); //----------------------------------------------------------------------------- // check the trigger path //----------------------------------------------------------------------------- fTriggerBlock->GetEntry(ientry); if (h->McFlag() == 0) { passed = TStntuple::CheckL3TriggerPath(fTriggerBlock,fListOfL3Triggers); } else { // do not check L3... passed = 1; } //----------------------------------------------------------------------------- // always fill histogram for the run number //----------------------------------------------------------------------------- fHist.fRunNumber->Fill(GetHeaderBlock()->RunNumber()); fHist.fFilterResult->Fill(passed); if (! passed) return 0; //----------------------------------------------------------------------------- // read the necessary data //----------------------------------------------------------------------------- fMuonBlock->GetEntry(ientry); fTrackBlock->GetEntry(ientry); fMetBlock->GetEntry(ientry); // do whatever you want nmu = fMuonBlock->NMuons(); for (int i=0; iMuon(i); FillMuonHistograms(fHist.fMuoHist,muo); } FillEventHistograms(fHist.fEventHist); //----------------------------------------------------------------------------- // 0. tight-(>=loose) pairs, opposite signs //----------------------------------------------------------------------------- for (int i1=0; i1Muon(i1); for (int i2=i1+1; i2Muon(i2); FillEScaleHistograms(fHist.fEScaleHist,mu1,mu2); if (mu1->Charge() != mu2->Charge()) { FillMuMuHistograms(fHist.fMuMuHist[0],mu1,mu2); } else { // e2 is loose FillMuMuHistograms(fHist.fMuMuHist[1],mu1,mu2); } } } return 0; } //_____________________________________________________________________________ int TJpsiMonModule::EndJob() { printf("----- end job: ---- %s\n",GetName()); return 0; } //_____________________________________________________________________________ void TJpsiMonModule::PlotHistograms(int run_number, int slide) { // plot slides char name[120], canvas_name[120], title[120], pad_name[120]; sprintf(name,"run_%i_%i",run_number,slide); sprintf(title,"run %i slide %i ",run_number, slide); sprintf(canvas_name,"%s_%s",GetName(),name); sprintf(pad_name,"%s_p1",canvas_name); if (slide == 1) { //----------------------------------------------------------------------------- // //----------------------------------------------------------------------------- // TPostScript ps("l1ana.ps",-111); TCanvas* c = NewSlide(name,title,2,3); TPad* p1 = (TPad*) c->GetPrimitive(pad_name); // p1->cd(1); // fHist.fEle.fEtVsEta->Draw(); gPad->Update(); } else if (slide == 2) { //----------------------------------------------------------------------------- // //----------------------------------------------------------------------------- // TPostScript ps("l1ana.ps",-111); TCanvas* c = NewSlide(name,title,2,3); TPad* p1 = (TPad*) c->GetPrimitive(pad_name); // p1->cd(1); // fHist.fEle.fIso->Draw(); gPad->Update(); } }