// class MacroCellGeometry
// $Id: MacroCellGeometry.cc,v 1.14 2005/02/24 11:00:27 pesce Exp $
// description of the geometry of a MacroCell
//
#include "MacroCellGeometry.hh"
#include "PmtGeometry.hh"
#include "OpticalAdaptor.hh"
#include "MacroCell.hh"
#include "ElementaryCell.hh"
ClassImp(MacroCellGeometry)
//______________________________________________________________________________
MacroCellGeometry::MacroCellGeometry( MacroCell *parent) : pCell(parent) {
//
// ctor
//
}
//______________________________________________________________________________
MacroCellGeometry::~MacroCellGeometry() {
//
// dtor
//
}
//______________________________________________________________________________
Double_t MacroCellGeometry::IsHit(const Photon &p) const {
//
//
//
// check if p is // to the MC plane
if ( p.dir*GetNormal() > -TOLERANCE ) return -1;
EVector dummy=p.pos-GetCenter();
EVector intPoint=dummy-p.dir*((dummy*GetNormal())/(p.dir*GetNormal()));
if ( intPoint*fXAxis < fXMax && intPoint*fXAxis > fXMin &&
intPoint*fYAxis < fYMax && intPoint*fYAxis > fYMin ){
#ifdef DEBUG
cout << "\nMacroCellGeometry::IsHit\n"
<< "p.pos = " << p.pos << " p.dir = " << p.dir << endl;
cout << "fCenter = " << GetCenter() << " GetNormal() = " << GetNormal() << endl;
cout << "fXAxis = " << fXAxis << " fYAxis = " << fYAxis << endl;
cout << "fXMax = " << fXMax << " fXMin = " << fXMin
<< " fYMax = " << fYMax << " fYMin = " << fYMin << endl;
cout << "intPoint = " << intPoint+GetCenter()<< " intPoint*fXAxis = " << intPoint*fXAxis << " intPoint*fYAxis = " << intPoint*fYAxis << endl;
cout << "DL = " << (intPoint-dummy).Mag() << endl;
#endif /* DEBUG */
return (intPoint-dummy).Mag();
}
return -1;
}
//______________________________________________________________________________
Photomultiplier* MacroCellGeometry::FindHitPmt(const Photon &p) const{
//
// finds which pmt has been hit, if exists
//
const vector<Photomultiplier*> &PMTs=*(pCell->Pmts());
Photomultiplier *HitPMT=NULL;
for(size_t j=0; j<PMTs.size(); j++) {
// bring p on pmt cathode
Double_t DL = PMTs[j]->Geometry()->GetOA()->IsHit(p);
if( DL > 0) {
HitPMT=PMTs[j];
break;
}
}
return HitPMT;
}
//______________________________________________________________________________
void MacroCellGeometry::Compute() {
//
// Compute internal infos
//
if ( !Cell() ) {
FatalError("Invalid MacroCell pointer in MacroCellGeometry");
}
// number of elementary cells
Int_t nECs = Cell()->ECs()->size();
EVector n(0.,0.,0.), y_axis(0., 0., 0.);
// compute normal vector and the y axis
for(Int_t i=0; i<nECs; i++) {
const ElementaryCell *pEC = Cell()->GetEC(i);
if ( pEC ) {
n += pEC->GetNormal();
y_axis += pEC->GetYAxis();
} else {
FatalError("Invalid Elementary Cell in MacroCellGeometry");
}
}
fNormal = n.Unit();
fXAxis = (y_axis.Cross( n )).Unit();
fYAxis = fNormal.Cross( fXAxis );
// compute center
EVector center(0.,0.,0.);
for(Int_t i=0; i < nECs; i++) {
ElementaryCell *ec = Cell()->GetEC(i);
center += ec->GetCenter();
}
center *= (1./(double)Cell()->GetNumECs());
// compute corners and radius
for(Int_t i=0; i < nECs; i++) {
EVector d = Cell()->GetEC(i)->GetCenter() - center;
// d projection on the MC plane
if ( !i ){
fXMin = d.Dot(fXAxis);
fXMax = d.Dot(fXAxis);
fYMin = d.Dot(fYAxis);
fYMax = d.Dot(fYAxis);
} else {
if ( d.Dot(fXAxis) < fXMin ) fXMin = d.Dot(fXAxis);
if ( d.Dot(fXAxis) > fXMax ) fXMax = d.Dot(fXAxis);
if ( d.Dot(fYAxis) < fYMin ) fYMin = d.Dot(fYAxis);
if ( d.Dot(fYAxis) > fYMax ) fYMax = d.Dot(fYAxis);
}
}
// shift center of 50 mm along normal to keep the square in front of the macrocell
fCenter = center + 50. * fNormal.Unit()*mm;
// shift points out by two full PMT sides
fXMin-=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();
fXMax+=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();
fYMin-=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();
fYMax+=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();
}
void MacroCellGeometry::Draw() const{
//FIXME: to be done (call EC->Draw)
}