#include #include #include IFgdCollinearCosmicControlSample::IFgdCollinearCosmicControlSample() : IControlSampleBase("FgdCollinearCosmic"){ AddValidTrigger(COMET::Trigger::kFGDCosmic); AddValidTrigger(COMET::Trigger::kTFBCosmic); // Get Prescale factor from parameters file int XZor3D = COMET::IOARuntimeParameters::Get().GetParameterI("oaControlSample.FgdCollinearCosmic.XZor3D"); SetXZor3D(XZor3D); //Get Collinearity cut double collCut = COMET::IOARuntimeParameters::Get().GetParameterD("oaControlSample.FgdCollinearCosmic.CollCut"); SetCollCut(collCut); // fPrescaleUy = COMET::IOARuntimeParameters::Get().GetParameterI("oaControlSample.FgdCollinearCosmic.PrescaleUy"); if(fPrescaleUy){ InitPrescaleUy(); } } void IFgdCollinearCosmicControlSample::InitHistograms(){ hUX1 = new TH1F("hUX1", ";U_{X1}; Events", 50, -1.0, 1.0); hUX2 = new TH1F("hUX2", ";U_{X2}; Events", 50, -1.0, 1.0); hUXAvg = new TH1F("hUXAvg", ";Average U_{X}; Events", 50, -1.0, 1.0); hUY1 = new TH1F("hUY1", ";U_{Y1}; Events", 50, -1.0, 1.0); hUY2 = new TH1F("hUY2", ";U_{Y2}; Events", 50, -1.0, 1.0); hUYAvg = new TH1F("hUYAvg", ";Average U_{Y}; Events", 50, -1.0, 1.0); hUZ1 = new TH1F("hUZ1", ";U_{Z1}; Events", 50, -1.0, 1.0); hUZ2 = new TH1F("hUZ2", ";U_{Z2}; Events", 50, -1.0, 1.0); hUZAvg = new TH1F("hUZAvg", ";Average U_{Z}; Events", 50, -1.0, 1.0); hUX1Fgd = new TH1F("hUX1Fgd", ";U_{X1}; Events", 50, -1.0, 1.0); hUX2Fgd = new TH1F("hUX2Fgd", ";U_{X2}; Events", 50, -1.0, 1.0); hUXAvgFgd = new TH1F("hUXAvgFgd", ";Average U_{X}; Events", 50, -1.0, 1.0); hUY1Fgd = new TH1F("hUY1Fgd", ";U_{Y1}; Events", 50, -1.0, 1.0); hUY2Fgd = new TH1F("hUY2Fgd", ";U_{Y2}; Events", 50, -1.0, 1.0); hUYAvgFgd = new TH1F("hUYAvgFgd", ";Average U_{Y}; Events", 50, -1.0, 1.0); hUZ1Fgd = new TH1F("hUZ1Fgd", ";U_{Z1}; Events", 50, -1.0, 1.0); hUZ2Fgd = new TH1F("hUZ2Fgd", ";U_{Z2}; Events", 50, -1.0, 1.0); hUZAvgFgd = new TH1F("hUZAvgFgd", ";Average U_{Z}; Events", 50, -1.0, 1.0); hUX1Tfb = new TH1F("hUX1Tfb", ";U_{X1}; Events", 50, -1.0, 1.0); hUX2Tfb = new TH1F("hUX2Tfb", ";U_{X2}; Events", 50, -1.0, 1.0); hUXAvgTfb = new TH1F("hUXAvgTfb", ";Average U_{X}; Events", 50, -1.0, 1.0); hUY1Tfb = new TH1F("hUY1Tfb", ";U_{Y1}; Events", 50, -1.0, 1.0); hUY2Tfb = new TH1F("hUY2Tfb", ";U_{Y2}; Events", 50, -1.0, 1.0); hUYAvgTfb = new TH1F("hUYAvgTfb", ";Average U_{Y}; Events", 50, -1.0, 1.0); hUZ1Tfb = new TH1F("hUZ1Tfb", ";U_{Z1}; Events", 50, -1.0, 1.0); hUZ2Tfb = new TH1F("hUZ2Tfb", ";U_{Z2}; Events", 50, -1.0, 1.0); hUZAvgTfb = new TH1F("hUZAvgTfb", ";Average U_{Z}; Events", 50, -1.0, 1.0); hUXAvgPrescaled = new TH1F("hUXAvgPrescaled", ";Average U_{X}; Events", 50, -1.0, 1.0); hUYAvgPrescaled = new TH1F("hUYAvgPrescaled", ";Average U_{Y}; Events", 50, -1.0, 1.0); hUZAvgPrescaled = new TH1F("hUZAvgPrescaled", ";Average U_{Z}; Events", 50, -1.0, 1.0); } void IFgdCollinearCosmicControlSample::InitPrescaleUy(){ // get the root file name std::string fFileName = std::string(getenv("OACONTROLSAMPLEROOT")) + std::string("/dat/") + COMET::IOARuntimeParameters::Get().GetParameterS("oaControlSample.FgdCollinearCosmic.PrescaleUyHistoFile"); fUyHistoFile = new TFile(fFileName.c_str()); if(!fUyHistoFile){ COMETError("\n***** COMET::IFgdCollinearCosmicControlSample *****\n" << "Cannot open Uy Histogram file '" << fFileName.c_str() << "'."); } // get the histogram name std::string fHistoName = COMET::IOARuntimeParameters::Get().GetParameterS("oaControlSample.FgdCollinearCosmic.PrescaleUyHisto"); fUyHisto = (TH1F*)fUyHistoFile->Get(fHistoName.c_str()); if(!fUyHisto){ COMETError("\n***** COMET::IFgdCollinearCosmicControlSample *****\n" << "Cannot find Uy Histogram '" << fHistoName.c_str() << "'."); } // get the cut value fUyPcut = COMET::IOARuntimeParameters::Get().GetParameterD("oaControlSample.FgdCollinearCosmic.PrescaleUyCut"); return; } bool IFgdCollinearCosmicControlSample::IsEventSelectedInternal(COMET::ICOMETEvent& event) { COMET::IHandle FgdResult = event.GetFit("ReconFGD"); if(!FgdResult) return false; COMET::IHandle TrackCont = FgdResult->GetResultsContainer("fittedFgdReconTracks_allFGD"); if(!TrackCont) return false; COMET::IHandle bin0result = fgdUtils::getTimeBinResult(event,0); COMET::IHandle recon = bin0result->GetResultsContainer("fittedFgdReconTracks_allFGD"); int nfgd1 =0, nfgd2 = 0; COMET::IHandle fgd1Track; COMET::IHandle fgd2Track; for( COMET::IReconObjectContainer::iterator tr = recon->begin(); tr < recon->end();tr++) { COMET::IHandle track = *tr; if(track->UsesDetector(COMET::IReconBase::kFGD1)){ nfgd1++; fgd1Track = track; } if(track->UsesDetector(COMET::IReconBase::kFGD2)){ nfgd2++; fgd2Track = track; } } // we want one track from the 0th time bin in each FGD if(nfgd1 != 1 || nfgd2 != 1) return false; // now we have one track in each FGD in the 0th time bin // determine collinearity COMET::IHandle fgd1State = fgd1Track->GetState(); COMET::IHandle fgd2State = fgd2Track->GetState(); if(!fgd1State || !fgd2State) return false; TVector3 u1 = fgd1State->GetDirection(); TVector3 u2 = fgd2State->GetDirection(); double u1u2 = -999; if(fXZor3D == 0){ // dot product in Ux/Uz Double_t u1x_2d = u1.X()/sqrt(u1.X()*u1.X() + u1.Z()*u1.Z()); Double_t u1z_2d = u1.Z()/sqrt(u1.X()*u1.X() + u1.Z()*u1.Z()); Double_t u2x_2d = u2.X()/sqrt(u2.X()*u2.X() + u2.Z()*u2.Z()); Double_t u2z_2d = u2.Z()/sqrt(u2.X()*u2.X() + u2.Z()*u2.Z()); u1u2 = u1x_2d * u2x_2d + u1z_2d * u2z_2d; } else { // dot product in 3D u1u2 = u1.Dot(u2); } bool passCollCut = (u1u2 > fCollCut); TVector3 uavg = 0.5*(u1+u2); bool fgdCosmic = (event.GetHeader().GetTriggerBits() & COMET::Trigger::kFGDCosmic); bool tfbCosmic = (event.GetHeader().GetTriggerBits() & COMET::Trigger::kTFBCosmic); if(passCollCut && fMakeHistoFile){ hUX1->Fill(u1.X()); hUY1->Fill(u1.Y()); hUZ1->Fill(u1.Z()); hUX2->Fill(u2.X()); hUY2->Fill(u2.Y()); hUZ2->Fill(u2.Z()); hUXAvg->Fill(uavg.X()); hUYAvg->Fill(uavg.Y()); hUZAvg->Fill(uavg.Z()); if(fgdCosmic){ hUX1Fgd->Fill(u1.X()); hUY1Fgd->Fill(u1.Y()); hUZ1Fgd->Fill(u1.Z()); hUX2Fgd->Fill(u2.X()); hUY2Fgd->Fill(u2.Y()); hUZ2Fgd->Fill(u2.Z()); hUXAvgFgd->Fill(uavg.X()); hUYAvgFgd->Fill(uavg.Y()); hUZAvgFgd->Fill(uavg.Z()); } else if(tfbCosmic){ hUX1Tfb->Fill(u1.X()); hUY1Tfb->Fill(u1.Y()); hUZ1Tfb->Fill(u1.Z()); hUX2Tfb->Fill(u2.X()); hUY2Tfb->Fill(u2.Y()); hUZ2Tfb->Fill(u2.Z()); hUXAvgTfb->Fill(uavg.X()); hUYAvgTfb->Fill(uavg.Y()); hUZAvgTfb->Fill(uavg.Z()); } } // prescale on Uy only for FGD Cosmics if(fPrescaleUy && fgdCosmic){ double uyEntries = fUyHisto->GetBinContent(fUyHisto->FindBin(uavg.Y())); double uyFrac = uyEntries/fUyHisto->Integral(); if(uyFrac != 0){ passCollCut = passCollCut && (fRandom.Rndm() < fUyPcut/uyFrac); if(passCollCut && fMakeHistoFile){ hUXAvgPrescaled->Fill(uavg.X()); hUYAvgPrescaled->Fill(uavg.Y()); hUZAvgPrescaled->Fill(uavg.Z()); } } } return passCollCut; }