001
002
003
004
005
006
007
008 #include "MuonDigitTest/CSCDigitsTestAlg.h"
009 #include "EventInfo/EventInfo.h"
010 #include "EventInfo/EventID.h"
011
012
013
014 #include "GaudiKernel/MsgStream.h"
015 #include "GaudiKernel/AlgFactory.h"
016 #include "GaudiKernel/ITHistSvc.h"
017
018
019 #include "StoreGate/StoreGateSvc.h"
020
021 #include "Identifier/Identifier.h"
022
023 #include "GeoAdaptors/GeoMuonHits.h"
024
025 #include "MuonGeoModel/MuonDetectorManager.h"
026 #include "MuonGeoModel/CscReadoutElement.h"
027
028 #include "MuonIdHelpers/CscIdHelper.h"
029 #include "MuonSimEvent/CscHitIdHelper.h"
030
031 #include "MuonDigitContainer/CscDigitContainer.h"
032 #include "MuonDigitContainer/CscDigitCollection.h"
033 #include "MuonDigitContainer/CscDigit.h"
034
035 #include "MuonDetCluster/CscClusterCollection.h"
036 #include "MuonDetCluster/CscCluster.h"
037
038 #include "GeneratorObjects/McEventCollection.h"
039
040 #include "CLHEP/Vector/LorentzVector.h"
041 #include "CLHEP/Vector/ThreeVector.h"
042
043 #include <TH2D.h>
044 #include "TTree.h"
045 #include <string>
046
047
048 #include "CLHEP/Vector/ThreeVector.h"
049
050 #include "HepMC/GenEvent.h"
051
052 #include "TrackRecord/TrackRecord.h"
053 #include "TrackRecord/TrackRecordCollection.h"
054
055 #include "CLHEP/Geometry/Point3D.h"
056 #include "CLHEP/Geometry/Vector3D.h"
057
058
059 #include "MuonSimData/MuonSimData.h"
060 #include "MuonSimData/MuonSimDataCollection.h"
061
062 #include <algorithm>
063 #include <cmath>
064 #include <fstream>
065
066 const int max_ValDigitNumber =100 ;
067
068 const int maxDigits = 200;
069
070 const int maxFired = 8192;
071
072 using namespace MuonGM;
073
074
075
076
077
078
079 CSCDigitsTestAlg::CSCDigitsTestAlg(const std::string& name, ISvcLocator* pSvcLocator):
080 Algorithm(name, pSvcLocator),
081 m_pMuonMgr (0),
082 m_pCscIdHelper (0),
083 TestAlgInterface(),
084 log( msgSvc(), name )
085 {
086 declareProperty("WriteNtuple", m_writeNtuple = true);
087 declareProperty("McEventKey", m_key="TruthEvent");
088 declareProperty("NtupleFileName", m_ntupleFileName = "MuonDigits");
089 declareProperty("NtupleDirectoryName", m_ntupleDirName = "CSCDigitsValidation");
090 declareProperty("NtupleTreeName", m_ntupleTreeName = "CSCTestDigits");
091 declareProperty("DoCSCDigits", m_DoCSCDigits = true);
092 declareProperty("DoCSCClusters", m_DoCSCClusters = true);
093 declareProperty("Threshold", m_threshold = 500);
094 declareProperty("DumpTrackRecord", m_dumpTrackRecord = true);
095 m_event = 0;
096
097 log << MSG::INFO << "++++++++++++ CSCDigitsTestAlg created ++++++++++++"
098 << endreq;
099
100 }
101
102 CSCDigitsTestAlg::~CSCDigitsTestAlg()
103 {
104 log << MSG::INFO << " deleting CSCDigitsTestAlg " << endreq;
105 }
106
107
108
109
110
111 StatusCode CSCDigitsTestAlg::initialize()
112 {
113 log << MSG::INFO << " initializing CSCDigitsTestAlg " << endreq;
114
115
116
117
118 StatusCode sc = service("StoreGateSvc", m_sgSvc);
119
120 if (!sc.isSuccess() || 0 == m_sgSvc) {
121 log << MSG::ERROR << "CSCDigitsTestAlg: Could not find StoreGateSvc" << endreq;
122 return StatusCode::FAILURE;
123 }
124 else {
125 log << MSG::DEBUG << "Retrieved StoreGateSvc" << endreq;
126 }
127
128
129
130
131 sc = service("ActiveStoreSvc",m_activeStore);
132 if (sc.isFailure()) {
133 log<<MSG::ERROR << "Could not retrieve StoreGate ActiveStoreSvc!" <<endreq;
134 } else {
135 log<<MSG::DEBUG<< "Retrieved StoreGate ActiveStoreSvc." << endreq;
136 }
137
138
139
140 StoreGateSvc* detStore = 0;
141 sc = service( "DetectorStore", detStore );
142 if ( sc.isSuccess() )
143 {
144 sc = detStore->retrieve( m_pMuonMgr );
145 if ( sc.isFailure() )
146 {
147 log << MSG::ERROR
148 << " CSCDigitsTestAlg::initialize() : Could not find MuonDetDescrMgr "
149 << endreq;
150 }
151 else
152 {
153 m_pCscIdHelper = m_pMuonMgr->cscIdHelper();
154 }
155 }
156 else
157 {
158 log << MSG::ERROR
159 << " CSCDigitsTestAlg::initialize() : MuonDetDescrMgr not found in DetectorStore "
160 << endreq;
161 }
162
163
164 if (sc.isFailure()){
165 log <<MSG::ERROR << "Could not get DetectorStore"<<endreq;
166 return( StatusCode::FAILURE );
167 }
168
169
170
171
172 if (m_writeNtuple) {
173
174
175 std::string cscStreamAndPath = "/"+m_ntupleFileName+"/"+m_ntupleDirName;
176 std::string csctreeName = m_ntupleTreeName;
177
178 std::string csctreePath = cscStreamAndPath+"/"+csctreeName;
179
180 csctree = new TTree(TString(csctreeName), "Muon CSC Digits output");
181
182 StatusCode status;
183 status=ToolRootHistSvc()->regTree(csctreePath, csctree);
184
185 csctree->Branch("m_Validation_CSC_NumberOfDigits",&m_Val_CSC_NumberOfDigits,"m_Validation_CSC_NumberOfDigits/I");
186 csctree->Branch("m_Validation_CSC_RunNumber",&m_Val_CSC_RunNumber,"m_Validation_CSC_RunNumber/I");
187 csctree->Branch("m_Validation_CSC_EventNumber",&m_Val_CSC_EventNumber,"m_Validation_CSC_EventNumber/I");
188 csctree->Branch("m_Validation_CSC_StationName",&m_Validation_CSC_StationName,"m_Validation_CSC_StationName/I");
189 csctree->Branch("m_Validation_CSC_StationEta",&m_Validation_CSC_StationEta,"m_Validation_CSC_StationEta/I");
190 csctree->Branch("m_Validation_CSC_StationPhi",&m_Validation_CSC_StationPhi,"m_Validation_CSC_Statio/I");
191 csctree->Branch("m_Validation_CSC_ChamberLayer",&m_Validation_CSC_ChamberLayer,"m_Validation_CSC_ChamberLayer/I");
192 csctree->Branch("m_Validation_CSC_WireLayer",&m_Validation_CSC_WireLayer,"m_Validation_CSC_WireLayer/I");
193 csctree->Branch("m_Validation_CSC_MeasuresPhi",&m_Validation_CSC_MeasuresPhi,"m_Validation_CSC_MeasuresPhi/I");
194 csctree->Branch("m_Validation_CSC_Strip",&m_Validation_CSC_Strip,"m_Validation_CSC_Strip/I");
195 csctree->Branch("m_Validation_CSC_PosX",&m_Validation_CSC_PosX,"m_Validation_CSC_PosX/D");
196 csctree->Branch("m_Validation_CSC_PosY",&m_Validation_CSC_PosY,"m_Validation_CSC_PosY/D");
197 csctree->Branch("m_Validation_CSC_PosZ",&m_Validation_CSC_PosZ,"m_Validation_CSC_PosZ/D");
198 csctree->Branch("m_Validation_CSC_Charge",&m_Validation_CSC_Charge,"m_Validation_CSC_Charge/D");
199
200
201
202
203
204
205 if(status.isFailure()) {
206 log << MSG::DEBUG << "CSCDigitsTestAlg:: Unable to register TTreeTuple : " << csctreePath << endreq;
207 return status;
208 }
209
210 }
211
212 log << MSG::INFO << "CSCDigitsTestAlg:: Initialisation ended " << endreq;
213 return StatusCode::SUCCESS;
214
215 }
216
217
218
219
220
221 StatusCode CSCDigitsTestAlg::execute()
222 {
223
224 StatusCode sc;
225 log << MSG::INFO << " CSCDigitsTestAlg:: execute " << endreq;
226
227
228
229
230 const EventInfo* pevt;
231
232
233 if (storeGateSvc()->retrieve(pevt).isFailure()) {
234 log << MSG::FATAL << "CSCDigitsTestAlg::Could not find event" << endreq;
235 return(StatusCode::FAILURE);
236 }
237 else {
238 log << MSG::INFO << "CSCDigitsTestAlg::Found EventInfo in SG" << endreq;
239 }
240
241 if(isVerbose())
242 {
243 log << MSG::VERBOSE <<"CSCDigitsTestAlg::Processing EventInfo event #"<<pevt->event_ID()->event_number() << " run: " << pevt->event_ID()->run_number() << endreq;
244 }
245
246 clear_variables();
247
248 m_Val_CSC_EventNumber = pevt->event_ID()->event_number();
249 m_Val_CSC_RunNumber = pevt->event_ID()->run_number();
250
251 if(m_DoCSCDigits) sc = CSCvalidateDigits();
252 if(m_DoCSCClusters) sc = CSCvalidateClusters();
253
254
255 return StatusCode::SUCCESS;
256
257 }
258
259 StatusCode CSCDigitsTestAlg::CSCvalidateDigits()
260 {
261 StatusCode sc;
262 static MsgStream log(msgSvc(), name());
263 log << MSG::DEBUG << "in execute(): validateDigits" << endreq;
264
265 std::string key = "CSC_DIGITS";
266 const DataHandle <CscDigitContainer> container;
267 sc = m_sgSvc->retrieve(container,key);
268 if (sc.isFailure())
269 {
270 log << MSG::ERROR << " Cannot retrieve CSC Container " << endreq;
271 return StatusCode::FAILURE;
272 }
273
274
275
276
277 m_Val_CSC_NumberOfDigits=0;
278 for (CscDigitContainer::const_iterator it1_coll=container->begin(); it1_coll!=container->end(); ++it1_coll)
279 {
280 const CscDigitCollection* cscCollection = *it1_coll;
281
282 Identifier moduleId = cscCollection->identify();
283 if (m_pCscIdHelper->validElement(moduleId))
284 {
285 for(CscDigitCollection::const_iterator it1_digit=cscCollection->begin(); it1_digit!=cscCollection->end();++it1_digit)
286 {
287
288 const CscDigit* cscDigit = *it1_digit;
289 if (cscDigit->is_valid(m_pCscIdHelper) )
290 {
291 Identifier channelId = cscDigit->identify();
292
293 const MuonGM::CscReadoutElement * descriptor = m_pMuonMgr->getCscReadoutElement(channelId);
294 HepPoint3D pos = descriptor->stripPos(channelId);
295
296 cscDigit->charge();
297 int stationName = m_pCscIdHelper->stationName(channelId);
298 int stationEta = m_pCscIdHelper->stationEta(channelId);
299 int stationPhi = m_pCscIdHelper->stationPhi(channelId);
300 int clayer = m_pCscIdHelper->chamberLayer(channelId);
301 int wlayer = m_pCscIdHelper->wireLayer(channelId);
302 int measuresPhi = m_pCscIdHelper->measuresPhi(channelId);
303 int strip = m_pCscIdHelper->strip(channelId);
304
305 m_Validation_CSC_StationName=stationName;
306 m_Validation_CSC_StationEta=stationEta;
307 m_Validation_CSC_StationPhi=stationPhi;
308 m_Validation_CSC_ChamberLayer=clayer;
309 m_Validation_CSC_WireLayer=wlayer;
310 m_Validation_CSC_MeasuresPhi=measuresPhi;
311 m_Validation_CSC_Strip=strip;
312 m_Validation_CSC_PosX=pos.x();
313 m_Validation_CSC_PosY=pos.y();
314 m_Validation_CSC_PosZ=pos.z();
315 m_Validation_CSC_Charge=cscDigit->charge();
316 m_Val_CSC_NumberOfDigits++;
317
318 csctree->Fill();
319 }
320 }
321 }
322 }
323 log << MSG::DEBUG <<" Number of CSC Collections Accessed " <<endreq;
324
325 return StatusCode::SUCCESS;
326
327 }
328
329 StatusCode CSCDigitsTestAlg::CSCvalidateClusters()
330 {
331 StatusCode sc;
332 static MsgStream log(msgSvc(), name());
333 log << MSG::DEBUG << "in execute(): validateClusters" << endreq;
334
335 const DataHandle <CscClusterCollection> collection;
336
337 std::string colKey = "cscClusters";
338 sc = m_sgSvc->retrieve(collection,colKey);
339 if (sc.isFailure())
340 {
341 log << MSG::ERROR << "CSCDigitsTestAlg::Cannot retrieve CSC Cluster Collection " << endreq;
342 return StatusCode::FAILURE;
343 }
344
345 log << MSG::INFO << "Size of the collection is " << collection->size() << endreq;
346
347 bool is_muon = false;
348 HepPoint3D vertex = truth_info(is_muon);
349
350 long icount = 0;
351 float clusterSize;
352
353
354
355 if (is_muon)
356 {
357 for(CscClusterCollection::const_iterator beginCluster=collection->begin(); beginCluster!=collection->end();++beginCluster)
358 {
359 const CscCluster* cluster = *beginCluster;
360 HepPoint3D position = cluster->position();
361 Identifier clusterId = cluster->identify();
362 int measuresPhi = m_pCscIdHelper->measuresPhi(clusterId);
363 icount = m_pCscIdHelper->wireLayer(clusterId);
364 float dist = (float) getDistance(position,vertex,measuresPhi, clusterId);
365 clusterSize = 0.0;
366
367 for (std::vector<CscDigit*>::const_iterator clusItr = cluster->digit_iterator_begin();
368 clusItr!=cluster->digit_iterator_end();
369 ++clusItr)
370 {
371 if ( (*clusItr)->charge() > m_threshold) clusterSize++;
372 }
373
374 m_Validation_CSC_ClusterPosX=position.x();
375 m_Validation_CSC_ClusterPosY=position.y();
376 m_Validation_CSC_ClusterPosZ=position.z();
377 m_Validation_CSC_ClusterSize=clusterSize;
378 m_Validation_CSC_ClusterSigma=cluster->sigma();
379 }
380 log << MSG::INFO << "Number of Clusters in the collection is " << clusterSize << endreq;
381 log << MSG::INFO << "***********************************************************" << endreq;
382 }
383
384 return StatusCode::SUCCESS;
385 }
386
387
388
389
390 HepPoint3D CSCDigitsTestAlg::truth_info(bool& is_muon) const
391 {
392 StatusCode status = StatusCode::SUCCESS;
393 MsgStream log(msgSvc(), name());
394
395 HepPoint3D pvert = HepPoint3D(-10000.0,-10000.0,-10000.0);
396
397 const DataHandle <TrackRecordCollection> collTR;
398
399 std::string location = "MuonEntryLayer";
400 status = m_sgSvc->retrieve(collTR,location);
401
402 if (status.isFailure() )
403 {
404 log << MSG::ERROR << "Could not find TrackRecordCollection at " << location << " trying MuonEntryRecord" << endreq;
405
406 location = "MuonEntryRecord";
407 status = m_sgSvc->retrieve(collTR,location);
408 if (status.isFailure() )
409 {
410 log << MSG::ERROR << "Could not find TrackRecordCollection at " << location <<endreq;
411 is_muon = false;
412 }
413 }
414 else if (status.isSuccess())
415 {
416 log<<MSG::DEBUG<<"TrackRecordCollection retrieved: size = "<<collTR->size()<<" TrackRecord hits!" <<endreq;
417
418 TrackRecordCollection::const_iterator iterTR;
419
420 for (iterTR=collTR->begin();iterTR!=collTR->end();++iterTR)
421 {
422 if (m_dumpTrackRecord) {
423 log << MSG::INFO << "Track found with pdg id= " << (*iterTR)->GetPDGCode()
424 << " with energy "<< (*iterTR)->GetEnergy() << "Muon EntryRecord: Pt = "
425 << ((*iterTR)->GetMomentum()).perp() << endreq;
426 }
427 int pdgCode = (*iterTR)->GetPDGCode();
428 if ( pdgCode == 13 || pdgCode == -13 || pdgCode == 0) {
429 double x = ((*iterTR)->GetPosition()).x();
430 double y = ((*iterTR)->GetPosition()).y();
431 double z = ((*iterTR)->GetPosition()).z();
432 double rho = sqrt(x*x+y*y);
433 double phi = ((*iterTR)->GetMomentum()).phi();
434 double theta = ((*iterTR)->GetMomentum()).theta();
435 double eta = ((*iterTR)->GetMomentum()).eta();
436 double pt = ((*iterTR)->GetMomentum()).perp();
437 if (m_dumpTrackRecord) {
438 log << MSG::DEBUG << "EntryRecord: The momentum is " << pt << " MeV" << endreq;
439 log << MSG::DEBUG << "EntryRecord: The eta is " << eta << endreq;
440 log << MSG::DEBUG << "EntryRecord: The vertex is (mm)" << x << " " << y << " "
441 << z << " PDG code = " << pdgCode << endreq;
442 }
443 pvert = HepPoint3D (rho,theta,phi);
444 is_muon = true;
445 }
446 }
447 }
448 return pvert;
449 }
450
451 double CSCDigitsTestAlg::getDistance(HepPoint3D maxPos, HepPoint3D vertex,
452 int measuresPhi, Identifier id) {
453
454 MsgStream log(messageService(), name());
455
456 double dist = -1000.0;
457 double thetaStrip = maxPos.theta();
458 double phiStrip = maxPos.phi();
459 double thetaMC = vertex.y();
460 double phiMC = vertex.z();
461
462 double alpha = 11.59 * 3.141592654 / 180;
463 double rhoS = sqrt( maxPos.x()*maxPos.x()+maxPos.y()*maxPos.y() );
464
465 const MuonGM::CscReadoutElement * descriptor = m_pMuonMgr->getCscReadoutElement(id);
466 HepVector3D theNorm = descriptor->normal();
467 if (measuresPhi == 0) theNorm = descriptor->normal(id);
468 HepPoint3D normDir = HepPoint3D(theNorm.x(),theNorm.y(),theNorm.z());
469 HepPoint3D thePoint = findThePoint(vertex,maxPos,normDir);
470
471 double rhoMC = thePoint.perp();
472
473
474 HepPoint3D localCenter = descriptor->localPos(maxPos);
475 HepPoint3D localPoint = descriptor->localPos(thePoint);
476
477 if (measuresPhi==0) {
478 dist = (localPoint.z()-localCenter.z());
479
480 } else {
481 dist = (localPoint.y()-localCenter.y());
482 }
483
484 if (m_writeNtuple) {
485 log << MSG::INFO <<" Comparing MC Truth to GeoModel Position "
486 << " rhoMC = " << std::setiosflags(std::ios::fixed) << std::setprecision(3) << std::setw(7) << rhoMC
487 << " thetaMC = " << std::setiosflags(std::ios::fixed) << std::setprecision(3) << std::setw(7) << thetaMC
488 << " phiMC = " << std::setiosflags(std::ios::fixed) << std::setprecision(3) << std::setw(7) << phiMC
489 << " rhoStrip = " << std::setiosflags(std::ios::fixed) << std::setprecision(3) << std::setw(7) << rhoS
490 << " thetaStrip = " << std::setiosflags(std::ios::fixed) << std::setprecision(3) << std::setw(7) << thetaStrip
491 << " phiStrip = " << std::setiosflags(std::ios::fixed) << std::setprecision(3) << std::setw(7) << phiStrip
492 << " distance from track (mm) = " << setiosflags(std::ios::fixed) << std::setprecision(3) << std::setw(7) << dist
493 << endreq;
494 }
495 return dist;
496 }
497
498 HepPoint3D CSCDigitsTestAlg::findThePoint(HepPoint3D vertex, HepPoint3D maxPos, HepPoint3D norm) {
499
500
501
502 double B = 1.0e+15;
503 double thetaMC = vertex.y();
504 double phiMC = vertex.z();
505
506 HepPoint3D b = HepPoint3D(B*sin(thetaMC)*cos(phiMC), B*sin(thetaMC)*sin(phiMC), B*cos(thetaMC));
507 HepPoint3D thePoint = ( maxPos.dot(norm) / b.dot(norm) ) * b;
508 return thePoint;
509 }
510
511 StatusCode CSCDigitsTestAlg::finalize()
512 {
513 log << MSG::ERROR << " CSCDigitsTestAlg :: !!!! in finalize()" << endreq;
514
515 return StatusCode::SUCCESS;
516 }
517
518 ITHistSvc* CSCDigitsTestAlg::ToolRootHistSvc()
519
520 {
521 MsgStream log(msgSvc(), name());
522 StatusCode sc = service("THistSvc",m_rootsvc, true);
523 if( sc.isFailure() )
524 {
525 log << MSG::WARNING << ">>> Unable to locate the CSCDigitsTestAlg Histogram service" << endreq;
526 }
527 return m_rootsvc;
528 }
529
530
531 void CSCDigitsTestAlg::clear_variables()
532 {
533 m_Val_CSC_RunNumber=-1;
534 m_Val_CSC_EventNumber=-1;
535 m_Val_CSC_NumberOfDigits=-1;
536 m_Validation_CSC_StationName=-1;
537 m_Validation_CSC_StationEta=-1;
538 m_Validation_CSC_StationPhi=-1;
539 m_Validation_CSC_ChamberLayer=-1;
540 m_Validation_CSC_WireLayer=-1;
541 m_Validation_CSC_MeasuresPhi=-1;
542 m_Validation_CSC_Strip=-1;
543 m_Validation_CSC_PosX=-1;
544 m_Validation_CSC_PosY=-1;
545 m_Validation_CSC_PosZ=-1;
546 m_Validation_CSC_Charge=-1;
547 m_Validation_CSC_ClusterPosX=-1;
548 m_Validation_CSC_ClusterPosY=-1;
549 m_Validation_CSC_ClusterPosZ=-1;
550 m_Validation_CSC_ClusterSize=-1;
551 m_Validation_CSC_ClusterSigma=-1;
552 }
553
554
| Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems |
|
This page was automatically generated by the
LXR engine.
|
|