loon: Bdtxt.cxx Source File

00001 #include "Bdtxt.h"
00002 #include "BmntUtil.h"
00003 
00004 #include <TFile.h>
00005 #include <TTree.h>
00006 
00007 #include <Validity/VldContext.h>
00008 #include <RawData/RawRecord.h>
00009 #include <RawData/RawBeamMonHeader.h>
00010 #include <RawData/RawBeamMonBlock.h>
00011 #include <RawData/RawDataBlock.h>
00012 #include <RawData/RawBeamData.h>
00013 #include <RawData/RawBeamSwicData.h>
00014 
00015 #include <cassert>
00016 #include <iostream>
00017 #include <vector>
00018 #include <string>
00019 #include <map>
00020 #include <fstream>
00021 #include <cmath>
00022 using namespace std;
00023 
00024 const char* simple_data[] = {
00025     "E:TOR101",
00026     "E:TORTGT",
00027     "E:HMRTD",
00028     "E:MM1RTD",
00029     "E:MM2RTD",
00030     "E:MM3RTD",
00031     "E:HMGPR",
00032     "E:HMGPD",
00033     "E:HMGF",
00034     "E:MM1GPR",
00035     "E:MM1GPD",
00036     "E:MM1GF",
00037     "E:MM2GPR",
00038     "E:MM2GPD",
00039     "E:MM2GF",
00040     "E:MM3GPR",
00041     "E:MM3GPD",
00042     "E:MM3GF",
00043     "E:HMHV1",
00044     "E:HMHV2",
00045     "E:HMHV3",
00046     "E:HMHV4",
00047     "E:HMHV5",
00048     "E:HMHV6",
00049     "E:HMHV7",
00050     "E:MM1HV1",
00051     "E:MM1HV2",
00052     "E:MM1HV3",
00053     "E:MM2HV1",
00054     "E:MM2HV2",
00055     "E:MM2HV3",
00056     "E:MM3HV1",
00057     "E:MM3HV2",
00058     "E:MM3HV3",
00059     "E:VP101",
00060     "E:HP101",
00061     "E:HP102",
00062     "E:VP103",
00063     "E:HP104",
00064     "E:HP105",
00065     "E:VP106",
00066     "E:HP107",
00067     "E:VP108",
00068     "E:HP109",
00069     "E:VP110",
00070     "E:VP111",
00071     "E:HP112",
00072     "E:VP113",
00073     "E:HP114",
00074     "E:HP115",
00075     "E:VP116",
00076     "E:HP117",
00077     "E:VP118",
00078     "E:HP119",
00079     "E:HP121",
00080     "E:VP121",
00081     "E:HPTGT",
00082     "E:VPTGT",
00083     "E:M101HM",
00084     "E:M101HS",
00085     "E:M101VM",
00086     "E:M101VS",
00087     "E:M105HM",
00088     "E:M105HS",
00089     "E:M105VM",
00090     "E:M105VS",
00091     "E:M107HM",
00092     "E:M107HS",
00093     "E:M107VM",
00094     "E:M107VS",
00095     "E:M108HM",
00096     "E:M108HS",
00097     "E:M108VM",
00098     "E:M108VS",
00099     "E:M112HM",
00100     "E:M112HS",
00101     "E:M112VM",
00102     "E:M112VS",
00103     "E:M114HM",
00104     "E:M114HS",
00105     "E:M114VM",
00106     "E:M114VS",
00107     "E:M115HM",
00108     "E:M115HS",
00109     "E:M115VM",
00110     "E:M115VS",
00111     "E:M117HM",
00112     "E:M117HS",
00113     "E:M117VM",
00114     "E:M117VS",
00115     "E:M121HM",
00116     "E:M121HS",
00117     "E:M121VM",
00118     "E:M121VS",
00119     "E:MTGTHM",
00120     "E:MTGTHS",
00121     "E:MTGTVM",
00122     "E:MTGTVS",
00123     0
00124 };
00125 const char* profile_monitor[] = {
00126     "E:M101DS",
00127     "E:M105DS",
00128     "E:M107DS",
00129     "E:M108DS",
00130     "E:M112DS",
00131     "E:M114DS",
00132     "E:M115DS",
00133     "E:M117DS",
00134     "E:M121DS",
00135     "E:MTGTDS",
00136     0
00137 };
00138 
00139 const char* hadmu_monitor[] = {
00140     "E:HADMDS",
00141     //    "E:MMA1DS",
00142     //    "E:MMA2DS",
00143     //    "E:MMA3DS",
00144     0
00145 };
00146 
00147 typedef map<string,vector<double> > ScaleMap;
00148 
00149 void init_simple_data(ostream& o)
00150 {
00151     for (int ind=0; simple_data[ind]; ++ind) {
00152         o << simple_data[ind] << " ";
00153     }
00154 }
00155 void init_profile_monitor(ostream& o)
00156 {
00157     for (int ind=0; profile_monitor[ind]; ++ind) {
00158         o << profile_monitor[ind] << "_xmean "
00159           << profile_monitor[ind] << "_xrms "
00160           << profile_monitor[ind] << "_xtot "
00161           << profile_monitor[ind] << "_ymean "
00162           << profile_monitor[ind] << "_yrms "
00163           << profile_monitor[ind] << "_ytot ";
00164     }
00165 }
00166 void init_hadmu_monitor(ostream& o)
00167 {
00168     for (int ind=0; hadmu_monitor[ind]; ++ind) {
00169         o << hadmu_monitor[ind] << "_daesec "
00170           << hadmu_monitor[ind] << "_daemsec "
00171           << hadmu_monitor[ind] << "_vmesec "
00172           << hadmu_monitor[ind] << "_vmemsec "
00173 
00174           << hadmu_monitor[ind] << "_xmean "
00175           << hadmu_monitor[ind] << "_ymean "
00176           << hadmu_monitor[ind] << "_tot "
00177           << hadmu_monitor[ind] << "_high ";
00178 
00179 #ifdef DUMP_HADMU_PIXELS
00180         for (int channel=1; channel <= 49; ++channel) {
00181             o << Form("%s_ch%02d ",hadmu_monitor[ind],channel);
00182         }
00183 #endif
00184     }
00185 }
00186 
00187 void dump_simple_data(ostream& o, const RawBeamMonBlock& rbmb)
00188 {
00189     for (int ind=0; simple_data[ind]; ++ind) {
00190         const RawBeamData* rbd = rbmb[simple_data[ind]];
00191         if (!rbd) o << "0.0 ";
00192         else o << rbd->GetData()[0] << " ";
00193     }
00194 
00195 }
00196 const double MVPERADC = -0.30518;
00197 void dump_one_profile(ostream& o, vector<int>& data, int start,
00198                       const vector<double>& scale)
00199 {
00200     double qx=0,q=0,q2x2=0,max=0;
00201     for (int ind=start; ind < start+44; ++ind) {
00202         double pos = (ind-(start+22.5))*0.5;
00203         double val = MVPERADC*data[ind] - scale[ind];
00204         q += val;
00205         if (val>max) max = val;
00206         double tmp = val*pos;
00207         qx += tmp;
00208         q2x2 += tmp*tmp;
00209     }
00210     if (q == 0.0) o << "0.0 0.0 0.0 ";
00211     else o << qx/q << " " << sqrt(q2x2/(q*q)) << " " << q << " ";
00212     
00213 }
00214 void dump_profile_monitor(ostream& o, const RawBeamMonBlock& rbmb,
00215                           ScaleMap& profmap)
00216 {
00217     for (int ind=0; profile_monitor[ind]; ++ind) {
00218         //cerr << "Getting " << profile_monitor[ind] << endl;
00219         const RawBeamData* rbd = rbmb[profile_monitor[ind]];
00220         assert(rbd);
00221         RawBeamSwicData swic(*rbd);
00222         vector<int> data;
00223         swic.UnscaledWireData(data);
00224         const vector<double>& scale = profmap[profile_monitor[ind]];
00225         dump_one_profile(o,data,2,scale);
00226         dump_one_profile(o,data,50,scale);
00227     }       
00228 
00229 }
00230 void dump_hadmu_monitor(ostream& o, const RawBeamMonBlock& rbmb,
00231                         ScaleMap& pedmap)
00232 {
00233     for (int ind=0; hadmu_monitor[ind]; ++ind) {
00234         //cerr << "Getting " << hadmu_monitor[ind] << endl;
00235         const RawBeamData* rbd = rbmb[hadmu_monitor[ind]];
00236         assert(rbd);
00237         RawBeamSwicData swic(*rbd);
00238         vector<int> data;
00239         swic.UnscaledWireData(data);
00240 
00241         o << rbd->GetSeconds() << " "
00242           << rbd->GetMsecs() << " "
00243           << swic.VmeSeconds() << " "
00244           << swic.VmeNanoseconds()/1000 << " ";
00245 
00246         const vector<double>& scale = pedmap[hadmu_monitor[ind]];
00247         double tot=0,max=0;
00248         double Qx=0, Qy=0, X=0, Y=0;
00249         for (int col=1; col <=7; ++col) {
00250             for (int row=1; row <=7; ++row) {
00251                 int index = hadron_index(hadron_cell(row,col));
00252                 double val = MVPERADC*data[index]-scale[index];
00253                 tot += val;
00254                 if (val>max) max=val;
00255                 Qx += val;
00256                 Qy += val;
00257                 X += val*hadron_position(col);
00258                 Y += val*hadron_position(row);
00259             }
00260         }
00261         if (Qx > 0.0) o << X/Qx << " ";
00262         else o << "0.0 ";
00263         if (Qy > 0.0) o << Y/Qy << " ";
00264         else o << "0.0 ";
00265         o << tot << " " << max << " ";
00266 
00267 #ifdef DUMP_HADMU_PIXELS
00268         for (int channel=1; channel <= 49; ++channel) {
00269             int index = hadron_index(channel);
00270             double val = MVPERADC*data[index]-scale[index];
00271 
00272             o << val << " ";
00273         }
00274 #endif
00275                 
00276     }
00277 }
00278 
00279 
00280 
00281 void load_scale_map(const char* file, ScaleMap &sm, int ndev)
00282 {
00283     ifstream fstr(file);
00284     assert (fstr);
00285     string blah;
00286     vector<string> devs;
00287 
00288     fstr >> blah;
00289     for (int ind=0; ind<ndev; ++ind) {
00290         string dev;
00291         fstr >> dev;
00292         devs.push_back(dev);
00293         cerr << dev << " ";
00294     }
00295     cerr << endl;
00296     for (int ch=0; ch<96; ++ch) {
00297         int channel;
00298         fstr >> channel;
00299         for (int dev=0; dev<ndev; ++dev) {
00300             double scale;
00301             fstr >> scale;
00302             string devname = devs[dev];
00303             //cerr << devname << ": " << scale << endl;
00304             if (!ch) {
00305                 vector<double> v;
00306                 sm[devname] = v;
00307             }
00308             sm[devname].push_back(scale);
00309         }
00310     }
00311 }
00312 
00313 
00314 void bd_text_dump(const char* rootfile,const char* textfile,
00315                   const char* pedfile, const char* swicfile)
00316 {
00317     TFile file(rootfile,"READ");
00318     TTree* tree = (TTree*)(file.Get("BeamMon"));
00319     RawRecord* record = 0;
00320 
00321     ScaleMap pedmap,swicmap;
00322     load_scale_map(pedfile,pedmap,4);
00323     load_scale_map(swicfile,swicmap,10);
00324 
00325     ofstream fstr(textfile);
00326     if (!fstr) {
00327         cerr << "Can't open \"" << textfile << "\" for writing\n";
00328         return;
00329     }
00330     fstr << "# ";
00331     init_simple_data(fstr);
00332     init_profile_monitor(fstr);
00333     init_hadmu_monitor(fstr);
00334     fstr << endl;
00335 
00336     for ( Int_t ient = 0; ient < tree -> GetEntries(); ient++ ) {
00337         tree -> SetBranchAddress("RawRecord",&record);
00338         tree->GetEntry(ient);
00339 
00340         if (!ient) {
00341             const VldContext* vld = record->GetVldContext();
00342             cout << "Reading " << rootfile << ":\n"
00343                  << *vld << endl;
00344         }
00345 
00346         TIter itr = record->GetRawBlockIter();
00347         const RawDataBlock* rdb = 0;
00348 
00349         cerr << "Entry " << ient << endl;
00350 
00351         // loop over blocks in record
00352         while ((rdb = dynamic_cast<RawDataBlock*>(itr()))) {
00353             if (! rdb->InheritsFrom("RawBeamMonBlock")) {
00354                 //cerr << "Doesn't inherit from RawBeamMonBlock" << endl;
00355                 continue;
00356             }
00357             const RawBeamMonBlock* rbmb = 
00358                 dynamic_cast<const RawBeamMonBlock*>(rdb);
00359             assert(rbmb);
00360 
00361             dump_simple_data(fstr,*rbmb);
00362             dump_profile_monitor(fstr,*rbmb,swicmap);
00363             dump_hadmu_monitor(fstr,*rbmb,pedmap);
00364             fstr << endl;
00365             break;
00366         }
00367 
00368         if (tree->GetEntries()-ient == 1) {
00369             const VldContext* vld = record->GetVldContext();
00370             cout << "finishing " << rootfile << ":\n"
00371                  << *vld << endl;
00372         }
00373         delete record; record = 0;
00374     }
00375     
00376 }