// COMET software includes #include "IMidasBank.hxx" #include "IMidasFile.hxx" #include "IMidasBankProxy.hxx" #include "IMidasBankProxyRegistry.hxx" #include "ICOMETRawEvent.hxx" // oaRawEvent includes #include "ITripTHistoBank.hxx" #include "IMidasTripTHistoItr.hxx" #include "IMidasTripTHisto.hxx" // ROOT includes #include "TApplication.h" #include "TFile.h" #include "TH1F.h" #include "TH2F.h" #include "TString.h" #include "TSystem.h" #include #include #include #include #include // Variables. If we turn this into a class, these will be private members static TH2F* gH2Hi = NULL; static TH2F* gH2Lo = NULL; static TFile* gMyFile = NULL; //____________________________________________________________________________ void Book() { gH2Hi = new TH2F("Hi","High gain LI histograms",30*72,0.,30*72.,1024,0.,1024.); gH2Hi->SetOption("COLZ"); gH2Hi->GetXaxis()->SetTitle("Channel + 18*(TRIP + 4*TFB)"); gH2Hi->GetYaxis()->SetTitle("ADC value"); gH2Hi->SetStats(kFALSE); gH2Lo = new TH2F("Lo","Low gain LI histograms",30*72,0.,30*72.,1024,0.,1024.); gH2Lo->SetOption("COLZ"); gH2Lo->GetXaxis()->SetTitle("Channel + 18*(TRIP + 4*TFB)"); gH2Lo->GetYaxis()->SetTitle("ADC value"); gH2Lo->SetStats(kFALSE); } //____________________________________________________________________________ void Event(COMET::ICOMETRawEvent* re) { // Loop over all banks of type ITripTHistoBank COMET::IHandle hBank; while ( hBank = re->GetMidasBank("",hBank) ) { const Char_t *name = hBank->GetName(); // hBank->Print("vv"); Int_t rmm = hBank->GetRMM(); Int_t capacitor = hBank->GetCapacitor(); // Useful Getter methods of hBank... // UInt_t GetRMM() const { // What's the RMM number - also in the bankname, the third character. // UInt_t GetAcquisitionStart() // Get unix time of start of aquisition of these histograms // UInt_t GetAcquisitionEnd() // Get unix time of end of aquisition of these histograms // UInt_t GetHiLo() // Returns 0=high or 1=low for whether these are high or low gain channels // UInt_t GetCapacitor() // Returns capacitor number 0-23. I forget what this returns if the // // capacitors are all in one histogram. 24 possibly? Let me know what you want. // UInt_t GetRecordSeqNo() // Sequence number. The dpt keeps writing as many histograms as will // // fit in each block (they are compressed, so are all different lengths) // // It will always start a new block for each capacitor and each of HiLo. // // and for each RMM separately. // UInt_t GetType() // Return the type of histogram 0=beam pedestal, 1=cosmic pedestal, 2=light injection // // (This is also in the last letter of the histogram bank name B, C or L) // // The record sequences are independent for each of B, C, L although you can // // get them coming interleaved. // GetStartChanIndex(), GetEndChanIndex() ... see blurb below // UInt_t GetNExtraChan() const { // UInt_t GetMaxBinDef() const { return RAWDPT_HVS_BINDEFMAX; } // void GetBinDef(UInt_t *bd) const; // UInt_t DecodeBins(UInt_t *le, UInt_t *w) const; // UInt_t DecodeBins(UInt_t *le, UInt_t *w, const UInt_t *bd) const; // // You have to assume that the histograms will come out in an almost random order: // The only thing which is sure is that for a given Type (B,C or L) and RMM, // it will write out the Hi histograms first with the ind=(tfb*4+trip)*nchanx+chan // variable increasing, until it gets to tfb=47,trip=3,chan=nchanx-1, then it // will do the Lo gain histos in the same order, then (if there is nore than one // capacitor) will do the capacitors in order. The GetRecordSeqNo() keeps incrementing // as you go through this sequence (until it pegs at 255 as it is packed in 8 bits // ). GetStartChanIndex() and GetEndChanIndex() give the first and // last values of ind in this bank (which get looped over in this order in the // inner loop. // // Now what was that about random order? ... Well, the sequence defined above // may get interleaved with other sequences. The RMM and the type ( B, C, L ) // are treated entirely independently within the DPT and the DPTs running // on each FPN are also independent of each other, so these all get interleaved // sometimes. // We should now use GetMaxBinDef and GetBinDef to check that the bin definition is what we expect... int bdmax = hBank->GetMaxBinDef(); // This returns a constant which is just the // size of array you need to allocate on the next line UInt_t bd[bdmax]; hBank->GetBinDef(bd); // This fills the array with the bin definition sequence // It is the same list of constants as in the ODB BinBorder variables // We could now use UInt_t DecodeBins(UInt_t *le, UInt_t *w, const UInt_t *bd) const; to make // arrays of the low edge and width of each bin. Instead, just print bd for now. printf("Got a bank %s Acquired between %d and %d HiLo %d Capacitor %d bdmax = %d:" ,hBank->GetName(),hBank->GetAcquisitionStart(),hBank->GetAcquisitionEnd(),hBank->GetHiLo(),capacitor,bdmax); for (int ii=0;ii= 1023) break; // Sequence terminates either when we have defined up to ADC 1023... if (ii != 0 && bd[ii] == 0 && bd[ii-1] == 0) break; // .. or when we get (0,0) } printf("\n"); // TODO: Now check bd is what we expected // Abbey+Alfons: For now, I have only carefully tested data where the number of bins is 1024 // and so the bin width is 1 everywhere in the histogram. This is straightforward // to fit. The way the DPT allows the histogram bin definition to be done is // simple but probably too flexible - i.e. I propose that you should not write // offline code which can fit a pedestal from every single configuration of the // BinBorder ODB variablke the user could choose - we should discuss some restriction // to make your coding easier. I suggest that we restrict the bin definitions so // you only ever have to deal with a region of bins where the width is 1, although // it may not be the whole range (e.g. 1 bin from 0 to 127, bins of width 1 from // 128 to 255, then a small number of wide bins from there onwards would be OK, // and we should let the ODB variable choose exactly where the region of width // 1 is). Let me know what you think once you have had a go with the histograms // with 1024 bins - we should probably decide a test which the DPT will do and will // warn if configured in a way which makes pedestal determination impossible. // Create an iterator over histograms COMET::IMidasTripTHistoItr hItr(hBank->GetMidasTripTHistoItr()); while ( ! hItr.EOD() ) { COMET::IMidasTripTHisto his(hItr.Get()); UInt_t tfb = his.GetTFB(); UInt_t trip = his.GetTrip(); UInt_t chan = his.GetChan(); const UInt_t *hbins = his.GetDataPtr(); // Note this pointer becomes invalid when the his object disappears Int_t tot = 0; double mean = 0; Int_t nbins = his.GetNBins(); for (int i=0; i 0) mean = mean/tot; if (tot > 0) { // This printout guesses a few bins where the pdestal may be printf("Got a histogram %c%c%c%c %d %d %d %2d cap %d tot %d mean%6.1f bins150-154[%d %d %d %d %d]\n" ,name[0],name[1],name[2],name[3] ,rmm,tfb,trip,chan,capacitor,tot,mean,hbins[150],hbins[151],hbins[152],hbins[153],hbins[154]); // hbins[0] ... [...] contains the contents of the histigram int nchanx = 16 + hBank->GetNExtraChan(); // There are 16 real channels on each Trip. There are two more channels which // are not externally cohnnected, so this measures the electronic pedestal // only (Alfons knows more). At the moment, nchanx will be 18, and extra channel#1 // will be in the high gain set of histograms as channel 16 (as returned by // his.GetChan();) extra channel#2's histogram appears in two places - either // in the high gain set of histograms as channel 17 and also in the low gain // set as channel number 16. This is historical - back in January we were only // just realising we needed the DPT to histogram the low gain channels as well // and memory looked tight. Now it is obvious that we will be doing the // low-gain channels inthe dpt, I suggest we just make nchanx be 17 to save // space and extra channel#2 will be only on channel 16 in the low gain histograms. // Abbey+Alfons: Let me know when you habe thouight bout it if this would be OK // The change is very simple in the DPT code. int ind = (tfb *4+trip)*nchanx+chan; for (int i=0;i<1024;i++) { // TODO: For now I know there are 1024 bins, should calculate from bindef if (hBank->GetHiLo() == 0) { // Hi gH2Hi->SetBinContent(ind,i,hbins[i]); } else { // Lo gH2Lo->SetBinContent(ind,i,hbins[i]); } } } } // End of loop over histograms in this bank. his object disappears here } // End of loop over banks of histograms in this event } //____________________________________________________________________________ void ProcessFile(const char *FileName) { // Start by listing the available Access Classes. COMETLog("Access Classes available to interpret the MIDAS banks"); COMET::IMidasBankProxyRegistry::Instance().Print(); // Make sure file exists FileStat_t fs; if ( gSystem->GetPathInfo(FileName,fs) ) { std::cerr << "Cannot find file: " << FileName << std::endl; return; } COMET::IMidasFile mf; mf.Open(FileName); // Loop over events in this file while ( COMET::ICOMETRawEvent* re = mf.ReadRawEvent() ) { re->PromoteMidasBanks(false); Event(re); delete re; } // End loop over events } //____________________________________________________________________________ int main(int argc, char **argv) { // Get the filename from argv[1] before TApplication modifies the // argument list!!! if (argc != 2) { std::cerr << "usage: " << argv[0] << " " << std::endl; } char FileName[1000]; strncpy(FileName,argv[1],1000); printf("Filename %s\n",FileName); // Start the root system TApplication *app = new TApplication("AccessHisto", &argc, argv); gROOT->SetBatch(); if (app) ; // This just removes 'unused variable app' compiler warning // Open an output root file. gMyFile = new TFile("AccessHisto.root","RECREATE","Demonstration ROOT file"); // Book the histograms Book(); // Analyse the data ProcessFile(FileName); // Close root histogram output file gMyFile->Write(); gMyFile->Close(); delete gMyFile; gMyFile = NULL; return 0; }