#include "TRandom.h"
#include "hrpcdetector.h"
#include "hrpcclusterf.h"
#include "rpcdef.h"
#include "hrpchit.h"
#include "hrpchitsim.h"
#include "hrpccluster.h"
#include "hrpcclustersim.h"
#include "hgeantrpc.h"
#include "hrpcgeomcellpar.h"
#include "hrpcclusfpar.h"
#include "hdebug.h"
#include "hades.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "hdetector.h"
#include "hevent.h"
#include "hgeantkine.h"
#include "hcategory.h"
#include <iostream>
#include <iomanip>
#include <time.h>
using namespace std;
void HRpcClusterF::initParContainer() {
pGeomCellPar = (HRpcGeomCellPar*)(gHades->getRuntimeDb()->getContainer("RpcGeomCellPar"));
pClusFPar = (HRpcClusFPar*)(gHades->getRuntimeDb()->getContainer("RpcClusFPar"));
}
HRpcClusterF::HRpcClusterF(void)
{
pClusterCat = NULL;
pHitCat = NULL;
pGeomCellPar = NULL;
pClusFPar = NULL;
doSingleCluster = kFALSE;
}
HRpcClusterF::HRpcClusterF(const Text_t *name, const Text_t *title) :
HReconstructor(name,title)
{
pClusterCat = NULL;
pHitCat = NULL;
pGeomCellPar = NULL;
pClusFPar = NULL;
doSingleCluster = kFALSE;
}
HRpcClusterF::~HRpcClusterF(void)
{
}
Bool_t HRpcClusterF::init(void)
{
initParContainer();
HRpcDetector *rpc;
rpc=(HRpcDetector *)gHades->getSetup()->getDetector("Rpc");
maxcells=rpc->getMaxCells();
if (!pGeomCellPar){
Error("init","No RpcGeomCellPar Parameters");
return kFALSE;
}
if (!pClusFPar){
Error("init","No RpcClusFPar Parameters");
return kFALSE;
}
if(!rpc){
Error("HRpcClusterF::init()","No Rpc Detector");
return kFALSE;
}
pHitCat=gHades->getCurrentEvent()->getCategory(catRpcHit);
if (!pHitCat) {
Error("HRpcClusterF::init()","No RpcHit Category");
return kFALSE;
}
HCategory* catKin=gHades->getCurrentEvent()->getCategory(catGeantKine);
if(catKin) simulation=kTRUE;
else simulation=kFALSE;
pClusterCat=gHades->getCurrentEvent()->getCategory(catRpcCluster);
if (!pClusterCat) {
if (simulation) pClusterCat=rpc->buildMatrixCategory("HRpcClusterSim",0.5);
else pClusterCat=rpc->buildMatrixCategory("HRpcCluster",0.5);
gHades->getCurrentEvent()->addCategory(catRpcCluster,pClusterCat,"Rpc");
}
if(simulation) {
pGeantRpcCat = gHades->getCurrentEvent()->getCategory(catRpcGeantRaw);
if (!pGeantRpcCat) {
Error("HRpcClusterF::init()","No HGeant RPC Category");
return kFALSE;
}
}
loc.set(3,0,0,0);
loc_cluster.set(1,0);
return kTRUE;
}
Int_t HRpcClusterF::execute(void)
{
HRpcHitSim *pHit = NULL;
HRpcClusterSim *pClusterSim = NULL;
HRpcCluster *pCluster = NULL;
Int_t cluster_type = -1;
Int_t sign = 0;
Float_t tof, charge, xmod, ymod, zmod, xsec, ysec, zsec;
Float_t xlab, ylab, zlab, theta, phi;
Float_t tof_neigb, charge_neigb, xmod_neigb, ymod_neigb, xsec_neigb, ysec_neigb, zsec_neigb;
Float_t xlab_neigb, ylab_neigb, zlab_neigb, theta_neigb, phi_neigb;
Float_t tof_neigb2;
Float_t sigma_X, sigma_Y, sigma_Z, sigma_T;
Float_t tof_cut, x_cut, DPlanes, tof_next;
Float_t v_light=300;
Int_t Track1L=-999, Track1R=-999, Track2L=-999, Track2R=-999, Track=-999;
Int_t Ref1L=-999, Ref1R=-999, Ref2L=-999, Ref2R=-999;
Bool_t LisAtBox1=kFALSE, LisAtBox2=kFALSE, RisAtBox1=kFALSE, RisAtBox2=kFALSE, IsAtBox=kFALSE;
Float_t NsigmasT, NsigmasX;
Int_t cell1=-1, column1=-1, sector1=-1, cell2=-1, column2=-1, sector2=-1;
Bool_t inFlag1=kFALSE, inFlag2=kFALSE;
Int_t weight[4], array[4], ref_array[4];
Bool_t isAtBox_array[4];
Bool_t BC1IsNotIn=kFALSE;
for (Int_t i=0;i<4;i++) weight[i] = 0;
for (Int_t i=0;i<4;i++) array[i] = -999;
for (Int_t i=0;i<4;i++) ref_array[i] = -999;
for (Int_t i=0;i<4;i++) isAtBox_array[i] = kFALSE;
DPlanes = pGeomCellPar->getDPlanes();
if(doSingleCluster && pHitCat && pClusterCat){
TIterator *pRpcHitIter = pHitCat->MakeIterator("native");
if (pRpcHitIter) {
HRpcHit* rpc_hit=0;
pRpcHitIter->Reset();
while ((rpc_hit = (HRpcHit*)pRpcHitIter->Next())) {
loc[0] = rpc_hit->getSector();
loc[1] = rpc_hit->getColumn();
loc[2] = rpc_hit->getCell();
Int_t index_cluster;
loc_cluster[0] = rpc_hit->getSector();
pClusterSim = (HRpcClusterSim*)pClusterCat->getNewSlot(loc_cluster,&index_cluster);
if ( !pClusterSim ) {
Error ("execute()", "Can't get slot in cluster finder: sec %i", loc_cluster[0]);
return EXIT_FAILURE;
}
if (simulation)pClusterSim = new(pClusterSim) HRpcClusterSim;
else pClusterSim = (HRpcClusterSim*) new(pClusterSim) HRpcCluster;
pClusterSim->setClusterType (1);
pClusterSim->setInsideCellFlag(rpc_hit->getInsideCellFlag());
pClusterSim->setDetID1 (rpc_hit->getSector(),rpc_hit->getColumn(),rpc_hit->getCell());
pClusterSim->setCluster(rpc_hit->getTof() ,rpc_hit->getCharge(),rpc_hit->getXMod(),rpc_hit->getYMod(),rpc_hit->getZMod());
pClusterSim->setRSEC (rpc_hit->getXSec() ,rpc_hit->getYSec() ,rpc_hit->getZSec());
pClusterSim->setRMS (rpc_hit->getTOFRMS(),rpc_hit->getXRMS() ,rpc_hit->getYRMS(),rpc_hit->getZRMS());
pClusterSim->setAddress(rpc_hit->getSector(),index_cluster);
pClusterSim->setTheta (rpc_hit->getTheta());
pClusterSim->setPhi (rpc_hit->getPhi());
rpc_hit ->getXYZLab (xlab,ylab,zlab);
pClusterSim->setXYZLab (xlab,ylab,zlab);
if (simulation) {
for (Int_t t=0;t<4;t++) weight[t] = 0;
for (Int_t t=0;t<4;t++) array[t] = -999;
for (Int_t t=0;t<4;t++) ref_array[t] = -999;
for (Int_t t=0;t<4;t++) isAtBox_array[t] = kFALSE;
array [0] = pHit->getTrackL();
array [1] = pHit->getTrackR();
ref_array [0] = pHit->getRefL();
ref_array [1] = pHit->getRefR();
isAtBox_array[0] = pHit->getLisAtBox();
isAtBox_array[1] = pHit->getRisAtBox();
for(Int_t m=0;m<4;m++){
for(Int_t n=0;n<4;n++){
if(array[m]==array[n]) weight[m]++;
}
}
Track = array [0];
IsAtBox = isAtBox_array[0];
for(Int_t m=0;m<3;m++) {
if (weight[m+1]>weight[m]) {
Track=array[m+1];
IsAtBox = isAtBox_array[m];}
}
pClusterSim->setTrack(Track);
pClusterSim->setIsAtBox(IsAtBox);
for(Int_t m=0;m<4;m++) {
pClusterSim->setTrackList (m,array[m]);
pClusterSim->setRefList (m,ref_array[m]);
pClusterSim->setIsAtBoxList(m,isAtBox_array[m]);
}
pClusterSim->setNTracksAtBox (howManyTracks (pClusterSim));
pClusterSim->setNTracksAtCells(howManyTracksAtCells(pClusterSim));
}
}
}
delete pRpcHitIter ;
return EXIT_SUCCESS;
}
for(Int_t i=0;i<6;i++) {
NsigmasX = pClusFPar->getNsigmasX();
NsigmasT = pClusFPar->getNsigmasT();
for(Int_t j=0;j<3;j++) {
for(Int_t k=0;k<(maxcells*2-1);k++){
for (Int_t t=0;t<4;t++) weight[t] = 0;
for (Int_t t=0;t<4;t++) array[t] = -999;
for (Int_t t=0;t<4;t++) ref_array[t] = -999;
for (Int_t t=0;t<4;t++) isAtBox_array[t] = kFALSE;
cell1 = column1 = sector1 = -1;
cell2 = column2 = sector2 = -1;
loc[0] = i;
loc[1] = 2*j+1-k%2;
loc[2] = k/2;
if(BC1IsNotIn && j==1) loc[1] = 2+k%2;
pHit = (HRpcHitSim*)pHitCat->getObject(loc);
if(!pHit) continue;
loc_cluster[0]=i;
Int_t index_cluster;
sigma_X= pHit->getXRMS();
sigma_Y= pHit->getYRMS();
sigma_Z= pHit->getZRMS();
sigma_T= pHit->getTOFRMS();
tof_cut = sigma_T*NsigmasT;
x_cut = sigma_X*NsigmasX;
if (simulation) {
pClusterSim = (HRpcClusterSim*)pClusterCat->getNewSlot(loc_cluster,&index_cluster);
if ( !pClusterSim ) {
Error ("execute()", "Can't get slot in cluster finder: sec %i", loc_cluster[0]);
return EXIT_FAILURE;
}
pClusterSim = new(pClusterSim) HRpcClusterSim;
}
else {
pCluster = (HRpcCluster*)pClusterCat->getNewSlot(loc_cluster,&index_cluster);
if ( !pCluster ) {
Error ("execute()", "Can't get slot in cluster finder: sec %i", loc_cluster[0]);
return EXIT_FAILURE;
}
pCluster = new(pCluster) HRpcCluster;
}
tof = pHit->getTof();
charge = pHit->getCharge();
xmod = pHit->getXMod();
ymod = pHit->getYMod();
zmod = DPlanes*(0.5-k%2);
if(BC1IsNotIn && j==1) zmod = DPlanes*(k%2-0.5);
xsec = pHit->getXSec();
ysec = pHit->getYSec();
zsec = pHit->getZSec();
theta = pHit->getTheta();
phi = pHit->getPhi();
sector1 = pHit->getSector();
column1 = pHit->getColumn();
cell1 = pHit->getCell();
inFlag1 = pHit->getInsideCellFlag();
pHit->getXYZLab(xlab,ylab,zlab);
cluster_type=1;
if (simulation){
Ref1L = pHit->getRefL();
Ref1R = pHit->getRefR();
Track1L = pHit->getTrackL();
Track1R = pHit->getTrackR();
LisAtBox1 = pHit->getLisAtBox();
RisAtBox1 = pHit->getRisAtBox();
}
loc[1] = 2*j + 1 - (k+1)%2;
loc[2] = (k+1)/2;
if(BC1IsNotIn && j==1) loc[1] = 2 + (k+1)%2;
pHit = (HRpcHitSim*)pHitCat->getObject(loc);
if (pHit) {
tof_neigb = pHit->getTof();
charge_neigb = pHit->getCharge();
xmod_neigb = pHit->getXMod();
ymod_neigb = pHit->getYMod();
xsec_neigb = pHit->getXSec();
ysec_neigb = pHit->getYSec();
zsec_neigb = pHit->getZSec();
theta_neigb = pHit->getTheta();
phi_neigb = pHit->getPhi();
sector2 = pHit->getSector();
column2 = pHit->getColumn();
cell2 = pHit->getCell();
inFlag2 = pHit->getInsideCellFlag();
pHit->getXYZLab(xlab_neigb,ylab_neigb,zlab_neigb);
sign = (2*loc[2]%2-1);
tof_next = sign*DPlanes/v_light;
if((fabs(xmod-xmod_neigb)<x_cut) && fabs(tof+tof_next-tof_neigb)<tof_cut){
tof = (tof + tof_neigb)/2;
charge = (charge + charge_neigb)/2;
xmod = (xmod + xmod_neigb)/2;
ymod = (ymod + ymod_neigb)/2;
zmod = 0.001;
xsec = (xsec + xsec_neigb)/2;
ysec = (ysec + ysec_neigb)/2;
zsec = (zsec + zsec_neigb)/2;
xlab = (xlab + xlab_neigb)/2;
ylab = (ylab + ylab_neigb)/2;
zlab = (zlab + zlab_neigb)/2;
theta = (theta + theta_neigb)/2;
phi = (phi + phi_neigb)/2;
cluster_type=2;
if (simulation) {
Ref2L = pHit->getRefL();
Ref2R = pHit->getRefR();
Track2L = pHit->getTrackL();
Track2R = pHit->getTrackR();
LisAtBox2 = pHit->getLisAtBox();
RisAtBox2 = pHit->getRisAtBox();
}
loc[1] = 2*j + 1 - (k+2)%2;
loc[2] = (k+2)/2;
if(BC1IsNotIn && j==1) loc[1] = 2 + (k+2)%2;
pHit = (HRpcHitSim*)pHitCat->getObject(loc);
if (pHit) {
tof_neigb2 = pHit->getTof();
sign = (2*loc[2]%2-1);
tof_next = sign*DPlanes/v_light;
if((fabs(xmod-xmod_neigb)>=x_cut) || fabs(tof+tof_next-tof_neigb2)>=tof_cut) k=k+1;
}
else k=k+1;
}
}
if (simulation) {
pClusterSim->setClusterType(cluster_type);
pClusterSim->setCluster(tof, charge, xmod, ymod, zmod);
pClusterSim->setRSEC(xsec,ysec,zsec);
pClusterSim->setXYZLab(xlab,ylab,zlab);
pClusterSim->setTheta(theta);
pClusterSim->setPhi(phi);
pClusterSim->setAddress(loc_cluster[0], index_cluster);
if(cluster_type==1) {
pClusterSim->setDetID1(sector1, column1, cell1);
array[0]=Track1L; array[1]=Track1R;
ref_array[0]=Ref1L; ref_array[1]=Ref1R;
isAtBox_array[0]=LisAtBox1; isAtBox_array[1]=RisAtBox1;
pClusterSim->setRMS(sigma_T, sigma_X, sigma_Y, sigma_Z);
}
else{
if(column1==0 || column1==2 || column1==4) {
pClusterSim->setDetID1(sector1, column1, cell1);
pClusterSim->setDetID2(sector2, column2, cell2);
array[0]=Track1L; array[1]=Track1R; array[2]=Track2L; array[3]=Track2R;
ref_array[0]=Ref1L; ref_array[1]=Ref1R; ref_array[2]=Ref2L; ref_array[3]=Ref2R;
isAtBox_array[0]=LisAtBox1; isAtBox_array[1]=RisAtBox1; isAtBox_array[2]=LisAtBox2; isAtBox_array[3]=RisAtBox2;
}
if(column1==1 || column1==3 || column1==5) {
pClusterSim->setDetID1(sector2, column2, cell2);
pClusterSim->setDetID2(sector1, column1, cell1);
array[0]=Track2L; array[1]=Track2R; array[2]=Track1L; array[3]=Track1R;
ref_array[0]=Ref2L; ref_array[1]=Ref2R; ref_array[2]=Ref1L; ref_array[3]=Ref1R;
isAtBox_array[0]=LisAtBox2; isAtBox_array[1]=RisAtBox2; isAtBox_array[2]=LisAtBox1; isAtBox_array[3]=RisAtBox1;
}
pClusterSim->setRMS(sigma_T/sqrt(2.), sigma_X/sqrt(2.), sigma_Y/sqrt(2.), sigma_Z/sqrt(2.));
}
for(Int_t m=0;m<4;m++){
for(Int_t n=0;n<4;n++){
if(array[m]==array[n]) weight[m]++;
}
}
Track = array[0];
IsAtBox = isAtBox_array[0];
for(Int_t m=0;m<3;m++) {
if (weight[m+1]>weight[m]) {
Track=array[m+1];
IsAtBox = isAtBox_array[m];}
}
pClusterSim->setTrack(Track);
pClusterSim->setIsAtBox(IsAtBox);
for(Int_t m=0;m<4;m++) {
pClusterSim->setTrackList(m,array[m]);
pClusterSim->setRefList(m,ref_array[m]);
pClusterSim->setIsAtBoxList(m,isAtBox_array[m]);
}
pClusterSim->setNTracksAtBox(howManyTracks(pClusterSim));
pClusterSim->setNTracksAtCells(howManyTracksAtCells(pClusterSim));
}
else {
pCluster->setClusterType(cluster_type);
pCluster->setCluster(tof, charge, xmod, ymod, zmod);
pCluster->setRSEC(xsec,ysec,zsec);
pCluster->setXYZLab(xlab,ylab,zlab);
pCluster->setTheta(theta);
pCluster->setPhi(phi);
pCluster->setAddress(loc_cluster[0], index_cluster);
if(cluster_type==1) {
pCluster->setDetID1(sector1, column1, cell1);
pCluster->setRMS(sigma_T, sigma_X, sigma_Y, sigma_Z);
if(inFlag1) pCluster->setInsideCellFlag(1);
}
else{
if(column1==0 || column1==2 || column1==4) {
pCluster->setDetID1(sector1, column1, cell1);
pCluster->setDetID2(sector2, column2, cell2);
}
if(column1==1 || column1==3 || column1==5) {
pCluster->setDetID1(sector2, column2, cell2);
pCluster->setDetID2(sector1, column1, cell1);
}
pCluster->setRMS(sigma_T/sqrt(2.), sigma_X/sqrt(2.), sigma_Y/sqrt(2.), sigma_Z/sqrt(2.));
if(inFlag1 && inFlag2) {
pCluster->setInsideCellFlag(1);
} else {
pCluster->setInsideCellFlag(0);
if(inFlag1) pCluster->setInsideCellFlag(2);
if(inFlag2) pCluster->setInsideCellFlag(3);
}
}
}
}
}
}
return EXIT_SUCCESS;
}
Int_t HRpcClusterF::howManyTracks(HRpcClusterSim* pClusterSim) {
HRpcHitSim *pHit1 = 0;
HRpcHitSim *pHit2 = 0;
HLocation loc1;
HLocation loc2;
loc1.set(3,0,0,0);
loc2.set(3,0,0,0);
Int_t TrackL1[10], TrackR1[10];
Int_t TrackL2[10], TrackR2[10];
Int_t TrackTot[40], TrackDifferent[40];
loc1[0] = pClusterSim->getSector1();
loc1[1] = pClusterSim->getColumn1();
loc1[2] = pClusterSim->getCell1();
pHit1 = (HRpcHitSim*)pHitCat->getObject(loc1);
pHit1->getTrackLArray(TrackL1);
pHit1->getTrackRArray(TrackR1);
for(Int_t i=0;i<10;i++) TrackL2[i] = -999;
for(Int_t i=0;i<10;i++) TrackR2[i] = -999;
if((pClusterSim->getClusterType())==2) {
loc2[0] = pClusterSim->getSector2();
loc2[1] = pClusterSim->getColumn2();
loc2[2] = pClusterSim->getCell2();
pHit2 = (HRpcHitSim*)pHitCat->getObject(loc2);
pHit2->getTrackLArray(TrackL2);
pHit2->getTrackRArray(TrackR2);
}
for(Int_t i=0;i<40;i++) TrackDifferent[i]=-999;
for(Int_t i=0;i<10;i++) TrackTot[i] = TrackL1[i];
for(Int_t i=0;i<10;i++) TrackTot[i+10] = TrackR1[i];
for(Int_t i=0;i<10;i++) TrackTot[i+20] = TrackL2[i];
for(Int_t i=0;i<10;i++) TrackTot[i+30] = TrackR2[i];
Int_t nTracks=0;
for(Int_t i=0;i<40;i++) {
if(TrackTot[i]<0) continue;
Int_t count=0;
for(Int_t j=0;j<40;j++) if(TrackDifferent[j]==TrackTot[i]) count++;
if(count==0) {
TrackDifferent[nTracks] = TrackTot[i];
nTracks++;
}
}
return nTracks;
}
Int_t HRpcClusterF::howManyTracksAtCells(HRpcClusterSim* pClusterSim)
{
HRpcHitSim *pHit1 = 0;
HRpcHitSim *pHit2 = 0;
HLocation loc1;
HLocation loc2;
loc1.set(3,0,0,0);
loc2.set(3,0,0,0);
Int_t TrackLDgtr1[10], TrackRDgtr1[10];
Int_t TrackLDgtr2[10], TrackRDgtr2[10];
Int_t TrackTot[40], TrackDifferent[40];
loc1[0] = pClusterSim->getSector1();
loc1[1] = pClusterSim->getColumn1();
loc1[2] = pClusterSim->getCell1();
pHit1 = (HRpcHitSim*)pHitCat->getObject(loc1);
pHit1->getTrackLDgtrArray(TrackLDgtr1);
pHit1->getTrackRDgtrArray(TrackRDgtr1);
for(Int_t i=0;i<10;i++) TrackLDgtr2[i] = -999;
for(Int_t i=0;i<10;i++) TrackRDgtr2[i] = -999;
if((pClusterSim->getClusterType())==2) {
loc2[0] = pClusterSim->getSector2();
loc2[1] = pClusterSim->getColumn2();
loc2[2] = pClusterSim->getCell2();
pHit2 = (HRpcHitSim*)pHitCat->getObject(loc2);
pHit2->getTrackLDgtrArray(TrackLDgtr2);
pHit2->getTrackRDgtrArray(TrackRDgtr2);
}
for(Int_t i=0;i<40;i++) TrackDifferent[i]=-999;
for(Int_t i=0;i<10;i++) TrackTot[i] = TrackLDgtr1[i];
for(Int_t i=0;i<10;i++) TrackTot[i+10] = TrackRDgtr1[i];
for(Int_t i=0;i<10;i++) TrackTot[i+20] = TrackLDgtr2[i];
for(Int_t i=0;i<10;i++) TrackTot[i+30] = TrackRDgtr2[i];
Int_t nTracks=0;
for(Int_t i=0;i<40;i++) {
if(TrackTot[i]<0) continue;
Int_t count=0;
for(Int_t j=0;j<40;j++) if(TrackDifferent[j]==TrackTot[i]) count++;
if(count==0) {
TrackDifferent[nTracks] = TrackTot[i];
nTracks++;
}
}
return nTracks;
}
ClassImp(HRpcClusterF)