#include <WriteAlignmentModule.h>
Inheritance diagram for WriteAlignmentModule:
Public Member Functions | |
| WriteAlignmentModule () | |
| virtual | ~WriteAlignmentModule () |
| void | BeginJob () |
| JobCResult | Ana (const MomNavigator *mom) |
| void | EndJob () |
|
|
Definition at line 49 of file WriteAlignmentModule.cxx. References DbiWriter< T >::Abort(), atoi(), DbiWriter< T >::Close(), PlexScintMdlId::GetModule(), DbiResultPtr< T >::GetNumRows(), CalHelpers::GetPlane(), PlexPlaneId::GetPlane(), DbiResultPtr< T >::GetRow(), UgliGeomHandle::GetScintPlnHandleVector(), VldTimeStamp::GetSec(), DbiResultPtr< T >::GetValidityRec(), MSG, UgliGeomHandle::Print(), DbiWriter< T >::SetDbNo(), PlexStripEndId::SetEnd(), and strlen(). 00050 {
00051 int iteration=0;
00052 string path;
00053 const char* strENV_ALIGN_ITER = getenv("ENV_ALIGN_ITER");
00054 if ( strENV_ALIGN_ITER == 0 || strlen(strENV_ALIGN_ITER) == 0 )
00055 {
00056 MSG("Align", Msg::kFatal) << "ENV_ALIGN_ITER is not set!" <<endl
00057 << " Program is terminating now." <<endl;
00058 abort();
00059 } else
00060 iteration = atoi(strENV_ALIGN_ITER);
00061
00062 const char* strENV_ALIGN_DIR = getenv("ENV_ALIGN_DIR");
00063 if ( strENV_ALIGN_DIR == 0 || strlen(strENV_ALIGN_DIR) == 0 )
00064 {
00065 MSG("Align", Msg::kFatal) << "ENV_ALIGN_DIR is not set!"
00066 << " Program is terminating now." <<endl;
00067 abort();
00068 } else
00069 path = string(strENV_ALIGN_DIR);
00070
00071 MSG("Align", Msg::kInfo) << endl << "Iteration = " << iteration << endl
00072 << "Alignment path = " << path << endl;
00073
00074
00075 map<PlexStripEndId, Double_t> StripTPos;
00076 map<PlexStripEndId, Double_t> StripZPos;
00077 map<PlexStripEndId, Double_t> StripZRot;
00078 map<PlexScintMdlId, Double_t> MdlTPos;
00079 map<PlexScintMdlId, Double_t> MdlZPos;
00080 map<PlexScintMdlId, Double_t> MdlZRot;
00081
00082 VldTimeStamp vts(2004, 10, 1, 0, 0, 0);
00083 VldContext vc(DetectorType::kNear, SimFlag::kData, vts);
00084 UgliGeomHandle ugh(vc);
00085 ugh.Print();
00086
00087 vector<UgliScintPlnHandle> uph = ugh.GetScintPlnHandleVector();
00088 MSG("Align", Msg::kInfo) << "uph.size() = " << uph.size() << endl;
00089
00090 for(vector<UgliScintPlnHandle>::iterator it = uph.begin(); it != uph.end(); ++it){
00091
00092 if(!it->IsValid()){
00093 MSG("Align", Msg::kInfo) << "Error! UgliScintPlnHandle is not valid." << endl;
00094 return;
00095 }
00096
00097 Double_t zpos = it -> GetZ0();
00098 vector<UgliScintMdlHandle> umhv = it -> GetScintMdlHandleVector();
00099 vector<UgliStripHandle> ushv = it -> GetStripHandleVector();
00100
00101 //Get strip tpos and zpos
00102 for(vector<UgliStripHandle>::iterator iter = ushv.begin();
00103 iter != ushv.end(); ++iter){
00104
00105 if(!iter -> IsValid()){
00106 MSG("Align", Msg::kInfo) << "Error! UgliStripHandle is not valid." << endl;
00107 return;
00108 }
00109
00110 //explicitly setting stripend id to whole
00111 //might cause problems for the Far Detector
00112 PlexStripEndId plexid = iter -> GetSEId();
00113 plexid.SetEnd(StripEnd::kWhole);
00114 StripTPos[plexid] = iter -> GetTPos();
00115 StripZPos[plexid] = zpos;
00116 StripZRot[plexid] = iter -> GetZRotRelMdlRad();
00117 }
00118
00119 //Get module zpos and tpos
00120 for(vector<UgliScintMdlHandle>::iterator iter = umhv.begin();
00121 iter != umhv.end(); ++iter){
00122
00123 if(!iter->IsValid()){
00124 cout << "Error! UgliScintMdlHandle is not valid." << endl;
00125 return;
00126 }
00127
00128 PlexScintMdlId plexid = iter -> GetPlexScintMdlId();
00129 MdlTPos[plexid] = iter -> GetTPosRelPln();
00130 MdlZPos[plexid] = zpos;
00131 MdlZRot[plexid] = iter -> GetZRotRelPlnRad();
00132 }
00133 }
00134
00135 cout << "StripTPos.size() = " << StripTPos.size() << endl;
00136 cout << "StripZPos.size() = " << StripZPos.size() << endl;
00137 cout << "StripZRot.size() = " << StripZRot.size() << endl;
00138
00139 cout << "MdlTPos.size() = " << MdlTPos.size() << endl;
00140 cout << "MdlZPos.size() = " << MdlZPos.size() << endl;
00141 cout << "MdlZRot.size() = " << MdlZRot.size() << endl;
00142
00143 stringstream ss;
00144 ss << path << "/alignment_" << iteration << ".root";
00145
00146 TFile *hfile = new TFile(ss.str().c_str());
00147
00148 if(!hfile){
00149 MSG("Align", Msg::kError) << ss.str().c_str() << " does not exit! Bailing out." <<endl;
00150 return;
00151 }
00152
00153 //_____________________________________________________________________________________
00154 TTree *MdlResTree = dynamic_cast<TTree*>(hfile->Get("MdlRes"));
00155
00156 if(!MdlResTree){
00157 cout << "Could not find MdlRes tree in the root file. Aborting now."<<endl;
00158 abort();
00159 }
00160
00161 Double_t hmean=0.0, hentries=0.0;
00162 Int_t plane=0, module=0;
00163 UInt_t plexmdlid=0;
00164
00165 TBranch *hmeanBranch = MdlResTree -> GetBranch("hmean");
00166 TBranch *hentriesBranch = MdlResTree -> GetBranch("hentries");
00167 TBranch *planeBranch = MdlResTree -> GetBranch("plane");
00168 TBranch *moduleBranch = MdlResTree -> GetBranch("module");
00169 TBranch *plexmdlidBranch = MdlResTree -> GetBranch("plexmdlid");
00170
00171 hmeanBranch -> SetAddress(&hmean);
00172 hentriesBranch -> SetAddress(&hentries);
00173 planeBranch -> SetAddress(&plane);
00174 moduleBranch -> SetAddress(&module);
00175 plexmdlidBranch -> SetAddress(&plexmdlid);
00176
00177 double nb = 0.0;
00178 map<PlexScintMdlId, Double_t> ModuleResidual;
00179 int nmdl = (int)MdlResTree->GetEntries();
00180 for (int i = 0; i < nmdl; ++i){
00181 nb += MdlResTree->GetEntry(i);
00182 PlexScintMdlId plexid(plexmdlid);
00183
00184 if(i % 10 == 0)
00185 cout << "Got new module entry: " << plexid << endl
00186 << "hmean = " << hmean << ", hentries = " << hentries
00187 <<", plane = " << plane << ", module = " << module <<endl;
00188
00189 if(hentries < 10000)
00190 continue;
00191
00192 if(plexid.GetPlane() != plane || plexid.GetModule() != module){
00193 cout << "PlexScintMdlId " <<plexid <<" is not constructed properly!" <<endl;
00194 continue;
00195 }
00196
00197 pair<PlexScintMdlId, Double_t> p(plexid, hmean);
00198 ModuleResidual.insert(p);
00199 }
00200
00201 cout << "Read " << nb <<" bytes for " <<ModuleResidual.size() << " module residuals"<< endl;
00202
00203 //_____________________________________________________________________________________
00204 TTree *MdlRotResTree = dynamic_cast<TTree*>(hfile->Get("MdlRotRes"));
00205
00206 if(!MdlResTree){
00207 cout << "Could not find MdlRes tree in the root file. Aborting now."<<endl;
00208 abort();
00209 }
00210
00211 Double_t par1=0.0, par1err=0.0;
00212
00213 TBranch *par1BranchRot = MdlRotResTree -> GetBranch("par1");
00214 TBranch *par1errBranchRot = MdlRotResTree -> GetBranch("par1err");
00215 TBranch *hmeanBranchRot = MdlRotResTree -> GetBranch("hmean");
00216 TBranch *hentriesBranchRot = MdlRotResTree -> GetBranch("hentries");
00217 TBranch *planeBranchRot = MdlRotResTree -> GetBranch("plane");
00218 TBranch *moduleBranchRot = MdlRotResTree -> GetBranch("module");
00219 TBranch *plexmdlidBranchRot = MdlRotResTree -> GetBranch("plexmdlid");
00220
00221
00222 par1BranchRot -> SetAddress(&par1);
00223 par1errBranchRot -> SetAddress(&par1err);
00224 hmeanBranchRot -> SetAddress(&hmean);
00225 hentriesBranchRot -> SetAddress(&hentries);
00226 planeBranchRot -> SetAddress(&plane);
00227 moduleBranchRot -> SetAddress(&module);
00228 plexmdlidBranchRot -> SetAddress(&plexmdlid);
00229
00230 map<PlexScintMdlId, Float_t> ModuleZRotOffset;
00231 map<int, vector<Float_t> > PlaneZRotOffsets;
00232
00233 nmdl = (int)MdlRotResTree->GetEntries();
00234 for (int i = 0; i < nmdl; ++i){
00235 nb += MdlRotResTree->GetEntry(i);
00236 PlexScintMdlId plexid(plexmdlid);
00237
00238 if(i % 10 == 0)
00239 cout << "Got new rotational module entry: " << plexid << endl
00240 << "par1 = " << par1 << ", hentries = " << hentries
00241 <<", plane = " << plane << ", module = " << module <<endl;
00242
00243 if(hentries < 10000)
00244 continue;
00245
00246 if(plexid.GetPlane() != plane || plexid.GetModule() != module){
00247 cout << "PlexScintMdlId " <<plexid <<" is not constructed properly!" <<endl;
00248 continue;
00249 }
00250
00251 if(fabs(par1) < 0.0001)
00252 continue;
00253
00254 if(fabs(par1err/par1) > 0.05)
00255 continue;
00256
00257 Float_t angle = (Float_t) par1;
00258 angle = atan(angle);
00259 pair<PlexScintMdlId, Float_t> p(plexid, angle);
00260 ModuleZRotOffset.insert(p);
00261
00262 PlaneZRotOffsets[plane].push_back(angle);
00263 }
00264
00265 cout << "Read " << nb <<" bytes for " << ModuleResidual.size() << " module residuals"<< endl;
00266
00267 hfile -> Close();
00268
00269 map<int, Float_t> PlaneZRotOffset;
00270 for(map<int, vector<Float_t> >::iterator pit = PlaneZRotOffsets.begin();
00271 pit != PlaneZRotOffsets.end(); ++pit)
00272 {
00273 vector<Float_t> &zv = pit -> second;
00274
00275 if(zv.size() < 3)
00276 continue;
00277
00278 Float_t tmp = 0.0;
00279 Float_t size = Float_t(zv.size());
00280 for(unsigned int i = 0; i < zv.size(); ++i)
00281 tmp += zv[i];
00282 tmp = tmp /size;
00283
00284 if(fabs(tmp) > 0.0001)
00285 PlaneZRotOffset[pit->first] = tmp;
00286 }
00287
00288
00289 vector<UgliDbiScintMdl> DbiScintMdlRow;
00290 vector<UgliDbiScintPln> DbiScintPlnRow;
00291
00292 VldTimeStamp timenow;
00293 VldTimeStamp begtime(timenow.GetSec() - 365*24*60*60, 0); //one year back from now
00294 VldTimeStamp endtime(timenow.GetSec() + 365*24*60*60, 0); //one year in future from now
00295
00296 VldRange vldrange(DetectorType::kNear, SimFlag::kData, begtime, endtime, "DetectorAlignment");
00297
00298 MSG("Align", Msg::kInfo) << "Database writer time range: " << vldrange << endl;
00299
00300 VldContext vldcnow(DetectorType::kNear, SimFlag::kData, timenow);
00301
00302 MSG("Align", Msg::kInfo) << "Reading UgliDbiScintMdl from " << vldcnow << endl;
00303 DbiResultPtr<UgliDbiScintMdl> uglimdlptr(vldcnow);
00304 DbiResultPtr<UgliDbiScintPln> ugliplnptr(vldcnow);
00305
00306 const DbiValidityRec* vldrec = uglimdlptr.GetValidityRec();
00307 const int nmdlrows = uglimdlptr.GetNumRows();
00308 const int nplnrows = ugliplnptr.GetNumRows();
00309
00310 MSG("Align", Msg::kInfo) << "UgliDbiScintMdl pointer has "<< uglimdlptr.GetNumRows() << "rows" << endl;
00311 MSG("Align", Msg::kInfo) << "UgliDbiScintPln pointer has "<< ugliplnptr.GetNumRows() << "rows" << endl;
00312 MSG("Align", Msg::kInfo) << "UgliDbiScintMdl pointer validity range: "<< vldrec -> GetVldRange() << endl;
00313
00314 if(uglimdlptr.GetNumRows() < 1){
00315 MSG("Align", Msg::kError) << "UgliDbiScintMdl has no rows. Bailing out now." <<endl;
00316 return;
00317 }
00318
00319 MSG("Align", Msg::kInfo) << "Write to database yes/no? ";
00320 string response;
00321 cin >> response;
00322
00323 int irow = 0;
00324 int nmovedtposmdl = 0;
00325 int nmovedzrotmdl = 0;
00326 int nmovedzrotpln = 0;
00327
00328 do{
00329
00330 UInt_t nextaggNo = 0;
00331 UInt_t currentaggNo = uglimdlptr.GetRow(irow) -> GetAggregateNo();
00332 Dbi::Task task = 0;
00333 DbiWriter<UgliDbiScintMdl> dbiwriter(vldrange, currentaggNo, task, timenow);
00334 dbiwriter.SetDbNo(0);
00335
00336 do{
00337
00338 if(uglimdlptr.GetRow(irow)){
00339 uglimdlptr.GetRow(irow) -> Print();
00340 PlexScintMdlId plexid = uglimdlptr.GetRow(irow) -> GetScintMdlId();
00341 Float_t width = uglimdlptr.GetRow(irow) -> GetWidth();
00342 Float_t tpos = uglimdlptr.GetRow(irow) -> GetTPosRelPln();
00343 Float_t lpos = uglimdlptr.GetRow(irow) -> GetLPosRelPln();
00344 Float_t zrot = uglimdlptr.GetRow(irow) -> GetZRotRelPlnRad();
00345 Float_t clear_east = uglimdlptr.GetRow(irow) -> GetClearLenEast();
00346 Float_t clear_west = uglimdlptr.GetRow(irow) -> GetClearLenWest();
00347 Float_t wls_east = uglimdlptr.GetRow(irow) -> GetWlsLenEast();
00348 Float_t wls_west = uglimdlptr.GetRow(irow) -> GetWlsLenWest();
00349
00350 if(fabs(tpos - MdlTPos[plexid]) > 0.0000001){
00351 MsgFormat ffmt("%13.12f");
00352 MSG("Align", Msg::kError) <<"For " << plexid << " tpos varies between UgliGeometry cache file and db by "
00353 << ffmt(tpos - MdlTPos[plexid]) << endl;
00354 abort();
00355 }
00356
00357 if(fabs(zrot - MdlZRot[plexid]) > 0.0000001){
00358 MsgFormat ffmt("%13.12f");
00359 MSG("Align", Msg::kError) <<"For " << plexid << " zrot varies between UgliGeometry cache file and db by "
00360 << ffmt(zrot - MdlZRot[plexid]) << endl;
00361 abort();
00362 }
00363
00364 map<PlexScintMdlId, Double_t>::const_iterator pit = ModuleResidual.find(plexid);
00365 //if this module should be moved - move it!
00366 if(pit != ModuleResidual.end() && plexid.GetPlane() < 140){
00367 tpos = tpos + pit -> second;
00368 nmovedtposmdl++;
00369 MSG("Align", Msg::kInfo) << "Adjusting tpos of module " << plexid << " by " << pit->second <<" meters" << endl;
00370 }
00371
00372 map<PlexScintMdlId, Float_t>::const_iterator zit = ModuleZRotOffset.find(plexid);
00373 //if this module should be moved - move it!
00374 if(false && zit != ModuleZRotOffset.end() && plexid.GetPlane() < 140){
00375 nmovedzrotmdl++;
00376 MSG("Align", Msg::kInfo) << "Adjusting zrot of module " << plexid << " by " << zit->second <<" radians" << endl;
00377 }
00378
00379 zrot = 0.0;
00380 UgliDbiScintMdl dbrow(plexid, width, tpos, lpos, zrot,
00381 clear_east, clear_west, wls_east, wls_west);
00382
00383 dbiwriter << dbrow;
00384
00385 }
00386 irow++;
00387 if(uglimdlptr.GetRow(irow))
00388 nextaggNo = uglimdlptr.GetRow(irow) -> GetAggregateNo();
00389 else
00390 nextaggNo = nmdlrows;
00391
00392 } while(currentaggNo == nextaggNo);
00393
00394 if(response == "yes_mdl") {
00395 if(dbiwriter.Close())
00396 MSG("Align", Msg::kInfo) << "DbiWriter successfully wrote aggNo " << currentaggNo << endl;
00397 else
00398 MSG("Align", Msg::kError) << "DbiWriter failed to write aggNo " << currentaggNo << endl;
00399 } else
00400 dbiwriter.Abort();
00401
00402 } while( uglimdlptr.GetRow(irow) );
00403
00404
00405 for(int i = 0; i < nplnrows; ++i)
00406 {
00407 if(!ugliplnptr.GetRow(i))
00408 continue;
00409
00410 UInt_t currentaggNo = ugliplnptr.GetRow(i) -> GetAggregateNo();
00411 Dbi::Task task = 0;
00412 DbiWriter<UgliDbiScintPln> dbiwriter(vldrange, currentaggNo, task, timenow);
00413 dbiwriter.SetDbNo(0);
00414
00415 Int_t plane = ugliplnptr.GetRow(i) -> GetPlane();
00416 PlexPlaneId plnid = ugliplnptr.GetRow(i) -> GetPlaneId();
00417 Float_t thick = ugliplnptr.GetRow(i) -> GetThickness();
00418 Float_t gap = ugliplnptr.GetRow(i) -> GetGap();
00419 Float_t x0rel = ugliplnptr.GetRow(i) -> GetX0RelSteel();
00420 Float_t y0rel = ugliplnptr.GetRow(i) -> GetY0RelSteel();
00421 Float_t zrot = ugliplnptr.GetRow(i) -> GetZRotRelSteelRad();
00422
00423 map<int, Float_t>::const_iterator zit = PlaneZRotOffset.find(plane);
00424 //if this module should be moved - move it!
00425 if(false && zit != PlaneZRotOffset.end() && plane < 140){
00426 zrot = zrot + zit->second;
00427 nmovedzrotpln++;
00428 MSG("Align", Msg::kInfo) << "Adjusting zrot of plane " << plane << " by " << zit->second <<" radians" << endl;
00429 }
00430
00431 if(zrot < 0.0)
00432 zrot = -0.785398;
00433 else
00434 zrot = 0.785398;
00435
00436 UgliDbiScintPln dbrow(plnid, thick, gap, x0rel, y0rel, zrot);
00437 dbiwriter << dbrow;
00438
00439 if(response == "yes_plane") {
00440 if(dbiwriter.Close())
00441 MSG("Align", Msg::kInfo) << "DbiWriter successfully wrote aggNo " << currentaggNo << endl;
00442 else
00443 MSG("Align", Msg::kError) << "DbiWriter failed to write aggNo " << currentaggNo << endl;
00444 } else
00445 dbiwriter.Abort();
00446 }
00447
00448 cout << "nmovedtposmdl = " << nmovedtposmdl << endl;
00449 cout << "nmovedzrotmdl = " << nmovedzrotmdl << endl;
00450 cout << "nmovedzrotpln = " << nmovedzrotpln << endl;
00451
00452 }
|
|
|
Definition at line 13 of file WriteAlignmentModule.h. 00013 {};
|
|
|
Implement this for read only access to the MomNavigator Reimplemented from JobCModule. Definition at line 455 of file WriteAlignmentModule.cxx. 00456 {
00457 //fake Ana function - just to supress complilation warning message
00458 if(mom)
00459 return JobCResult::kAOK;
00460 return JobCResult::kAOK;
00461 }
|
|
|
Implement for notification of begin of job Reimplemented from JobCModule. Definition at line 14 of file WriteAlignmentModule.h. 00014 {};
|
|
|
Implement for notification of end of job Reimplemented from JobCModule. Definition at line 16 of file WriteAlignmentModule.h. 00016 {};
|
1.3.9.1