loon: TimeCalibratorSRModule.cxx Source File

00001 
00002 // $Id: TimeCalibratorSRModule.cxx,v 1.8 2007/11/11 04:53:36 rhatcher Exp $
00003 //
00004 // TimeCalibratorSRModule.cxx
00005 //
00006 // Begin_Html<img src="../../pedestrians.gif" align=center>
00007 // <a href="../source_warning.html">Warning for beginners</a>.<br> 
00008 //
00009 // A JobControl Module for calculating T0 time calibration constants.
00010 //
00011 // Author:  R. Lee 2002.03.01
00012 //
00014 
00015 #include <cassert>
00016 #include <cmath>
00017 #include <fstream>
00018 #include <iostream>
00019 
00020 
00021 #include "Algorithm/AlgConfig.h"
00022 #include "Algorithm/AlgFactory.h"
00023 #include "Algorithm/AlgHandle.h"
00024 #include "CandData/CandHeader.h"
00025 #include "CandData/CandRecord.h"
00026 #include "CandDigit/CandDigitListHandle.h"
00027 #include "CandDigit/CandDigitHandle.h"
00028 #include "CandTrackSR/CandTrackSRListHandle.h"
00029 #include "CandTrackSR/CandTrackSRHandle.h"
00030 #include "CandTrackSR/TrackSRListModule.h"
00031 #include "Candidate/CandContext.h"
00032 #include "Conventions/Munits.h"
00033 #include "Conventions/StripEnd.h"
00034 #include "JobControl/JobCModuleRegistry.h"
00035 #include "JobControl/JobCommand.h"
00036 #include "MessageService/MsgService.h"
00037 #include "MinosObjectMap/MomNavigator.h"
00038 #include "Plex/PlexStripEndId.h"
00039 #include "RawData/RawDigit.h"
00040 #include "RawData/RawHeader.h"
00041 #include "RawData/RawRecord.h"
00042 #include "RawData/RawChannelId.h"
00043 #include "RawData/RawDaqSnarlHeader.h"
00044 #include "RawData/RawDaqHeaderBlock.h"
00045 #include "RawData/RawDigitDataBlock.h"
00046 #include "RawData/RawVarcErrorInTfBlock.h"
00047 #include "RecoBase/CandSliceListHandle.h"
00048 #include "RecoBase/CandStripHandle.h"
00049 #include "RecoBase/CandStripListHandle.h"
00050 #include "RecoBase/CandTrackHandle.h"
00051 #include "RecoBase/CandTrackListHandle.h"
00052 #include "TimeCalibratorSR/TimeCalibratorSRModule.h"
00053 #include "UgliGeometry/UgliGeomHandle.h"
00054 #include "UgliGeometry/UgliStripHandle.h"
00055 #include "Validity/VldContext.h"
00056 
00057 ClassImp(TimeCalibratorSRModule)
00058 
00059 //......................................................................
00060 CVSID("$Id: TimeCalibratorSRModule.cxx,v 1.8 2007/11/11 04:53:36 rhatcher Exp $");
00061 JOBMODULE(TimeCalibratorSRModule, "TimeCalibratorSRModule",
00062          "Builds CandTimeCalibratorSRList from CandSliceList");
00063 
00064 //......................................................................
00065 TimeCalibratorSRModule::TimeCalibratorSRModule()
00066 {
00067   MSG("TimeCalibratorSR", Msg::kVerbose) << "TimeCalibratorSRModule::Constructor\n";
00068 
00069   vec_dtime = new vector<Float_t>;
00070   vec_weight = new vector<Float_t>;
00071   vec_vachip = new vector<Short_t>;
00072   vec_vachip0 = new vector<Short_t>;
00073   vec_plnstripendkey = new vector<Int_t>;
00074 
00075   Int_t nsize = 4000000; // about one day's worth of data
00076 
00077   vec_dtime->reserve(nsize);
00078   vec_weight->reserve(nsize);
00079   vec_vachip->reserve(nsize);
00080   vec_vachip0->reserve(nsize);
00081   vec_plnstripendkey->reserve(nsize);
00082 
00083   lastplane=0;
00084 
00085 }
00086 
00087 //......................................................................
00088 TimeCalibratorSRModule::~TimeCalibratorSRModule() 
00089 {
00090   MSG("TimeCalibratorSR", Msg::kVerbose) << "TimeCalibratorSRModule::Destructor\n";
00091 
00092   if (vec_dtime) delete vec_dtime;
00093   if (vec_weight) delete vec_weight;
00094   if (vec_vachip) delete vec_vachip;
00095   if (vec_vachip0) delete vec_vachip0;
00096   if (vec_plnstripendkey) delete vec_plnstripendkey;
00097 }
00098 
00099 //......................................................................
00100 void TimeCalibratorSRModule::BeginJob() 
00101 {
00102   MSG("TimeCalibratorSR", Msg::kVerbose) << "TimeCalibratorSRModule::BeginJob\n";
00103 
00104   for (Int_t i=0; i<1728; i++) {
00105     for (Int_t j=0; j<1728; j++) {
00106       wsum[i][j] = 0.;
00107       delt[i][j] = 0.;
00108     }
00109   }
00110 
00111   for (Int_t i=0; i<384*512; i++) {
00112     timeoffset[i] = 0.;
00113   }
00114 
00115   fFile = new TFile("timecalibratorsr.root","RECREATE");
00116 
00117   th2_dtime = new TH2F("t0vsvachip","T0 vs VAChip",1728,0.,1728.,100,-50.,50.);
00118   th2_dtime->SetXTitle("VA Chip");
00119   th2_dtime->SetYTitle("T0  (ns)");
00120 
00121   th2_dtimeend = new TH2F("dtvsvachip","#Delta T vs VAChip",1728,0.,1728.,100,-50.,50.);
00122   th2_dtimeend->SetXTitle("VA Chip");
00123   th2_dtimeend->SetYTitle("#Delta T  (ns)");
00124 
00125 }
00126 
00127 JobCResult TimeCalibratorSRModule::Ana(const MomNavigator *mom)
00128 {
00129 
00130   JobCResult result(JobCResult::kPassed);
00131 
00132 
00133   MSG("TimeCalibratorSR", Msg::kVerbose) << "TimeCalibratorSRModule::Ana\n";
00134 
00135 
00136 
00137 
00138   CandRecord* candrec = dynamic_cast<CandRecord*>
00139     (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
00140   if (candrec==0) {
00141     MSG("TimeCalibratorSR", Msg::kWarning)
00142       << "No PrimaryCandidateRecord in MOM." << endl;
00143     result.SetWarning().SetFailed();
00144     return result;
00145   }
00146 
00147   CandTrackListHandle *tracklist = dynamic_cast<CandTrackListHandle*>
00148     (candrec->FindCandHandle("CandTrackListHandle"));
00149   if (!tracklist) {
00150     MSG("TimeCalibratorSR", Msg::kVerbose)
00151       << "No CandTrackListHandle in CandRecord." << endl;
00152     result.SetFailed();
00153     return result;
00154   }
00155 
00156 
00157 
00158   RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
00159   if (rr == 0) {
00160     MSG("TimeCalibratorSR", Msg::kWarning) << "No RawRecord in MOM." << endl;
00161     result.SetWarning().SetFailed();
00162     return result;
00163   }
00164 
00165   const RawDaqSnarlHeader* snarlHdr =
00166      dynamic_cast<const RawDaqSnarlHeader*>(rr->GetRawHeader());
00167   if (snarlHdr) {
00168   }
00169 
00170 
00171   const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *>
00172                         (rr->FindRawBlock("RawDigitDataBlock"));
00173   if (rddb == 0) {
00174     MSG("TimeCalibratorSR", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
00175     result.SetWarning().SetFailed();
00176     return result;
00177 
00178   }
00179 
00180   TIter trackItr(tracklist->GetDaughterIterator());
00181 
00182   CandTrackHandle *track = dynamic_cast<CandTrackHandle*>(trackItr());
00183 
00184   if (!track) {
00185     MSG("TimeCalibratorSR", Msg::kVerbose) << "No tracks in CandTrackListHandle." << endl;
00186     result.SetFailed();
00187     return result;
00188   }
00189 
00190   CandTrackSRHandle *tracksr = 0;
00191   if (track->InheritsFrom("CandTrackSRHandle")) {
00192     tracksr = dynamic_cast<CandTrackSRHandle *>(track);
00193   }
00194   Double_t zdircos = fabs(track->GetDirCosZ());
00195   if (tracksr && tracksr->GetHoughResid2()>=0.) {
00196     zdircos = fabs(tracksr->GetHoughDirCosZ());
00197   }
00198   if (zdircos<=0.) {
00199     MSG("TimeCalibratorSR",Msg::kWarning) << "Z Direction Cosine <=0, setting to 1" << endl;
00200     zdircos = 1.;
00201   } 
00202 
00203 
00204 
00205   TIter stripItr(track->GetDaughterIterator());
00206 
00207 // find the first vachip
00208   Int_t ivachip0 = -1;
00209   CandDigitHandle *digit0 = 0;
00210   Double_t time0 = -1;
00211   Double_t weight0 = 0.;
00212   Double_t z0 = 0.;
00213   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
00214     TIter digitItr(strip->GetDaughterIterator());
00215     while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())) {
00216       const RawDigit *rd = rddb->At(digit->GetRawDigitIndex());
00217       if (rd) {
00218         RawChannelId rawid = rd->GetChannel();
00219         Int_t ivachip = rawid.GetCrate()*108+rawid.GetVarcId()*36+rawid.GetVmm()*6+rawid.GetVaAdcSel()*3+rawid.GetVaChip();
00220         if (GoodDigit(track,digit) && (ivachip0<0 || ivachip<ivachip0)) {
00221           ivachip0 = ivachip;
00222           digit0 = digit;
00223           time0 = track->GetT(strip->GetPlane(),digit->GetPlexSEIdAltL().GetBestSEId().GetEnd());
00224           weight0 = 1.;
00225           z0 = strip->GetZPos();
00226         }
00227       }
00228     }
00229   }
00230 
00231 
00232   stripItr.Reset();
00233   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
00234     TIter digitItr(strip->GetDaughterIterator());
00235     Float_t dtimeend=0.;
00236     Int_t ndigit=0;
00237     Int_t ivachip = -1;
00238     Int_t thisplane = strip->GetPlane();
00239     if (thisplane>lastplane) {
00240       lastplane = thisplane;
00241     }
00242     while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())) {
00243       const RawDigit *rd = rddb->At(digit->GetRawDigitIndex());
00244       if (rd) {
00245         RawChannelId rawid = rd->GetChannel();
00246         ivachip = rawid.GetCrate()*108+rawid.GetVarcId()*36+rawid.GetVmm()*6+rawid.GetVaAdcSel()*3+rawid.GetVaChip();
00247         Float_t timesign = 1.;
00248         if (track->GetTimeSlope()<0.) timesign = -1.;
00249         Float_t corrtime = track->GetT(strip->GetPlane(),digit->GetPlexSEIdAltL().GetBestSEId().GetEnd());
00250         Float_t dtime0 = corrtime-timesign*(strip->GetZPos()-z0)/zdircos/(3.e8)-time0;
00251         if (digit->GetPlexSEIdAltL().GetBestSEId().GetEnd()==StripEnd::kNegative) {
00252           dtimeend -= dtime0;
00253           ndigit++;
00254         }
00255         else if (digit->GetPlexSEIdAltL().GetBestSEId().GetEnd()==StripEnd::kPositive) {
00256           dtimeend += dtime0;
00257           ndigit++;
00258         }
00259         if (ivachip!=ivachip0 && GoodDigit(track,digit)) {
00260           vec_dtime->push_back(dtime0);
00261           vec_weight->push_back(weight0);
00262           vec_vachip->push_back(ivachip);
00263           vec_vachip0->push_back(ivachip0);
00264           vec_plnstripendkey->push_back(digit->GetPlexSEIdAltL().GetBestSEId().BuildPlnStripEndKey());
00265           delt[ivachip][ivachip0] += weight0*dtime0;
00266           wsum[ivachip][ivachip0] += weight0;
00267         }
00268       }
00269     }
00270     dtimeend *= 1.e9;
00271     if (ndigit==2) {
00272       th2_dtimeend->Fill(float(ivachip),dtimeend);
00273     }
00274   }
00275 
00276   return result;
00277 }
00278 
00279 //......................................................................
00280 Bool_t TimeCalibratorSRModule::GoodDigit(CandTrackHandle * /* track */,
00281                                          CandDigitHandle *digit) 
00282 {
00283 // return true if this strip should be used in timing calibration
00284   if (digit->GetCharge()>200.) {
00285     return kTRUE;
00286   }
00287   return kFALSE;
00288 }
00289 
00290 //......................................................................
00291 void TimeCalibratorSRModule::HandleCommand(JobCommand *command) 
00292 {
00293 
00294 // Process configuration commands
00295   MSG("TimeCalibratorSR", Msg::kDebug) << "TimeCalibratorSRModule::HandleCommand" << endl;
00296 
00297   TString cmd = command->PopCmd();
00298   if (cmd == "Set") {
00299     TString opt = command->PopOpt();
00300     MSG("TimeCalibratorSR", Msg::kWarning)
00301            << "TimeCalibratorSRModule::HandleCommand: Unrecognized option "
00302                                                        << opt << endl;
00303   } else {
00304       MSG("TimeCalibratorSR", Msg::kWarning)
00305             << "TimeCalibratorSRModule::HandleCommand: Unrecognized command "
00306                                                          << cmd << endl;
00307   }
00308 }
00309 
00310 void TimeCalibratorSRModule::EndJob()
00311 {
00312 //======================================================================
00313 // At the end of the job print some stuff out and save the histogram
00314 // to a file
00315 //======================================================================
00316 
00317 
00318 
00319   for (int ivachip=1; ivachip<1728; ivachip++) {
00320     Double_t t0 = 0.;
00321     Double_t w0 = 0.;
00322     for (int ivachip1=ivachip-1; ivachip1>=0; ivachip1--) {
00323       if (ivachip1==0) {
00324         t0 += delt[ivachip][ivachip1];
00325         w0 += wsum[ivachip][ivachip1];
00326       }
00327       else if (wsum[ivachip][ivachip1]>0. && wsum[ivachip1][0]>0.) {
00328         t0 += (delt[ivachip][ivachip1]/wsum[ivachip][ivachip1]+delt[ivachip1][0])/(1./wsum[ivachip][ivachip1]+1./wsum[ivachip1][0]);
00329         w0 += 1./(1./wsum[ivachip][ivachip1]+1./wsum[ivachip1][0]);
00330       }
00331     }
00332     if (w0>0.) {
00333       t0 /= w0;
00334       delt[ivachip][0] = t0;
00335       wsum[ivachip][0] = w0;
00336     }
00337     MSG("TimeCalibratorSR",Msg::kVerbose) << "T0 " << ivachip << " " << delt[ivachip][0]*1.e9 << " " << wsum[ivachip][0] << "\n";
00338     char histname[80];
00339     sprintf(histname,"t0_%04d",ivachip);
00340   }
00341 
00342 
00343   vector<Float_t>::iterator iter_dtime;
00344   vector<Float_t>::iterator iter_weight;
00345   vector<Short_t>::iterator iter_vachip;
00346   vector<Short_t>::iterator iter_vachip0;
00347   vector<Int_t>::iterator iter_plnstripendkey;
00348   for (iter_dtime = vec_dtime->begin(), iter_weight = vec_weight->begin(), iter_vachip = vec_vachip->begin(), iter_vachip0 = vec_vachip0->begin(), iter_plnstripendkey = vec_plnstripendkey->begin(); iter_dtime!= vec_dtime->end();) {
00349     Int_t ivachip = (Int_t)(*iter_vachip);
00350     Int_t ivachip0 = (Int_t)(*iter_vachip0);
00351     Int_t plnstripendkey = (Int_t)(*iter_plnstripendkey);
00352     Float_t weight = *iter_weight;
00353     Float_t dtime = *iter_dtime;
00354     Float_t yhist = (dtime+delt[ivachip0][0])*1.e9;
00355     Float_t xhist = (Float_t)ivachip;
00356     Float_t whist = 0.;
00357     if (ivachip0==0) {
00358       whist = weight;
00359     }
00360     else if (weight>0. && wsum[ivachip0][0]>0.) {
00361       whist = 1./(1./weight+1./wsum[ivachip0][0]);
00362     }
00363     th2_dtime->Fill(xhist,yhist,whist);
00364     if (plnstripendkey>=0 && plnstripendkey<512*384) {
00365       timeoffset[plnstripendkey] = delt[ivachip][0];
00366     }
00367     iter_dtime++;
00368     iter_weight++;
00369     iter_vachip++;
00370     iter_vachip0++;
00371     iter_plnstripendkey++;
00372   }
00373 
00374   UpdateTable();
00375 
00376   th2_dtime->Write();
00377   th2_dtimeend->Write();
00378 
00379 
00380   fFile->Write();
00381   fFile->Close();
00382 
00383 }
00384 
00385 
00386 void TimeCalibratorSRModule::UpdateTable()
00387 {
00388 
00389   Int_t SeqNo = 200000000;  // for far detector
00390   Detector::Detector_t det = Detector::kFar; // or kNear
00391   SimFlag::SimFlag_t simflag = SimFlag::kData;  // or kReroot
00392   VldTimeStamp time(2009,2,1,0,0,0);
00393   char timestamp1[21];
00394   sprintf(timestamp1,"'2000-01-01 00:00:00'");
00395   char timestamp2[21];
00396   sprintf(timestamp2,"'2010-01-01 00:00:00'");
00397   char timestamp3[21];
00398   sprintf(timestamp3,"'2001-11-07 00:00:00'");
00399   
00400   int mode=0;
00401   int aggnum=-1;
00402 
00403 
00404 /* time */
00405 
00406   ofstream TimeCalibrationTable("TimeCalibrationTable.far");
00407   TimeCalibrationTable << "# SeqNo  \t SEIDkey \t SEID   Chi2 \t Scale \t Offset  "  << endl; 
00408 
00409   ofstream TimeCalibrationValidity("TimeCalibrationValidity.far");
00410   TimeCalibrationValidity << "# SeqNo \t TimeStart \t TimeEnd \t DetMask \t SimMask \t Mode \t AggrNo \t CreationTime \t InsertTime" << endl;
00411 
00412 
00413   TimeCalibrationValidity << SeqNo << '\t' << timestamp1 << '\t' << timestamp2 << '\t' <<(Int_t) det << '\t' << (Int_t) simflag << '\t' << mode << '\t' << aggnum << '\t' <<  timestamp3 << '\t' << timestamp3 << endl;
00414 
00415 
00416 
00417 
00418   cout << " UgliLoanPool::SetAlwaysUseDbi" << endl;
00419  //  UgliLoanPool::SetAlwaysUseDbi(true);
00420  //  SimFlag::SimFlag_t simflag = SimFlag::kReroot;
00421  //  VldTimeStamp time(2001,11,1,0,0,0);
00422 
00423    // if using kReroot for kFar avoid using 2001-09-18 to 2001-10-31
00424    // for timestamp except for "compressed detector"
00425    //
00426    // if using kData for kFar a linear schedule of 1 steel plane per
00427    // day is assumed starting with plane 0 on 2001-08-01 and ending
00428    // with steel plane 485 on 2002-11-29
00429 
00430    // create the appropriate context
00431    VldContext vldc(det,simflag,time);
00432 
00433    // decide now whether we have two-sided or one-sided readout
00434    StripEnd::StripEnd_t one_sided[] = { StripEnd::kWest };
00435    StripEnd::StripEnd_t two_sided[] = { StripEnd::kEast,
00436                                         StripEnd::kWest };
00437    StripEnd::StripEnd_t *side_list = two_sided;
00438    int nsides = 2;
00439    if (det==Detector::kNear) { side_list = one_sided; nsides = 1; }
00440    UgliGeomHandle ugh(vldc);
00441 
00442 
00443    for (int ipln=0; ipln <= lastplane; ipln++) {   // note: "<=" not "<"
00444       PlexPlaneId plnid(det,ipln,false);
00445       // skip planes that don't have an strips
00446       PlaneCoverage::PlaneCoverage_t cover = plnid.GetPlaneCoverage();
00447       if (cover == PlaneCoverage::kNoActive ||
00448           cover == PlaneCoverage::kUnknown)    continue;
00449 
00450 
00451       // get a handle for the plane and ask it how many strips it has
00452          UgliScintPlnHandle usph = ugh.GetScintPlnHandle(plnid);
00453          int nstrips = usph.NumberOfStrips();
00454          for (int istrip=0; istrip < nstrips; istrip++) {
00455            
00456            
00457            for (int indx_side=0; indx_side<nsides; indx_side++) {
00458              
00459              PlexStripEndId oneend(plnid,istrip,side_list[indx_side]);
00460              
00461              // what the heck, check that the geometry isn't screwed up
00462              UgliStripHandle ush = ugh.GetStripHandle(oneend);
00463              if (!ush.IsValid()) cout << "invalid handle?" << endl;
00464              
00465              Int_t SEID = oneend.GetEncoded();
00466              
00467              Int_t SEIDkey = oneend.BuildPlnStripEndKey();
00468              
00469 
00470              TimeCalibrationTable << SeqNo  << '\t' << SEIDkey  << '\t' << SEID << '\t';  
00471              if (SEIDkey>=0 && SEIDkey<384*512) {
00472                TimeCalibrationTable << 1.0 << '\t' << timeoffset[SEIDkey] << '\t';
00473              }
00474              else {
00475                TimeCalibrationTable << 1.0 << '\t' << 0.0 << '\t';
00476              }
00477              TimeCalibrationTable << 0 << '\t' << 0 << '\t';
00478              TimeCalibrationTable << 0 << '\t' << 0 << endl;
00479 
00480            } // loop over legal ends
00481            
00482          } // loop over strips in plane
00483          
00484          
00485    } // loop over planes
00486 }
00487 
00488 
00489 
00490 
00491