loon: CheckGC.cxx Source File

00001 
00002 // $Id: CheckGC.cxx,v 1.6 2006/04/19 19:59:06 rhatcher Exp $
00003 //
00004 // read GC data from database and make gc_summary ntuples
00005 //
00006 // S. Murgia and T. Yang
00008 #include "CheckGC.h"
00009 
00010 #include "PulserCalibration/PulserXScale.h"
00011 #include "PulserCalibration/PulserGain.h"
00012 #include "PulserCalibration/PulserGainPin.h"
00013 
00014 #include "TH2F.h"
00015 #include "TF1.h"
00016 #include "TSystem.h"
00017 #include "TTree.h"
00018 #include "TFile.h"
00019 
00020 #include "MessageService/MsgService.h"
00021 #include "Validity/VldContext.h"
00022 #include "Validity/VldRange.h"
00023 #include "Plex/PlexPinDiodeId.h"
00024 #include "Plex/PlexHandle.h"
00025 #include "RawData/RawChannelId.h"
00026 #include "DatabaseInterface/DbiWriter.h"
00027 #include "DatabaseInterface/DbiResultPtr.h"
00028 #include "DatabaseInterface/DbiTableProxy.h"
00029 #include "DatabaseInterface/DbiDBProxy.h"
00030 #include "DatabaseInterface/DbiCache.h"
00031 #include "DatabaseInterface/DbiSqlContext.h"
00032 #include "DatabaseInterface/DbiValidityRec.h"
00033 #include "HistMan/HistMan.h"
00034 
00035 ClassImp(CheckGC)
00036 
00037 CVSID("$Id: CheckGC.cxx,v 1.6 2006/04/19 19:59:06 rhatcher Exp $");
00038 
00039 //......................................................................
00040 
00041 CheckGC::CheckGC()
00042 {
00043   gc = new GCSummary();
00044 }
00045 
00046 //......................................................................
00047 
00048 CheckGC::~CheckGC()
00049 {
00050 }
00051 
00052 char* digit(UInt_t num){
00053   char* str = new char[2];
00054   if (num<10){
00055     sprintf(str,"0%d",num);
00056     return str;
00057   }
00058   else if (num>=10&&num<100){
00059     sprintf(str,"%d",num);
00060     return str;
00061   }
00062   return str;
00063 }
00064 
00065 void FillArray(float* array1, const Float_t* array2, int numpoints)
00066 {
00067   for (int i = 0; i< numpoints; i++){
00068     array1[i] = array2[i];
00069   }
00070 }
00071 
00072 void CheckGC::SettsStart(UInt_t year, UInt_t month,
00073                     UInt_t day,  UInt_t hour,
00074                     UInt_t min,  UInt_t sec,
00075                     UInt_t nsec, 
00076                     Bool_t isUTC, Int_t secOffset)
00077 {
00078   VldTimeStamp ts(year, month, day, hour, min, sec, nsec, isUTC, secOffset);
00079   ts.Copy(tsStart);
00080 }
00081 
00082 void CheckGC::SettsEnd(UInt_t year, UInt_t month,
00083                     UInt_t day,  UInt_t hour,
00084                     UInt_t min,  UInt_t sec,
00085                     UInt_t nsec, 
00086                     Bool_t isUTC, Int_t secOffset)
00087 {
00088   VldTimeStamp ts(year, month, day, hour, min, sec, nsec, isUTC, secOffset);
00089   ts.Copy(tsEnd);
00090 }
00091 
00092 void CheckGC::SetDetType(Int_t dettype) 
00093 {
00094   if (dettype==1) fDetector = Detector::kNear;
00095   else if (dettype==2) fDetector = Detector::kFar;
00096   else fDetector = Detector::kUnknown;
00097   //fDetector = dettype;
00098 }
00099 
00100 void CheckGC::checklin()
00101 {
00102 
00103   //create the new TFile for holding the GC tree
00104   char filename[100];
00105 
00106   cout<<"tsStart: "<<tsStart<<endl;
00107   cout<<"tsEnd:   "<<tsEnd<<endl;
00108 
00109   UInt_t year1,month1,day1;
00110   UInt_t year2,month2,day2;
00111   tsStart.GetDate(true,0,&year1,&month1,&day1);
00112   tsEnd.GetDate(true,0,&year2,&month2,&day2);
00113 
00114   if (fDetector == Detector::kNear){
00115     cout<<"Checking near detector gain curves"<<endl;
00116     sprintf(filename,"ngc%s%s%s.root",digit(year1%100),digit(month1),digit(day1));
00117   }
00118   else if (fDetector == Detector::kFar){
00119     cout<<"Checking far detector gain curves"<<endl;
00120     sprintf(filename,"fgc%s%s%s.root",digit(year1%100),digit(month1),digit(day1));
00121   }
00122 
00123   TFile *fNtpFile = new  TFile(filename,"RECREATE");
00124   fNtpFile->cd();
00125   cout<<"File opened "<<filename<<endl;
00126   TTree *fNtuple = new TTree("gc_tree","GC Analysis Tree");
00127   fNtuple->Branch("","GCSummary",&gc,64000,1);
00128 
00129   //Set validity context
00130   VldContext vldc(fDetector,SimFlag::kData,tsStart);
00131   PlexHandle plex(vldc);
00132 
00133   DbiSqlContext extContext(DbiSqlContext::kStarts,tsStart,tsEnd,fDetector,SimFlag::kData);
00134 
00135   string extClause = "stripend=124017 || stripend=124016";
00136 
00137   // perform the query
00138   DbiResultPtr<PulserGainPin> rsPin("PULSERGAINPIN",extContext,Dbi::kAnyTask); 
00139   DbiResultPtr<PulserGain> rsStripEnd("PULSERGAIN",extContext,Dbi::kAnyTask);
00140   //DbiResultPtr<PulserGain> rsStripEnd("PULSERGAIN",extContext,Dbi::kAnyTask,extClause.c_str()); 
00141   
00142   const int numPin = rsPin.GetNumRows();
00143   const int numStripEnd = rsStripEnd.GetNumRows();
00144   
00145   cout << "got " << numPin << " PINs" << endl;
00146   cout << "got " << numStripEnd << " StripEnds"<<endl;
00147   map<int,PulserGainPin*> hgpinmap;
00148   map<int,PulserGainPin*> lgpinmap;
00149   map<int,VldTimeStamp> pintimeStart; //Currently PulserGain and PulserGainPin share the same timestamp for entries with the same aggregateno. We'll record timestamps to synchronize CALPULSERFITS with PulserGain and PulserGainPin.
00150   map<int,VldTimeStamp> pintimeEnd;
00151 
00152   //  HistMan *hm = new HistMan("allres");
00153 
00154   for (int ir=0; ir<numPin; ir++) {
00155     
00156     //Pointer to row
00157     const PulserGainPin* rs = rsPin.GetRow(ir);
00158     
00159     int pingain = rs->GetPinDiodeId().GetGain();
00160     int aggregateno = rs->GetAggregateNo();
00161     PlexLedId led(aggregateno);
00162 
00163     // Get the DbiValidityRec for the specfic row and hence its start and end time.
00164     // Unfortunately, because of overlaps, it may not be exactly what you
00165     // want.  The DBI is interested in effective validities i.e. for what
00166     // period of time an entry is the highest priority.  Normally a new
00167     // calibration will be made before the old one has expired and the two
00168     // will overlap with the later entry taking priority.  So in the normal
00169     // situation, timeStart as derived above will be O.K., but the timeEnd
00170     // is likely to be too early.  --Nick West
00171     // The following code was suggested by Nick
00172     const DbiValidityRec* vrec  = rsPin.GetValidityRec(rs);
00173     UInt_t seqno                = vrec->GetSeqNo();
00174     DbiResultPtr<PulserGainPin> rsPin2("PULSERGAINPIN",seqno,0);  //0 = cascade entry
00175     const DbiValidityRec* vrec2 = rsPin2.GetValidityRec();
00176     const VldRange& vrng        = vrec2->GetVldRange();
00177     VldTimeStamp timeStart      = vrng.GetTimeStart();
00178     VldTimeStamp timeEnd        = vrng.GetTimeEnd();
00179     //insert timeStart and timeEnd
00180     pintimeStart.insert(pair<int,VldTimeStamp>(aggregateno,timeStart));
00181     pintimeEnd.insert(pair<int,VldTimeStamp>(aggregateno,timeEnd));
00182 
00183     //If low gain (gain = 1) save first
00184     if (pingain == 1) {
00185       lgpinmap.insert(pair<int,PulserGainPin*>(aggregateno,const_cast<PulserGainPin*>(rs) ));
00186     }
00187     //If high gain (gain = 0)
00188     if (fDetector==Detector::kNear && led.GetPulserBox()==0) continue;//no high gain pin diode chennles in the near detector
00189     //if (fDetector==1 && led.GetPulserBox()==0) continue;//no high gain pin diode chennles in the near detector
00190     if (pingain == 0) { 
00191       hgpinmap.insert(pair<int,PulserGainPin*>(aggregateno,const_cast<PulserGainPin*>(rs) ));
00192     }
00193   }
00194 
00195   cout<<"Finish reading PULSERGAINPIN"<<endl;
00196 
00197   int counter = 0;
00198   int countled = 0;
00199   int aggregateno_old = -1;
00200 
00201   gc->dettype = int(fDetector);
00202 
00203   for (int ir=0; ir<numStripEnd; ir++) {
00204     if (ir%1000==0) cout<<ir<<"/"<<numStripEnd<<endl;
00205     
00206     gc->Reset();
00207     
00208     //Pointer to row
00209     const PulserGain* rs = rsStripEnd.GetRow(ir);
00210     int aggregateno = rs->GetAggregateNo();
00211     int stripendid = rs->GetStripEnd();
00212     PlexLedId ledid0(aggregateno);
00213     PlexLedId ledid(fDetector,ledid0.GetPulserBox(),ledid0.GetLedInBox());
00214     gc->aggno = aggregateno;
00215     gc->led = ledid.GetLedInBox();
00216     gc->pb = ledid.GetPulserBox();
00217     gc->stripend = stripendid;
00218     (pintimeStart.find(aggregateno)->second).Copy(gc->timeStart);
00219     (pintimeEnd.find(aggregateno)->second).Copy(gc->timeEnd);
00220     // Get mean and entries to calculate residual
00221     int numpoints = rs->GetNumPoints();
00222     
00223     gc->numPointsPmt = rs->GetNumPoints();
00224     //if (ledid.GetPulserBox()!=0) continue;
00225     
00226     if (numpoints<20) {
00227       cout << "Less than 20 points in GC! Number of entries: " << numpoints << endl;
00228     }
00229 
00230     FillArray(gc->meanPmtraw,rs->GetMean(),gc->numPointsPmt);
00231     FillArray(gc->errorPmt,rs->GetError(),gc->numPointsPmt);
00232     FillArray(gc->numEntriesPmt,rs->GetNumEntries(),gc->numPointsPmt);
00233     FillArray(gc->numTrigsPmt,rs->GetNumTriggers(),gc->numPointsPmt);
00234     
00235     if ( aggregateno_old != aggregateno) {
00236       
00237       aggregateno_old = aggregateno;
00238       counter = 0;
00239       countled++;
00240     } else {
00241       //if (counter>5) continue;
00242       
00243     }
00244     
00245     PlexStripEndId seid = PlexStripEndId::UnbuildPlnStripEndKey(fDetector,stripendid);
00246     int planeno = seid.GetPlane();
00247     int stripno = seid.GetStrip();
00248     
00249     gc->plane = planeno;
00250     gc->strip = stripno;
00251 
00252 
00253     // Get validity of row
00254     const DbiValidityRec* dbival =  rsStripEnd.GetValidityRec(rs);
00255     
00256     // near end: task = 1
00257     // far end: task = 2
00258     int task = dbival->GetTask();
00259     if (task==1) {
00260       gc->nearfar = 0;
00261       if (lgpinmap.find(aggregateno)!=lgpinmap.end()){//low gain
00262         gc->highlow = 1;
00263         gc->fnhlcode = gc->nearfar*2 + gc->highlow;
00264         gc->numPointsPin = (lgpinmap.find(aggregateno)->second)->GetNumPoints();
00265         FillArray(gc->meanPinraw,(lgpinmap.find(aggregateno)->second)->GetMean(),gc->numPointsPin);
00266         FillArray(gc->errorPin,(lgpinmap.find(aggregateno)->second)->GetError(),gc->numPointsPin);
00267         FillArray(gc->numEntriesPin,(lgpinmap.find(aggregateno)->second)->GetNumEntries(),gc->numPointsPin);
00268         FillArray(gc->numTrigsPin,(lgpinmap.find(aggregateno)->second)->GetNumTriggers(),gc->numPointsPin);
00269         this->ZeroCorr();
00270         fNtuple->Fill();
00271       }
00272       if (hgpinmap.find(aggregateno)!=hgpinmap.end()){//high gain
00273         gc->highlow = 0;
00274         gc->fnhlcode = gc->nearfar*2 + gc->highlow;
00275         gc->numPointsPin = (hgpinmap.find(aggregateno)->second)->GetNumPoints();
00276         FillArray(gc->meanPinraw,(hgpinmap.find(aggregateno)->second)->GetMean(),gc->numPointsPin);
00277         FillArray(gc->errorPin,(hgpinmap.find(aggregateno)->second)->GetError(),gc->numPointsPin);
00278         FillArray(gc->numEntriesPin,(hgpinmap.find(aggregateno)->second)->GetNumEntries(),gc->numPointsPin);
00279         FillArray(gc->numTrigsPin,(hgpinmap.find(aggregateno)->second)->GetNumTriggers(),gc->numPointsPin);
00280         this->ZeroCorr();
00281         fNtuple->Fill();
00282       }
00283     }
00284 
00285     else if (task==2) {
00286       gc->nearfar = 1;
00287       if (lgpinmap.find(aggregateno)!=lgpinmap.end()){//low gain
00288         gc->highlow = 1;
00289         gc->fnhlcode = gc->nearfar*2 + gc->highlow;
00290         gc->numPointsPin = (lgpinmap.find(aggregateno)->second)->GetNumPoints();
00291         FillArray(gc->meanPinraw,(lgpinmap.find(aggregateno)->second)->GetMean(),gc->numPointsPin);
00292         FillArray(gc->errorPin,(lgpinmap.find(aggregateno)->second)->GetError(),gc->numPointsPin);
00293         FillArray(gc->numEntriesPin,(lgpinmap.find(aggregateno)->second)->GetNumEntries(),gc->numPointsPin);
00294         FillArray(gc->numTrigsPin,(lgpinmap.find(aggregateno)->second)->GetNumTriggers(),gc->numPointsPin);
00295         this->ZeroCorr();
00296         fNtuple->Fill();
00297       }
00298       if (hgpinmap.find(aggregateno)!=hgpinmap.end()){//high gain
00299         gc->highlow = 0;
00300         gc->fnhlcode = gc->nearfar*2 + gc->highlow;
00301         gc->numPointsPin = (hgpinmap.find(aggregateno)->second)->GetNumPoints();
00302         FillArray(gc->meanPinraw,(hgpinmap.find(aggregateno)->second)->GetMean(),gc->numPointsPin);
00303         FillArray(gc->errorPin,(hgpinmap.find(aggregateno)->second)->GetError(),gc->numPointsPin);
00304         FillArray(gc->numEntriesPin,(hgpinmap.find(aggregateno)->second)->GetNumEntries(),gc->numPointsPin);
00305         FillArray(gc->numTrigsPin,(hgpinmap.find(aggregateno)->second)->GetNumTriggers(),gc->numPointsPin);
00306         this->ZeroCorr();
00307         fNtuple->Fill();
00308       }
00309     }
00310   }
00311   
00312   fNtpFile->Write();
00313   fNtpFile->Close();
00314 }
00315 
00316 void CheckGC::ZeroCorr()
00317 {
00318   for (int i = 0; i<gc->numPointsPin; i++){
00319     if (gc->numTrigsPin[i]) 
00320       gc->meanPin[i] = gc->meanPinraw[i]*gc->numEntriesPin[i]/gc->numTrigsPin[i];
00321     else gc->meanPin[i] = gc->meanPinraw[i];
00322   }
00323 
00324   for (int i = 0; i<gc->numPointsPmt; i++){
00325     if (gc->numTrigsPmt[i])
00326       gc->meanPmt[i] = gc->meanPmtraw[i]*gc->numEntriesPmt[i]/gc->numTrigsPmt[i];
00327     else gc->meanPmt[i] = gc->meanPmtraw[i];
00328   }
00329 }