#include <IOADatabase.hxx> #include <cometEventLoop.hxx> #include <ITFBChannelId.hxx> #include <ITripTDigitBank.hxx> #include <ITripTXTimeBank.hxx> #include <IMidasTripTDigitItr.hxx> #include <IMidasTripTDigit.hxx> #include <IRawDataHeader.hxx> #include <ITripTInfoBank.hxx> #include <IGeometryDatabase.hxx> #include <IMCMBank.hxx> #include "P0D_GainAnalysis.hxx" using namespace COMET; bool HasGoodTDC(const COMET::IMidasTripTDigit& mDigit) { UInt_t time_offset = mDigit.GetTimeOffset(); Bool_t err = (( time_offset == RAWDPT_DVX_FMT0_TIMEOVERFLOW) || (time_offset == RAWDPT_DVX_FMT0_TIMENOHIT) || (time_offset == RAWDPT_DVX_FMT0_TIMEMAX) ); return !err; } class IP0dHealth: public COMET::ICOMETEventLoopFunction { public: IP0dHealth() :fDB(COMET::IGeometryDatabase::Get()), fTFBs(300, -1), fReadEventCount(0){ } virtual ~IP0dHealth() {}; void Usage(void) { std::cout << " -O <option> Set the option for TObject::ls()" << std::endl; std::cout << " -O filter=<value> Sets event filter type" << std::endl; } void Initialize(); void Finalize(ICOMETOutput * const output); using COMET::ICOMETEventLoopFunction::Finalize; virtual bool SetOption(std::string option,std::string value="") { if (option == "filter"){ fFilterType=atoi(value.c_str()); fFilter=true; return true; } if (option == "outputfile" ){ fOutputFileName = value; return true; } return false; } bool operator () (COMET::ICOMETEvent& event); private: COMET::IGeometryDatabase& fDB; std::vector<int> fTFBs; int fReadEventCount; bool fAllocated; //Raw ADC Map Histograms TH2D* fBarHighADCMap; TH2F* fBarLowADCMap; //Timing Histograms TH2I* fBarTime; TH1I* fBarTrig; TH2I* fBarTimeADCHigh; TH2I* fBarTimeADCLow; TH2F* fTDCHitHigh; TH2F* fTDCHitLow; //Histograms based on electronic names (not p0dule and channel) TH2F* fTripTADCHigh; TH2F* fTripTADCLow; TH2F* fTripTADCHighTDC; TH2F* fTripTADCLowTDC; //Histogram of trigger types TH1I* fTrigger; bool fFilter; bool fStart; int fRunId; unsigned int fFilterType; int fEventCounter; int fGroupCounter; std::vector<int> fGroupEvents; std::vector<unsigned int> fGroupTime; TFile* fOutputFile; std::string fOutputFileName; }; void IP0dHealth::Initialize() { TH1::AddDirectory(false); fTrigger = new TH1I("Triggers","Triggers",256,0,256); fEventCounter=0; fGroupCounter=0; fStart=true; } void IP0dHealth::Finalize(ICOMETOutput* const output){ if(fEventCounter<1000 && fEventCounter != 0){ COMETLog("Last Data point has " << fEventCounter << " events in it. Use with caution!!"); fGroupEvents.push_back(fEventCounter); fBarHighADCMap->Write(); fBarLowADCMap->Write(); fBarTime->Write(); fBarTrig->Write(); fBarTimeADCHigh->Write(); fBarTimeADCLow->Write(); fTDCHitHigh->Write(); fTDCHitLow->Write(); fTripTADCHigh->Write(); fTripTADCLow->Write(); fTripTADCHighTDC->Write(); fTripTADCLowTDC->Write(); fTrigger->Write(); } int time[fGroupCounter],events[fGroupCounter]; for(int i=0;i<fGroupCounter;i++){ std::cout<<fGroupTime.at(i)<<"\t"<<fGroupEvents.at(i)<<"\n"; time[i]=fGroupTime.at(i); events[i]=fGroupEvents.at(i); //call gain analysis for each group GainAnalysis(fOutputFile,time[i]); } TGraph* GroupTiming = new TGraph(fGroupCounter,time,events); GroupTiming->SetMarkerStyle(8); GroupTiming->Draw("AP"); GroupTiming->GetXaxis()->SetTitle("Seconds since first group"); GroupTiming->GetYaxis()->SetTitle("Number of events in the group"); GroupTiming->SetTitle("Grouping Timing vs Events per group"); GroupTiming->Write("Grouping Timing vs Events per group"); fOutputFile->Close(); } bool IP0dHealth::operator () (COMET::ICOMETEvent& event) { unsigned int icurTime; const COMET::IP0DChannelMap& P0DChanMap = fDB.P0D(); if(fFilter!=true){ COMETLog("If you are seeing this message please note you must choose a trigger type to filter on using the -O filter=<value> where the following values are used:\n1=beam\n2=pedestal\n8=P0DLED\n128=Cosmic"); } fRunId = event.GetContext().GetRun(); const int eventId = event.GetContext().GetEvent(); if(eventId%50 == 0){ COMETLog("\t Run " << fRunId << " event "<<eventId); } if(fStart==true && fRunId!=eventId){ COMETVerbose("Filter type is " << fFilterType ); fOutputFile = new TFile( fOutputFileName != "" ? fOutputFileName.c_str() : Form("p0d-health-%d-filter-%d.root", fRunId,fFilterType), "RECREATE"); fStart=false; } fTrigger->Fill(event.GetHeader().GetTriggerBits()); if(event.GetHeader().GetTriggerBits()==fFilterType){ COMETVerbose("Got correct trigger:" << event.GetHeader().GetTriggerBits() ); if(fEventCounter==0 || fEventCounter == 1000){ icurTime=event.GetContext().GetTimeStamp(); fGroupTime.push_back(icurTime); fBarHighADCMap = new TH2D(Form("HighADCMap-%d",icurTime),Form("High ADC Map-%d",fGroupCounter), 13002, -2, 13000, 1050, 0, 1050); fBarHighADCMap->GetXaxis()->SetTitle("f(bar,layer,p0dule,capacitor)"); fBarHighADCMap->GetYaxis()->SetTitle("High Gain ADC"); fBarLowADCMap = new TH2F(Form("LowADCMap-%d",icurTime),Form("Low ADC Map-%d",fGroupCounter), 13002, -2, 13000, 1050, 0, 1050); fBarLowADCMap->GetXaxis()->SetTitle("f(bar,layer,p0dule,capacitor)"); fBarLowADCMap->GetYaxis()->SetTitle("Low Gain ADC"); fBarTime = new TH2I(Form("TFB vs Time-%d",icurTime),Form("TFBvsTime-%d",fGroupCounter),12000,-6000,6000,400,0,400); fBarTime->GetXaxis()->SetTitle("Time (2.5ns per bin)"); fBarTime->GetYaxis()->SetTitle("48*RMM+TFB"); fBarTrig = new TH1I(Form("Hit Timing-%d",icurTime),Form("Hit Timing-%d",fGroupCounter),12000,-6000,6000); fBarTrig->GetXaxis()->SetTitle("Time (2.5ns per bin"); fBarTrig->GetYaxis()->SetTitle("N"); fBarTimeADCHigh = new TH2I(Form("High Gain ADC vs Time-%d",icurTime),Form("ADCvsTime-%d",fGroupCounter),12000,-6000,6000,1050,0,1050); fBarTimeADCHigh->GetXaxis()->SetTitle("Time (2.5ns per bin)"); fBarTimeADCHigh->GetYaxis()->SetTitle("High Gain ADC"); fBarTimeADCLow = new TH2I(Form("Low Gain ADC vs Time-%d",icurTime), Form("ADCvsTime-%d",fGroupCounter),12000,-6000,6000,1050,0,1050); fBarTimeADCLow->GetXaxis()->SetTitle("Time (2.5ns per bin)"); fBarTimeADCLow->GetYaxis()->SetTitle("Low Gain ADC"); fTDCHitHigh = new TH2F(Form("hiADC with TDC-%d",icurTime),Form("hiADC with TDC-%d",fGroupCounter),13002,-2,13000,1050,0,1050); fTDCHitHigh->GetXaxis()->SetTitle("f(bar,layer,p0dule,capacitor)"); fTDCHitHigh->GetYaxis()->SetTitle("High Gain ADC"); fTDCHitLow = new TH2F(Form("loADC with TDC-%d",icurTime),Form("loADC with TDC-%d",fGroupCounter),13002,-2,13000,1050,0,1050); fTDCHitLow->GetXaxis()->SetTitle("f(bar,layer,p0dule,capacitor)"); fTDCHitLow->GetYaxis()->SetTitle("Low Gain ADC"); fTripTADCHigh = new TH2F(Form("TripT raw ADC high gain-%d",icurTime),Form("TripT raw ADC high gain-%d",fGroupCounter),2000,0,2000,1050,0,1050); fTripTADCHigh->GetXaxis()->SetTitle("RMM*48*4+port*4+trip"); fTripTADCHigh->GetYaxis()->SetTitle("PE"); fTripTADCLow = new TH2F(Form("TripT raw ADC low gain-%d",icurTime),Form("TripT raw ADC low gain-%d",fGroupCounter),2000,0,2000,1050,0,1050); fTripTADCLow->GetXaxis()->SetTitle("RMM*48*4+port*4+trip"); fTripTADCLow->GetYaxis()->SetTitle("PE"); fTripTADCHighTDC = new TH2F(Form("TripT raw ADC high gain valid TDC-%d",icurTime),Form("TripT raw ADC high gain valid TDC-%d",fGroupCounter),2000,0,2000,1050,0,1050); fTripTADCHighTDC->GetXaxis()->SetTitle("RMM*48*4+port*4+trip"); fTripTADCHighTDC->GetYaxis()->SetTitle("PE"); fTripTADCLowTDC = new TH2F(Form("TripT raw ADC low gain valid TDC-%d",icurTime),Form("TripT raw ADC low gain valid TDC-%d",fGroupCounter),2000,0,2000,1050,0,1050); fTripTADCLowTDC->GetXaxis()->SetTitle("RMM*48*4+port*4+trip"); fTripTADCLowTDC->GetYaxis()->SetTitle("PE"); fGroupCounter +=1; fEventCounter=0; } fEventCounter += 1; IHandle<ICOMETRawEvent> rawEvt(event.GetRaw()); assert(rawEvt); const int first = 0; const int last = 22; int trigtime = 0XDEADBEEF; //int timesincebusy_mcm,evttime; IHandle<IMCMBank> mcmInfo; while(mcmInfo = rawEvt->GetMidasBank<COMET::IMCMBank>("", mcmInfo)){ std::stringstream rr; rr << mcmInfo->GetName(); std::string name(rr.str()); std::string p0d="P"; if(name[0]!=p0d) continue; trigtime = mcmInfo->GetUnixTimeSSecTrig(); //timesincebusy_mcm = mcmInfo->GetTimeSinceBusy(); //evttime = mcmInfo->GetEventTime(); } IHandle<ITripTDigitBank> triptBank; while (triptBank = rawEvt->GetMidasBank<ITripTDigitBank>("",triptBank) ) { if (triptBank->GetSubDetector() != 1){ // Not a P0D Bank, ignore it COMETVerbose("Not a P0D Bank, ignore it"); continue; } IMidasTripTDigitItr itr(triptBank->GetMidasTripTDigitItr()); while (!itr.EOD() ) { IMidasTripTDigit digit(itr.Get()); ITFBChannelId tfbChannelId(digit.GetChannelId()); const int rmm = digit.GetRMMNum(); const int port = digit.GetTFBNum(); const int tripT = digit.GetTripTNum(); //const int channel = digit.GetChannelNum(); //Unused const int capacitor = digit.GetIntegrationNum(); const int time_tfb = digit.GetTime(); //const int time_offset = digit.GetTimeOffset(); //Unused //ignore first and last integration cycle as they tend to shift the pedestal values (in error) if (capacitor == first || capacitor == last) continue; //const int id = (rmm*48 + port)*64 + tripT*16 + channel; //const int id = tfbChannelId.GetCableId(); COMET::IGeometryId geometryId; if (!fDB.GetGeometryId(geometryId, tfbChannelId)) { COMETVerbose("Skip: !GetGeometryId()"); continue; } int p0dule, layer, bar; if(!P0DChanMap.GeomIDToLocation(p0dule, layer, bar, geometryId)) { COMETVerbose("Skip: !GeomIDToLocation()"); continue; } int position = p0dule*140+50+bar+layer*6000; int timeRelToTrig_tfb = (time_tfb)%400000000-trigtime*4; //Fill raw ADC maps fBarHighADCMap->Fill(position, digit.GetHighGainADC()); fBarLowADCMap->Fill(position, digit.GetLowGainADC()); fTripTADCHigh->Fill(rmm*48*4+port*4+tripT,digit.GetHighGainADC()); fTripTADCLow->Fill(rmm*48*4+port*4+tripT,digit.GetLowGainADC()); //Fill timing histograms fBarTime->Fill(timeRelToTrig_tfb,rmm*48+port); fBarTrig->Fill(timeRelToTrig_tfb); fBarTimeADCHigh->Fill(timeRelToTrig_tfb,digit.GetHighGainADC()); fBarTimeADCLow->Fill(timeRelToTrig_tfb,digit.GetLowGainADC()); //Fill valid TDC ADC histograms if(HasGoodTDC(digit)){ fTDCHitHigh->Fill((p0dule*140+50+bar+6000*layer),digit.GetHighGainADC()); fTDCHitLow->Fill((p0dule*140+50+bar+6000*layer),digit.GetLowGainADC()); fTripTADCHighTDC->Fill((rmm*48*4+port*4+tripT),digit.GetHighGainADC()); fTripTADCLowTDC->Fill((rmm*48*4+port*4+tripT),digit.GetLowGainADC()); } } } if(fEventCounter==1000){ fGroupEvents.push_back(fEventCounter); fBarHighADCMap->Write(); fBarLowADCMap->Write(); fBarTime->Write(); fBarTrig->Write(); fBarTimeADCHigh->Write(); fBarTimeADCLow->Write(); fTDCHitHigh->Write(); fTDCHitLow->Write(); fTripTADCHigh->Write(); fTripTADCLow->Write(); fTripTADCHighTDC->Write(); fTripTADCLowTDC->Write(); delete fBarHighADCMap; delete fBarLowADCMap; delete fBarTime; delete fBarTrig; delete fBarTimeADCHigh; delete fBarTimeADCLow; delete fTDCHitHigh; delete fTDCHitLow; delete fTripTADCHigh; delete fTripTADCLow; delete fTripTADCHighTDC; delete fTripTADCLowTDC; } } ++fReadEventCount; return true; } int main(int argc, char **argv) { IP0dHealth userCode; cometEventLoop(argc,argv,userCode,1); }