using namespace std;
# include <fstream>
# include <math.h>
# include "TFile.h"
# include "TH1.h"
# include "TMatrix.h"
# include "TNtuple.h"
# include "hades.h"
# include "hdebug.h"
# include "hdetector.h"
# include "hevent.h"
# include "hlocation.h"
# include "hmatrixcategory.h"
# include "hmatrixcatiter.h"
# include "hmdccal2.h"
# include "hmdccal3.h"
# include "hmdccaltable.h"
# include "hmdcdetector.h"
# include "hmdcmodulegeometry.h"
# include "hmdcgeomstruct.h"
# include "hmdchit.h"
# include "hmdchitaux.h"
# include "hmdchitf.h"
# include "hmdchitfpar.h"
# include "hmdchitsfilter.h"
# include "hmdctargetgeometry.h"
# include "hreconstructor.h"
# include "hruntimedb.h"
# include "hspectrometer.h"
# include "mdcsdef.h"
ClassImp(HMdcHitF)
void HMdcHitF :: fit(HMdcHitAux* data){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::fit");
#endif
Double_t a11 = 0.0, a12 = 0.0, a13 = 0.0, a14 = 0.0;
Double_t a22 = 0.0, a23 = 0.0, a24 = 0.0;
Double_t a33 = 0.0, a34 = 0.0, a44 = 0.0;
Double_t b11, b12, b13, b14, b22, b23, b24, b33, b34, b44;
Double_t tt1 = 0.0, tt2 = 0.0, tt3 = 0.0, tt4 = 0.0;
Double_t s1, s2, s3, s4, aa, bb, th;
Double_t thx, xx, sss1, sss2, d1, d2, r311, r312, r321, r322;
Double_t r511, r512, r522, sig2, error4[10], xg, weight;
Double_t thfit, chi;
Int_t nhit=0;
Double_t theta[6];
const Int_t nplane = 6;
const Int_t numHits = data->getNumHits();
Int_t j=0,i;
Int_t listPlanes[6];
for(i=0; i< nplane; i++){
if (data->getCal(i)!=NULL){
theta[j] = data->getCal(i)->getPos();
listPlanes[j] = i;
j++;
}
}
#if DEBUG_LEVEL>3
ofstream fitfile("fit.txt",ios::app);
fitfile << "suceso " << fEventId << endl;
for(i=0; i<numHits; i++){
fitfile << i << '\t' << listPlanes[i] << '\t' << theta[i] << endl;
}
fitfile << endl;
fitfile.close();
#endif
Int_t sector = fLoc[0];
Int_t module = fLoc[1];
Double_t zpoint = parContainer->getZHit(sector,module);
Double_t cosn[nplane], sinn[nplane], zplane[nplane], sgpln[nplane];
for(i=0; i<numHits; i++){
zplane[i] = geo->getZPlane(sector,module,listPlanes[i]);
cosn[i] = geo->getCosn(sector,module,listPlanes[i]);
sinn[i] = geo->getSinn(sector,module,listPlanes[i]);
sgpln[i] = geo->getSgpln(sector,module,listPlanes[i]);
}
#if DEBUG_LEVEL>3
ofstream output("fit.txt",ios::app);
output << "number of planes=" << nplane << '\n';
output << "number" << '\t' <<"zplane"<< '\t'<< "cos[i]" <<'\t' <<
"sin[i]"<< '\t'<< "hits [i]" <<'\t' << "rms[i]" << '\n';
for (i=0; i<nplane; i++){
output << i+1 <<'\t' << zplane[i] <<'\t' << cosn[i] <<'\t' << sinn[i] <<
'\t' <<theta[i] << '\t'<< sgpln[i] << '\n';
}
# endif
if (zpoint > -1000.0){
nhit = numHits;
data->setPointZ(zpoint);
}
else {
sss1 = 0.0;
sss2 = 0.0;
for(i=0; i < nplane; i++){
nhit = nhit + 1;
weight = 1.0/(sgpln[i]*sgpln[i]);
sss1 = sss1 + zplane[i] * weight;
sss2 = sss2 + weight;
data->setPointZ(sss1/sss2);
}
}
xg = data->getPointZ();
data->setChi(1.0E05);
for (i=0; i<numHits; i++){
aa = -cosn[i];
bb = sinn[i];
xx = zplane[i] - xg;
th = theta[i];
sig2 = 1.0/(sgpln[i] * sgpln[i]);
thx = th * xx;
s1 = aa * aa * sig2;
s2 = bb * bb * sig2;
s3 = aa * bb * sig2;
aa = aa * sig2;
bb = bb * sig2;
s4 = xx * xx;
a11 = a11 + s1;
a12 = a12 + s1 * xx;
a13 = a13 + s3;
a14 = a14 + s3 * xx;
a22 = a22 + s4 * s1;
a24 = a24 + s3 * s4;
a33 = a33 + s2;
a34 = a34 + s2 * xx;
a44 = a44 + s2 * s4;
tt1 = tt1 + aa * th;
tt2 = tt2 + aa * thx;
tt3 = tt3 + bb * th;
tt4 = tt4 + thx * bb;
}
a23 = a14;
d1 = 1.0 / (a12*a12 - a11*a22);
r311 = (a22*a13 - a12*a23) * d1;
r312 = (a22*a14 - a12*a24) * d1;
r321 = (a11*a23 - a12*a13) * d1;
r322 = (a11*a24 - a12*a14) * d1;
r511 = a13*r311 + a23*r321 + a33;
r512 = a13*r312 + a23*r322 + a34;
r522 = a14*r312 + a24*r322 + a44;
d2 = 1.0 / (r511*r522 - r512*r512);
b33 = r522 * d2;
b34 = -r512 * d2;
b44 = r511 * d2;
b13 = r311 * b33 + r312 * b34;
b14 = r311 * b34 + r312 * b44;
b23 = r321 * b33 + r322 * b34;
b24 = r321 * b34 + r322 * b44;
b11 = r311*b13 + r312*b14 - a22 * d1;
b12 = r311*b23 + r312*b24 + a12 * d1;
b22 = r321*b23 + r322*b24 - a11 * d1;
error4[0] = b44;
error4[1] = b24;
error4[2] = b34;
error4[3] = b14;
error4[4] = b22;
error4[5] = b23;
error4[6] = b12;
error4[7] = b33;
error4[8] = b13;
error4[9] = b11;
data->setErrorX(b33);
data->setErrorY(b11);
data->setErrorSlopeX(b44);
data->setErrorSlopeY(b22);
data->setPointY(b11*tt1 + b12*tt2 + b13*tt3 + b14*tt4);
data->setSlopeY(b12*tt1 + b22*tt2 + b23*tt3 + b24*tt4);
data->setPointX(b13*tt1 + b23*tt2 + b33*tt3 + b34*tt4);
data->setSlopeX(b14*tt1 + b24*tt2 + b34*tt3 + b44*tt4);
data->setCovElement(0,1,b13);
data->setCovElement(0,2,b34);
data->setCovElement(0,3,b23);
data->setCovElement(1,2,b14);
data->setCovElement(1,3,b12);
data->setCovElement(2,3,b24);
if(nhit == 4) data->setChi(0.0);
else{
chi = 0.0;
for (i=0; i<numHits; i++){
aa = -cosn[i];
bb = sinn[i];
xx = zplane[i] - xg;
thfit = aa*(data->getPointY() + data->getSlopeY()*xx) +
bb*(data->getPointX() + data->getSlopeX()*xx);
chi = chi + pow((thfit-theta[i]),2)/(pow(sgpln[i],2));
}
data->setChi(chi/(nhit-4));
}
# if DEBUG_LEVEL>3
ofstream output("fit.txt",ios::app);
output << "event " << fEventId << endl;
data->print(output);
output << '\n' << "Check results >>>" << '\n';
Double_t thetth, deviat, xochi2;
for (i=0; i<nplane; i++){
thetth = sinn[i] * (data->getPointX() + (zplane[i]-data->getPointZ())*
data->getSlopeX()) -
cosn[i] * (data->getPointY() + (zplane[i]-data->getPointZ())*
data->getSlopeY());
deviat = fabs((thetth-theta[i])/thetth)*100;
xochi2 = pow((thetth-theta[i]),2)/(pow(sgpln[i],2));
output << "Plane " << (i+1) << '\n';
output << "Fitted Hit " << thetth << '\n';
output << "Real Hit " << theta [i] << '\n';
output << "Deviation " << deviat << '\n';
output << "chi2 contribution " << xochi2 << '\n' << '\n';
}
output.close();
# endif
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::fit");
#endif
}
HMdcHitF :: HMdcHitF(void){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::HMdcHitF()");
#endif
mdcCal = NULL;
mdcSegments = NULL;
fCalCat = NULL;
fList = NULL;
target = kTRUE;
slopeCorrection = kFALSE;
fLoc.setNIndex(2);
fEventId = 0;
fHitCat = NULL;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::HMdcHitF()");
#endif
}
HMdcHitF :: HMdcHitF(const Text_t* name,const Text_t* title) : HReconstructor(name, title){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::HMdcHitF(name,title)");
#endif
mdcCal = NULL;
mdcSegments = NULL;
fCalCat = NULL;
fList = NULL;
fHitCat = NULL;
target = kTRUE;
slopeCorrection = kFALSE;
fLoc.setNIndex(2);
fEventId = 0;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::HMdcHitF(name,title)");
#endif
}
HMdcHitF :: HMdcHitF(HMatrixCategory* data, HMdcModuleGeometry* mdcGeo, HCategory* seg){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::HMdcHitF(HMatrixCategory*, HMdcModuleGeometry* , HCategory*)");
#endif
geo = mdcGeo;
mdcCal = data;
fCalCat = NULL;
mdcSegments = seg;
fLoc.setNIndex(2);
fList=NULL;
fHitCat = NULL;
target = kTRUE;
slopeCorrection = kFALSE;
fEventId=0;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::HMdcHitF(HMatrixCategory*, HMdcModuleGeometry* , HCategory*)");
#endif
}
HMdcHitF :: HMdcHitF(HCategory* data, HMdcModuleGeometry* mdcGeo, HCategory* seg){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::HMdcHitF(HCategory*, HMdcModuleGeometry* , HCategory*)");
#endif
geo = mdcGeo;
mdcCal = data;
mdcSegments = seg;
fCalCat = NULL;
fLoc.setNIndex(2);
fList=NULL;
fHitCat = NULL;
target = kTRUE;
slopeCorrection = kFALSE;
fEventId=0;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::HMdcHitF(HCategory*, HMdcModuleGeometry* , HCategory*)");
#endif
}
void HMdcHitF :: setParContainer(void){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::setParContainer");
#endif
geo=(HMdcModuleGeometry*)gHades->getRuntimeDb()->getContainer("MdcModuleGeometry");
parContainer=(HMdcHitFPar*)gHades->getRuntimeDb()->getContainer("MdcHitFPar");
if(target)
targetGeo=(HMdcTargetGeometry*)gHades->getRuntimeDb()->getContainer("MdcTargetGeometry");
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::setParContainer");
#endif
}
Bool_t HMdcHitF :: init(void){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::init");
#endif
setParContainer();
fCalCat = gHades->getCurrentEvent()->getCategory(catMdcCal2);
if (!fCalCat) {
fCalCat=gHades->getSetup()->getDetector("Mdc")->buildCategory(catMdcCal2);
if (!fCalCat) return kFALSE;
else gHades->getCurrentEvent()->addCategory(catMdcCal2,fCalCat,"Mdc");
}
fHitCat = gHades -> getCurrentEvent() -> getCategory (catMdcHit);
if (! fHitCat ){
fHitCat = gHades->getSetup()->getDetector("Mdc")->buildCategory(catMdcHit);
if (!fHitCat) return kFALSE;
else gHades->getCurrentEvent()->addCategory(catMdcHit,fHitCat,"Mdc");
}
TArrayI* ind = new TArrayI(4);
HMdcGeomStruct* p=
(HMdcGeomStruct*)gHades->getRuntimeDb()->getContainer("MdcGeomStruct");
if (!p) return 0;
p->getMaxIndices(ind);
Int_t nSizes=0;
Int_t* sizes;
mdcCal = gHades->getCurrentEvent()->getCategory(catMdcCalHit);
if(!mdcCal){
nSizes=ind->GetSize();
sizes=new Int_t[nSizes];
for (Int_t i=0;i<nSizes-1;i++) sizes[i]= ind->At(i)+1;
sizes[nSizes-1] = 2 * (ind->At(nSizes-1)+1);
mdcCal = new HMatrixCategory(getCalHitClassName(),nSizes,sizes,0.5);
gHades->getCurrentEvent()->addCategory(catMdcCalHit,mdcCal,"Mdc");
delete [] sizes;
}
mdcSegments = gHades->getCurrentEvent()->getCategory(catMdcSegment);
if(!mdcSegments){
sizes=new Int_t[nSizes-1];
for (Int_t i=0;i<nSizes-2;i++) sizes[i] = ind->At(i)+1;
sizes[nSizes-2] = 1200;
mdcSegments = new HMatrixCategory("HMdcHitAux",nSizes-1,sizes,0.5);
mdcSegments->setPersistency(kFALSE);
gHades->getCurrentEvent()->addCategory(catMdcSegment,mdcSegments,"Mdc");
delete [] sizes;
}
delete ind;
fActive = kTRUE;
Int_t layers = ((HMdcDetector* )(gHades->getSetup()->getDetector("Mdc")))->getMaxComponents();
fList = new HMdcCalLinkList*[layers];
for(Int_t i=0; i< layers; i++) fList[i] = new HMdcCalLinkList;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::init");
#endif
return kTRUE;
}
void HMdcHitF :: setLoc(HLocation &location){
fLoc.setNIndex(2);
fLoc.readIndexes(location);
if(fLoc[1] > 1) target = kFALSE;
}
HMdcHitF :: ~HMdcHitF(void){
#if DEBUG_LEVEL >2
gDebuger->enterFunc("HMdcHitF::~HMdcHitF");
#endif
Int_t layers = 6;
for(Int_t i = 0; i < layers; i++){
if(fList[i]) {
fList[i]->deleteList();
delete fList[i];
}
}
if(fList) delete[] fList;
#if DEBUG_LEVEL >2
gDebuger->leaveFunc("HMdcHitF::~HMdcHitF");
#endif
}
Bool_t HMdcHitF :: checkFit(HMdcHitAux* segment){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::checkFit");
#endif
Int_t sector = fLoc[0];
Int_t module = fLoc[1];
Float_t maxChi = parContainer->getMaxChi(sector, module);
Int_t minNumHits = parContainer->getMinNumHits(sector,module);
if(segment->getNumHits() < minNumHits || segment->getChi() > maxChi) return kFALSE;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::checkFit");
#endif
return kTRUE;
}
Bool_t HMdcHitF :: firstCheck(TMatrix param){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::firstcheck");
#endif
Int_t sector = fLoc[0];
Int_t module = fLoc[1];
Float_t maxSlopeX = parContainer->getMaxSlopeX(sector, module);
Float_t minSlopeX = parContainer->getMinSlopeX(sector, module);
Float_t maxSlopeY = parContainer->getMaxSlopeY(sector, module);
Float_t minSlopeY = parContainer->getMinSlopeY(sector, module);
Float_t zplane0 = parContainer->getZHit(sector,module);
if( (param(2,0) > maxSlopeX) || (param(2,0) < minSlopeX) ||
(param(3,0) > maxSlopeY) || (param(3,0) < minSlopeY)) return kFALSE;
if(target){
#if 0
Double_t x,y;
Double_t zdist = targetGeo->getZPos(sector, module)- zplane0;
x = param(0,0) + param(2,0) * zdist;
y = param(1,0) + param(3,0) * zdist;
if(fabs(x) > roadTargetX | fabs(y)> roadTargetY) return kFALSE;
#endif
Float_t ztarget = targetGeo->getZPos(sector,module);
Float_t u1[3];
u1[0] = param(2,0);
u1[1] = param(3,0);
u1[2] = 1;
Float_t ut[3];
ut[0] = 0;
if (module == 0) {
ut[1] = -0.697109;
ut[2] = 0.716965;
}
else if(module==1){
ut[1] = -0.808962;
ut[2] = 0.587861;
}
else{
Error("firstCheck","Unknown module %i",module);
return kFALSE;
}
Float_t points[3];
points[0] = param(0,0);
points[1] = param(1,0);
points[2] = zplane0 - ztarget;
Float_t normcross[3];
TMath::NormCross(u1,ut,normcross);
Float_t distance = 0;
for(Int_t l=0; l< 3;l++){
distance = distance + points[l]*normcross[l];
}
distance = fabs(distance);
Float_t maxDistance = parContainer->getDistZAxis(fLoc[0],fLoc[1]);
if(distance > maxDistance) return kFALSE;
}
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::firstcheck");
#endif
return kTRUE;
}
Int_t HMdcHitF :: find(Int_t listPlanes[6], Int_t nWantedHits, Int_t seg){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::find");
#endif
Int_t sector = fLoc[0];
Int_t module = fLoc[1];
Float_t road = parContainer->getRoad(sector,module);
TMatrix m(4,4);
TMatrix param(4,1);
TMatrix T(4,1);
Double_t z0 = parContainer->getZHit(sector,module);
HMdcCal3* listCal[6];
HMdcCal3* data;
HLocation locSeg;
HMdcHitAux *segment= new HMdcHitAux;
HMdcHitAux *segment2;
Double_t thetaint, thetaint2;
locSeg.setNIndex(3);
locSeg.setIndex(0,fLoc.getIndex(0));
locSeg.setIndex(1,fLoc.getIndex(1));
for(Int_t i=0; i<4; i++){
m(i,0) = geo->getSinn(sector, module, listPlanes[i]);
m(i,1) = -geo->getCosn(sector, module,listPlanes [i]);
m(i,2) = m(i,0) * (geo->getZPlane(sector, module,listPlanes[i])-z0);
m(i,3) = m(i,1) * (geo->getZPlane(sector, module,listPlanes[i])-z0);
}
m.Invert();
Int_t layers = geo->getNLayers();
HMdcCalLinkListIter **iter = new HMdcCalLinkListIter*[layers];
for (Int_t i=0;i<layers;i++) iter[i]=0;
HMdcCalLink* link;
for(Int_t i = 0; i< nWantedHits; i++){
iter[listPlanes[i]] = (HMdcCalLinkListIter* )fList[listPlanes[i]]->MakeIterator();
}
iter[listPlanes[0]]->Reset();
while((data=(HMdcCal3 *)iter[listPlanes[0]]->Next())!=NULL){
T(0,0) = data->getPos();
listCal[0]= data;
iter[listPlanes[1]]->Reset();
while((data=(HMdcCal3 *)iter[listPlanes[1]]->Next())!=NULL){
T(1,0) = data->getPos();
listCal[1] = data;
iter[listPlanes[2]]->Reset();
while((data=(HMdcCal3 *)iter[listPlanes[2]]->Next())!=NULL){
T(2,0) = data->getPos();
listCal[2] = data;
iter[listPlanes[3]]->Reset();
while((data=(HMdcCal3 *)iter[listPlanes[3]]->Next())!=NULL){
T(3,0) = data->getPos();
listCal[3] = data;
param.Mult(m,T);
if(firstCheck(param)){
if (nWantedHits > 4){
thetaint = geo->getSinn(sector, module,listPlanes[4])*
(param(0,0) + param(2,0)*
(geo->getZPlane(sector, module,listPlanes[4])-z0))-
geo->getCosn(sector, module,listPlanes[4])*
(param(1,0)+param(3,0)*
(geo->getZPlane(sector, module,listPlanes[4])-z0));
Int_t cellNum = calculateCell(sector,module,listPlanes[4],thetaint);
if(cellNum >=0){
Float_t pitch = geo->getPitch(sector,module,listPlanes[4]);
Int_t roadCells = Int_t(road/pitch)+1;
Int_t minCell1 = cellNum - roadCells;
if (minCell1 < 0) minCell1 = 0;
Int_t maxCell1 = cellNum + roadCells;
for(Int_t m1 = minCell1; m1 <= maxCell1; m1++){
HMdcCell* cell1 = fList[listPlanes[4]]->getCell(m1);
Int_t nHits1 = cell1->getNHits();
for(Int_t n1=0; n1<nHits1; n1++){
link = cell1->getHit(n1);
if(link!=NULL){
data = link->getData();
Double_t coord;
coord = data->getPos();
if(fabs(coord-thetaint) <= road){
listCal[4]= data;
if (nWantedHits > 5){
thetaint2 = geo->getSinn(sector,module,listPlanes[5])*
(param(0,0)+param(2,0)*
(geo->getZPlane(sector, module,
listPlanes[5])-z0))-
geo->getCosn(sector, module,listPlanes[5])*
(param(1,0)+ param(3,0)*
(geo->getZPlane(sector, module,
listPlanes[5])-z0));
cellNum = calculateCell(sector,module,listPlanes[5],thetaint2);
if(cellNum>=0){
pitch = geo->getPitch(sector,module,listPlanes[5]);
roadCells = Int_t(road/pitch)+1;
Int_t minCell2 = cellNum - roadCells;
if (minCell2 < 0) minCell2 = 0;
Int_t maxCell2 = cellNum + roadCells;
for(Int_t m2 = minCell2; m2 <= maxCell2; m2++){
HMdcCell* cell2 = fList[listPlanes[5]]->getCell(m2);
Int_t nHits2 = cell2->getNHits();
for(Int_t n2=0; n2<nHits2; n2++){
link = cell2->getHit(n2);
if(link!=NULL){
data = link->getData();
coord = data->getPos();
if(fabs(coord-thetaint2) <= road){
listCal[5] = data;
for(Int_t m=0; m< nWantedHits; m++){
segment->addCal(listCal[m],listPlanes[m]);
}
fit(segment);
if(checkFit(segment)){
if(!checkCommon(segment)){
locSeg.setIndex(2,seg);
segment2 = (HMdcHitAux*)mdcSegments->getSlot(locSeg);
if(segment2){
segment2 = new(segment2) HMdcHitAux;
*segment2 = *segment;
}
seg++;
segment->clear();
}
}
}
}
}
}
}
}
else{
for(Int_t m=0; m< nWantedHits; m++){
segment->addCal(listCal[m],listPlanes[m]);
}
fit(segment);
if(checkFit(segment)){
if(!checkCommon(segment)){
locSeg.setIndex(2,seg);
segment2 = (HMdcHitAux*)mdcSegments->getSlot(locSeg);
if(segment2){
segment2= new(segment2) HMdcHitAux;
*segment2 = *segment;
}
seg++;
segment->clear();
}
}
}
}
}
}
}
}
}
else{
for(Int_t i=0; i<nWantedHits; i++){
segment->addCal(listCal[i], listPlanes[i]);
}
segment->setPointX(param(0,0));
segment->setPointY(param(1,0));
segment->setPointZ(z0);
segment->setErrorX(0.);
segment->setErrorY(0.);
segment->setSlopeX(param(2,0));
segment->setSlopeY(param(3,0));
segment->setErrorSlopeX(0.);
segment->setErrorSlopeY(0.);
segment->setChi(0.);
locSeg.setIndex(2,seg);
segment2 = (HMdcHitAux*)mdcSegments->getSlot(locSeg);
if(segment2){
segment2= new(segment2) HMdcHitAux;
*segment2 = *segment;
}
seg++;
segment->clear();
#if DEBUG_LEVEL>3
HMdcHitAux* temp;
HMatrixCatIter* iter = NULL;
iter=(HMatrixCatIter *)mdcSegments->MakeIterator();
iter->Reset();
iter->gotoLocation(fLoc);
while((temp=(HMdcHitAux *)iter->Next())!=NULL) {
if(temp->isDel()==0) temp->print(full);
}
delete iter; iter = 0;
full << endl << endl;
full << endl << endl;
#endif
}
}
}
}
}
}
delete segment;
for(Int_t i=0; i< nWantedHits; i++){
delete iter[listPlanes[i]];
iter[listPlanes[i]] = 0;
}
delete[] iter;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::find");
#endif
return seg;
}
void HMdcHitF :: markHits(Int_t first, Int_t last){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::markHits");
#endif
HMdcHitAux* temp;
HMdcCal3* hit, *hit2;
HMdcCalLink* partner, *link;
HLocation loc;
loc.setNIndex(3);
loc.setIndex(0,fLoc.getIndex(0));
loc.setIndex(1,fLoc.getIndex(1));
loc.setIndex(2,first);
Int_t layers = ((HMatrixCategory*)mdcCal)->getSize(2);
for(; first <= last; first++){
loc.setIndex(2,first);
temp = (HMdcHitAux*)mdcSegments->getObject(loc);
if(temp!=NULL){
temp->markHits();
for(Int_t i=0; i<layers; i++){
hit = temp->getCal(i);
if(hit!=NULL){
link = hit->getLink();
if(link != NULL){
partner = hit->getLink()->getPartner();
if(partner != NULL){
hit2 = partner->getData();
if(hit2!=NULL) hit2->incrFlag();
}
}
fList[i]->remove(hit);
hit->setLink(NULL);
}
}
}
}
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::markHits");
#endif
}
Bool_t HMdcHitF:: checkCommon(HMdcHitAux* thissegment){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::checkCommon");
#endif
HMdcHitAux *temp;
Int_t numCommons = 0;
Int_t maxNumCommonHits = parContainer->getMaxNumCommonHits(fLoc[0],fLoc[1]);
HMatrixCatIter* iter = NULL;
iter=(HMatrixCatIter *)mdcSegments->MakeIterator();
iter->Reset();
iter->gotoLocation(fLoc);
while((temp=(HMdcHitAux *)iter->Next())!=NULL) {
if(temp->isDel()==0) {
if(temp->getNumHits() == thissegment->getNumHits()){
for(Int_t i=0; i< 6; i++){
if(temp->getCal(i)!=NULL && thissegment->getCal(i)!=NULL){
HMdcCalLink *link=thissegment->getCal(i)->getLink();
HMdcCalLink *tlink=temp->getCal(i)->getLink();
if (link && tlink) {
HMdcCalLink *partner=link->getPartner();
HMdcCalLink *tpartner=tlink->getPartner();
if (partner && tpartner) {
if ((temp->getCal(i)==thissegment->getCal(i) ||
temp->getCal(i)==thissegment->getCal(i)->getLink()->
getPartner()->getData() ||
temp->getCal(i)->getLink()->getPartner()->getData()==
thissegment->getCal(i) ||
temp->getCal(i)->getLink()->getPartner()->getData()==
thissegment->getCal(i)->getLink()->getPartner()->getData())){
numCommons++;
}
} else Warning("checkCommon","Partner not found");
} else Warning("checkCommom","Inconsistency in link list");
}
}
if(numCommons > maxNumCommonHits){
if(thissegment->getChi()< temp->getChi()){
temp->setDel();
}
else {
delete iter;
return kTRUE;
}
}
}
}
numCommons=0;
}
delete iter; iter = 0;
#if DEBUG_LEVEL >2
gDebuger->leaveFunc("HMdcHitF::checkCommon");
#endif
return kFALSE;
}
void HMdcHitF :: printSegments(ofstream& file){
HMdcHitAux* temp;
HMatrixCatIter* iter = NULL;
HLocation loc;
iter=(HMatrixCatIter *)mdcSegments->MakeIterator();
iter->Reset();
iter->gotoLocation(fLoc);
while((temp=(HMdcHitAux *)iter->Next())!=NULL) {
loc = iter->getLocation();
file << "pos" << '\t' << loc[0] << '\t' << loc[1] << '\t' <<loc[2] << endl;
temp->print(file);
}
delete iter;
iter = 0;
}
void HMdcHitF :: buildList(void){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::buildList");
#endif
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::buildList");
#endif
HMdcCal2* cal;
Int_t nlayers = geo->getNLayers();
Float_t pitch[nlayers];
Float_t centralWire[nlayers];
Int_t sector = fLoc[0];
Int_t module = fLoc[1];
Int_t layer, cell, nHits;
Int_t offset[nlayers];
for(Int_t i=0; i<nlayers; i++){
offset[i] = 0;
pitch[i] = geo ->getPitch(sector,module,i);
centralWire[i] = geo->getCentralWire(sector,module,i)-1;
}
HLocation loc;
loc.set(4,sector,module,0,0);
HMdcCal3 *hit1, *hit2;
Float_t pos1, pos2, coordwire;
HIterator* iter=0;
iter = (HIterator*)fCalCat->MakeIterator();
iter->Reset();
iter->gotoLocation(fLoc);
while((cal=(HMdcCal2*)iter->Next())!=NULL){
nHits = cal->getNHits();
if(nHits != 0){
layer = cal->getLayer();
cell = cal->getCell();
loc[2] = layer;
loc[3] = offset[layer];
hit1 = (HMdcCal3*)mdcCal->getSlot(loc);
offset[layer]++;
loc[3] = offset[layer];
hit2 = (HMdcCal3*)mdcCal->getSlot(loc);
offset[layer]++;
coordwire = (cell-centralWire[layer])*pitch[layer];
if(hit1 && hit2){
hit1 = new(hit1) HMdcCal3;
hit1->setAddress(sector,module,layer,cell,0);
hit1->setError(cal->getErrDist1());
if(layer==0 || layer == 4) coordwire = -coordwire;
pos1 = coordwire + cal->getDist1();
pos2 = coordwire - cal->getDist1();
hit1->setPos(pos1);
hit2 = new(hit2) HMdcCal3(hit2);
hit2->setAddress(sector,module,layer,cell,1);
hit2->setPos(pos2);
hit2->setError(cal->getErrDist1());
fList[layer]->addHits(hit1,hit2);
}
if(nHits<-1){
loc[3] = offset[layer];
hit1 = (HMdcCal3*)mdcCal->getSlot(loc);
offset[layer]++;
loc[3] = offset[layer];
hit2 = (HMdcCal3*)mdcCal->getSlot(loc);
offset[layer]++;
if(hit1 && hit2){
hit1 = new(hit1) HMdcCal3;
hit1->setAddress(sector,module,layer,cell,2);
hit1->setError(cal->getErrDist2());
pos1 = coordwire + cal->getDist2();
pos2 = coordwire - cal->getDist2();
hit1->setPos(pos1);
hit2 = new(hit2) HMdcCal3(hit2);
hit2->setAddress(sector,module,layer,cell,3);
hit2->setPos(pos2);
hit2->setError(cal->getErrDist2());
fList[layer]->addHits(hit1,hit2);
}
}
}
}
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::buildList");
#endif
}
Float_t HMdcHitF :: calculatePosPartner(HMdcCal3* hit1){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::calculatePosPartner");
#endif
Float_t coordwire, distwire;
Float_t posPartner;
Int_t layer = hit1->getLayer();
Int_t cell = hit1->getCell();
if(layer==0 || layer ==4) cell = -cell;
Float_t coord = hit1->getPos();
coordwire = geo->getFirstWirePos(fLoc[0],fLoc[1],layer) +
cell*geo->getPitch(fLoc[0],fLoc[1],layer);
distwire = coord - coordwire;
posPartner = coord - 2*distwire;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::calculatePosPartner");
#endif
return posPartner;
}
Int_t HMdcHitF :: execute(void){
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HMdcHitF::execute");
#endif
buildList();
#if DEBUG_LEVEL>3
ofstream file ("calhit.txt",ios::app);
file << "suceso " << fEventId << endl;
printCalHitCat(file);
file.close();
#endif
Int_t nWantedHits, nMdcPlanes, nPlanesScanned;
Int_t nPossibleFirstPlanes, nMissingHitsTolerated;
Int_t minNumHits = parContainer->getMinNumHits(fLoc[0],fLoc[1]);
nMdcPlanes = geo->getNLayers();
Int_t segPos =0, posIni=0;
Int_t listPlanes[6];
HMdcHitsFilter* filtering = new HMdcHitsFilter;
for(nWantedHits=nMdcPlanes; nWantedHits>=minNumHits;nWantedHits--){
for(nPlanesScanned=nMdcPlanes; nPlanesScanned>=nWantedHits;
nPlanesScanned--){
nMissingHitsTolerated = nMdcPlanes - nWantedHits;
nPossibleFirstPlanes = nMdcPlanes - nPlanesScanned +1;
Int_t firstScannedPlane = 1;
Int_t lastFirstScannedPlane = firstScannedPlane + nPossibleFirstPlanes-1;
for(; firstScannedPlane<=lastFirstScannedPlane; firstScannedPlane++){
Int_t lastScannedPlane= firstScannedPlane+nPlanesScanned-1;
listPlanes[0] = firstScannedPlane-1;
listPlanes[3] = lastScannedPlane-1;
for(Int_t secondScannedPlane = firstScannedPlane+1;
secondScannedPlane<= firstScannedPlane+nMissingHitsTolerated+1;
secondScannedPlane++){
listPlanes[1] = secondScannedPlane-1;
for(Int_t thirdScannedPlane = lastScannedPlane-1;
thirdScannedPlane-secondScannedPlane>= 3-nMissingHitsTolerated;
thirdScannedPlane--){
listPlanes[2] = thirdScannedPlane-1;
Int_t nIntermPlanes = 2 - nMissingHitsTolerated;
Int_t firstIntermPlane = secondScannedPlane +1;
Int_t lastFirstIntermPlane = thirdScannedPlane-nIntermPlanes;
if(nWantedHits==4){
segPos=find(listPlanes, nWantedHits,segPos);
}
else{
for (Int_t intermPlane=firstIntermPlane;
intermPlane <= lastFirstIntermPlane;
intermPlane++){
for (Int_t index=0; index < nIntermPlanes; index++){
listPlanes[4+index] = intermPlane+index-1;
}
segPos=find(listPlanes, nWantedHits,segPos);
}
}
}
}
}
}
mdcSegments->filter(*filtering);
markHits(posIni,segPos-1);
posIni = segPos;
}
delete filtering;
for (Int_t i=0; i< 6; i++){
if (fList[i]) fList[i]->deleteList();
}
if(slopeCorrection) makeSlopeCorrection();
fillHitCat();
fEventId++;
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HMdcHitF::execute");
#endif
return 0;
}
Int_t HMdcHitF :: calculateCell(Int_t sector, Int_t module, Int_t layer, Float_t coordinate){
Float_t distance;
Float_t diff;
Int_t cell;
Float_t pitch;
Float_t offset;
Int_t maxCells;
pitch = geo->getPitch(sector,module,layer);
offset = geo->getFirstWirePos(sector,module,layer);
maxCells = geo->getMaxNumCells(sector,module,layer);
diff = coordinate-offset;
distance = fabs(diff)+pitch/2.;
if(layer==0 || layer==4) diff=-diff;
cell = Int_t(distance/pitch);
if(diff < 0 || cell > maxCells) return -1;
else return cell;
}
void HMdcHitF :: makeSlopeCorrection(void){
#if DEBUG_LEVEL> 2
gDebuger->enterFunc("HMdcHitF::makeSlopeCorrection");
#endif
Int_t sector = fLoc[0];
Int_t module = fLoc[1];
Int_t nlayers = geo->getNLayers();
Float_t pitch[nlayers];
Float_t wiredist,wirepos;
Float_t centralWire[nlayers];
for(Int_t layer=0;layer<nlayers;layer++){
pitch[layer] = geo->getPitch(sector, module, layer);
centralWire[layer] = geo->getCentralWire(sector, module, layer) -1;
}
HMdcCal3 *hit;
HMdcHitAux* seg;
Float_t pos[nlayers];
Float_t newpos, newSlopeY, correction;
Float_t mycosn[6] = {0.76604444,0.93969262,1.,1.,0.93969262,0.76604444};
Float_t mysinn[6] = {0.64278761,-0.34202014,0.,0.,0.34202014,-0.64278761};
Float_t chi;
HMatrixCatIter* iter=0;
iter = (HMatrixCatIter *)mdcSegments->MakeIterator();
iter->Reset();
iter->gotoLocation(fLoc);
while((seg=(HMdcHitAux*)iter->Next())!=NULL){
chi = seg->getChi();
for(Int_t i=0; i<nlayers; i++){
hit = seg->getCal(i);
if(hit!=NULL){
pos[i] = hit->getPos();
wirepos = (hit->getCell() - centralWire[i]) * pitch[i];
if(i == 0 || i == 4) wirepos = -wirepos;
wiredist = pos[i] - wirepos;
newSlopeY = seg->getSlopeY() * mycosn[i] - seg->getSlopeX() * mysinn[i];
correction = (wiredist/cos(atan(newSlopeY)))-wiredist;
newpos = pos[i] + correction;
hit->setPos(newpos);
}
}
fit(seg);
if(seg->getChi()>chi){
Float_t bestchi = seg->getChi();
HMdcHitAux* seg2 = new HMdcHitAux();
*seg2 = *seg;
for(Int_t i=0;i<nlayers;i++){
hit = seg->getCal(i);
if(hit!=NULL){
pos[i] = hit->getPos();
wirepos = (hit->getCell() - centralWire[i]) * pitch[i];
if(i == 0 || i == 4) wirepos = -wirepos;
wiredist = pos[i] - wirepos;
newpos = wirepos - wiredist;
hit->setPos(newpos);
fit(seg);
if(seg->getChi()>chi){
if(seg->getChi()<bestchi){
*seg2 = *seg;
bestchi = seg2->getChi();
}
else hit->setPos(pos[i]);
if(i==5) {
*seg = *seg2;
fit(seg);
}
}
else {
break;
}
}
}
}
}
#if DEBUG_LEVEL> 2
gDebuger->leaveFunc("HMdcHitF::makeSlopeCorrection");
#endif
}
void HMdcHitF :: fillHitCat(void) {
#if DEBUG_LEVEL > 2
gDebuger->enterFunc("HMdcHitF::fillHitCat");
#endif
Int_t nlayers = geo->getNLayers();
HMdcCal3 *calhit;
HMdcHitAux* seg;
HMdcHit *hit;
HLocation loc;
loc.set(3,fLoc[0],fLoc[1],0);
Int_t offset = 0;
HMatrixCatIter* iter=0;
iter = (HMatrixCatIter *)mdcSegments->MakeIterator();
iter->Reset();
iter->gotoLocation(fLoc);
while((seg=(HMdcHitAux*)iter->Next())!=NULL){
loc[2] = offset;
hit = (HMdcHit*)fHitCat->getSlot(loc);
if (hit){
hit = new (hit) HMdcHit;
hit->setX( seg->getPointX(), seg->getXError());
hit->setY( seg->getPointY(), seg->getYError());
Float_t xslope = seg->getSlopeX();
Float_t yslope = seg->getSlopeY();
Float_t det = sqrt(xslope*xslope + yslope*yslope + 1);
Float_t xdir = xslope/det;
Float_t ydir = yslope/det;
hit->setXDir(xdir, seg->getErrorSlope1());
hit->setYDir(ydir, seg->getErrorSlope2());
hit->setChi2(seg->getChi());
hit->setFlag(seg->getIsUsed());
hit->setSecMod(fLoc[0], fLoc[1]);
Int_t id=0;
Int_t hitNum = 0;
for(Int_t i=0; i<nlayers; i++){
calhit = seg->getCal(i);
if(calhit!=NULL){
hitNum = calhit->getHitNumber();
if (hitNum == 0 || hitNum == 1) id = 1;
else if (hitNum == 2 || hitNum == 3) id = 2;
else id = 0;
hit->setSignId(i, calhit->getCell(), id);
}
}
offset++;
}
}
delete iter;
#if DEBUG_LEVEL > 2
gDebuger->leaveFunc("HMdcHitF::fillHitCat");
#endif
}
void HMdcHitF :: printCalHitCat(ofstream &file){
#if DEBUG_LEVEL > 2
gDebuger->enterFunc("HMdcHitF::printCalHitCat");
#endif
HMdcCal3* calhit;
HIterator* iter = NULL;
iter = (HMatrixCatIter *)mdcCal->MakeIterator();
iter->Reset();
while((calhit=(HMdcCal3 *)iter->Next())!=NULL ){
calhit->print(file);
}
#if DEBUG_LEVEL > 2
gDebuger->enterFunc("HMdcHitF::printCalHitCat");
#endif
}
Last change: Sat May 22 13:02:20 2010
Last generated: 2010-05-22 13:02
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.