496 vector<uint8_t> macs;
497 map< uint8_t, vector< pair<int,float> > > pesmap;
516 vector<infoA> informationA;
517 vector<infoA> informationB;
525 double hitpoint[3], hitpointerr[3];
526 TVector3 hitpos (0.,0.,0.);
529 float petot = 0., pemax = 0., pex=0., pey=0.;
530 int adsid_max = -1, nabove=0, nx=0, ny=0, nz = 0, ntrig = 0;
534 vector<uint64_t> ttrigs;
535 vector<uint64_t> t1trigs;
536 vector<TVector3> tpos;
537 double zmin=DBL_MAX, zmax = -DBL_MAX;
538 double ymin=DBL_MAX, ymax = -DBL_MAX;
539 double xmin=DBL_MAX, xmax = -DBL_MAX;
540 std::vector<int> layID;
541 std::vector<int> febA;
542 std::vector<int> febB;
545 uint64_t southt0_v = -999, southt0_h =-999;
548 for(
auto const& data : coinData) {
551 febA.push_back(data->fMac5);
553 febB.push_back(data->fMac5);
558 std ::cout <<
"line 451: size of febA: \t" << (int)febA.size()
559 <<
" size of febB: " << (int)febB.size() <<
'\n';
563 for(
auto const& data : coinData) {
566 macs.push_back(data->fMac5);
571 layID.push_back(layer);
573 auto idx = &data - coinData.data();
574 if(
fVerbose)
std :: cout <<
"index: " << idx <<
" , mac5: " << (int)macs.back() <<
" ,adid: " << adid<<
'\n';
577 if ((
int)febA.size() == 0 or (
int)febB.
size() == 0)
continue;
581 if (febA.size() < febB.size()) size = febB.size();
582 else size = febA.size();
585 for (
int i = 0; i < size; i++){
586 if (macs.back() == (int)febA[i]) {
589 for(
int chan=0; chan<32; chan++) {
591 float pe = (data->fAdc[chan]-chg_cal.second)/chg_cal.first;
597 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"\nfebA (mac5, channel, gain, pedestal, adc, pe) = (" << (int)data->fMac5 <<
", " << chan <<
", "
598 << chg_cal.first <<
", " << chg_cal.second <<
"," << data->fAdc[chan] <<
"," << pe <<
")\n";
602 informationA.push_back ({macs.back(), chan, data->fTs0, postmp, adsidA});
606 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"looking for mac " << (int)macs.back()
607 <<
" and matching with febA : " << (int)febA[i]<<
'\n';
608 }
else if ( macs.back() == (int)febB[i]){
611 for(
int chan=0; chan<32; chan++) {
613 float pe = (data->fAdc[chan]-chg_cal.second)/chg_cal.first;
618 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"\nfebB (mac5, channel, gain, pedestal, adc, pe) = (" << (int)data->fMac5 <<
", " << chan <<
", "
619 << chg_cal.first <<
", " << chg_cal.second <<
"," << data->fAdc[chan] <<
"," << pe <<
")\n";
623 informationB.push_back ({macs.back(), chan, data->fTs0, postmp, adsidB});
626 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"else if looking for mac "<< (int)macs.back()
627 <<
" and matching with febB : " << (int)febB[i]<<
'\n';
632 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"In CRTHitRecoAlg::MakeSideHit 1st feb is " << (int)
fCrtutils->
ADToMac(adid).first
634 <<
", time "<< data->fTs0
635 <<
" with module number " << adid <<
", no. of FEB " <<
'\n';
638 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"In CRTHitRecoAlg::MakeSideHit functions mac is " << (int)macs.back()
639 <<
" with module number " << adid <<
", no. of FEB " <<
'\n';
642 for(
int chan=0; chan<32; chan++) {
645 float pe = (data->fAdc[chan]-chg_cal.second)/chg_cal.first;
650 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"\nfebC (mac5, channel, gain, pedestal, adc, pe) = (" << (int)data->fMac5 <<
", " << chan <<
", "
651 << chg_cal.first <<
", " << chg_cal.second <<
"," << data->fAdc[chan] <<
"," << pe <<
")\n";
654 pesmap[macs.back()].push_back(std::make_pair(chan,pe));
659 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"feb: " << (int)macs.back() <<
" ,chan : \t" << chan
660 <<
" ,pe: \t"<< pe <<
", adc:\t" << data->fAdc[chan] <<
", time: \t"<< data->fTs0
661 <<
" ,x: \t"<< postmp.X() <<
" ,y: \t" << postmp.Y() <<
" ,z: \t" << postmp.Z() <<
" petotal : "<< petot <<
'\n';
666 if(!(region==
"South" && layer==1)) {
669 hitpos.SetY(1.0*postmp.Y()+hitpos.Y());
671 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"!(region==South && layer==1) : \t" <<
" feb: " << (int)macs.back() <<
" ,chan : \t" << chan
672 <<
" ,pe: \t"<< pe <<
", adc:\t" << data->fAdc[chan] <<
", time: \t"<< data->fTs0
673 <<
" ,x: \t"<< postmp.X() <<
" ,y: \t" << postmp.Y() <<
" ,z: \t" << postmp.Z() <<
" petotal : "<< petot<<
'\n';
681 if(region!=
"South") {
683 hitpos.SetX(1.0*postmp.X()+hitpos.X());
696 hitpos.SetX(1.0*postmp.X()+hitpos.X());
698 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"vertical strips in South wall : \t" <<
" feb: " << (int)macs.back() <<
" ,chan : \t" << chan
699 <<
" ,pe: \t"<< pe <<
", adc:\t" << data->fAdc[chan] <<
", time: \t"<< data->fTs0
700 <<
" ,x: \t"<< postmp.X() <<
" ,y: \t" << postmp.Y() <<
" ,z: \t" << postmp.Z() <<
" petotal : "<< petot<<
'\n';
712 hitpos.SetZ(1.0*postmp.Z()+hitpos.Z());
716 mf::LogInfo(
"CRTHitRecoAlg: ") <<
" South wall z: \t" <<
" feb: " << (int)macs.back() <<
" ,chan : \t" << chan
717 <<
" ,pe: \t"<< pe <<
", adc:\t" << data->fAdc[chan] <<
", time: \t"<< data->fTs0
718 <<
" ,x: \t"<< postmp.X() <<
" ,y: \t" << postmp.Y() <<
" ,z: \t" << postmp.Z()<<
" petotal : "<< petot<<
'\n';
737 auto const& adsGeo = adGeo.SensitiveVolume(adsid_max);
746 ttrigs.push_back(data->fTs0);
747 t1trigs.push_back(data->fTs1);
748 tpos.push_back(postrig);
752 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"raw T0: " << data->fTs0 <<
" ,half lenth : \t" << adsGeo.HalfLength() <<
" ,delay: \t" <<
fPropDelay
753 <<
" ,corrected time: " << data->fTs0 - uint64_t(adsGeo.HalfLength()*
fPropDelay) <<
'\n';
756 if(region==
"South" && layer==1) {
757 southt0_h = data->fTs0;
759 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"southt0_h : " << layer <<
"\t" << southt0_h <<
'\n';
760 }
else if (region==
"South" && layer!=1){
761 southt0_v = data->fTs0;
763 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"southt0_v : " << layer <<
"\t"<< southt0_v <<
'\n';
774 std::sort(layID.begin(), layID.end());
775 layID.resize(
std::distance(layID.begin(), std::unique(layID.begin(), layID.end())));
778 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"size of layer ID : \t" <<layID.size() <<
'\n';
784 if(nabove==0 || layID.size() < 2) {
785 if(nabove==0 &&
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ") <<
"no channels above threshold!" <<
'\n';
786 if(layID.size()<2 &&
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ") <<
"no coincidence found" <<
'\n';
787 return FillCRTHit({},{},0,0,0,0,0,0,0,0,0,0,
"");
815 uint64_t t0_1=-5; uint64_t t0_2=-5;
816 uint64_t t1_1=-5; uint64_t t1_2=-5;
817 uint64_t t2_1=-5; uint64_t t2_2=-5;
818 int mac5_1=-5;
int mac5_2 = -5;
821 for(
auto const& infn: informationA) {
822 auto i = &infn - informationA.data();
823 auto const& adsGeo = adGeo.SensitiveVolume(infn.strip);
824 center = adsGeo.GetCenter();
827 mf::LogInfo(
"CRTHitRecoAlg: ")<<
"A type ----------> time: " << (
long long int)infn.t0 <<
" ,macs "<< (
int)infn.mac5s
828 <<
" ,chal "<< infn.channel
829 <<
" , position "<<infn.pos[2] <<
'\n';
831 if ((
int)infn.mac5s != (int)informationA[i+1].mac5s
and i < (
int)informationA.size()-1){
834 if ((
int)infn.mac5s % 2 == 0) t1_1 = infn.t0;
835 else t1_1 = informationA[i+1].t0;
836 if ((
int)informationA[i+1].mac5s % 2 != 0) t1_2 = informationA[i+1].
t0;
839 mf::LogInfo(
"CRTHitRecoAlg: ")<<
"t1: " << t1_1 <<
", t2:"<< t1_2 <<
", deltat : "<< int64_t(t1_1 - t1_2) <<
'\n';
841 float zaxixpos = 0.5*(int64_t(t1_1 - t1_2)/
fPropDelay);
843 posA = adsGeo.GetCenter() +
geo::Zaxis() * zaxixpos;
844 if (
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ")<<
"posA (==0): "<< posA<<
'\n';
848 mf::LogInfo(
"CRTHitRecoAlg: ")<< i <<
" ,1st mac5: "<< (int)infn.mac5s <<
" 1st time: " << (
long long int)infn.t0
849 <<
" ,2nd mac5: "<<(int)informationA[i+1].mac5s <<
", 2nd time " << (
long long int)informationA[i+1].t0 <<
" , deltaT: "
850 << int64_t(t1_1 - t1_2) <<
" , length: " << adsGeo.Length()
851 <<
" ,propagation delay: "<<
fPropDelay <<
" , pos z: " << zaxixpos
852 <<
" , center: " << adsGeo.GetCenter() <<
" , zaxis: "<<
geo::Zaxis() <<
" , half length: " << adsGeo.HalfLength()
853 <<
" , actual pos w.rt. z: " << posA <<
'\n';
857 for(
auto const& infn: informationB) {
861 mf::LogInfo(
"CRTHitRecoAlg: ")<<
" B type ----------> time: " << infn.t0 <<
" ,macs "<< (int)infn.mac5s
862 <<
" ,chal "<< infn.channel
863 <<
" , position "<<infn.pos[2] <<
'\n';
865 auto i = &infn - informationB.data();
866 auto const& adsGeo = adGeo.SensitiveVolume(infn.strip);
869 if ((
int)infn.mac5s != (
int)informationB[i+1].mac5s
and i < (
int)informationB.size()-1){
872 mac5_2 = (int)infn.mac5s;
873 t0_2 = (uint64_t)infn.t0;
875 if ((
int)infn.mac5s % 2 == 0) t2_1 = infn.t0;
876 else t2_1 = informationB[i+1].t0;
877 if ((
int)informationB[i+1].mac5s % 2 != 0) t2_2 = informationB[i+1].
t0;
881 mf::LogInfo(
"CRTHitRecoAlg: ")<<
"t1: " << t2_1 <<
", t2:"<< t2_2 <<
", deltat : "<< int64_t(t2_1 - t2_2) <<
'\n';
883 float zaxixpos = 0.5*(int64_t(t2_1 - t2_2)/
fPropDelay);
885 posB = adsGeo.GetCenter() +
geo::Zaxis() * zaxixpos;
886 if (
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ")<<
"posB (== 0): "<< posB<<
'\n';
890 mf::LogInfo(
"CRTHitRecoAlg: ")<< i <<
" ,1st mac5: "<< (int)infn.mac5s <<
" 1st time: " << (
long long int)infn.t0
891 <<
" ,2nd mac5: "<<(int)informationB[i+1].mac5s <<
", 2nd time " << (
long long int)informationB[i+1].t0 <<
" , deltaT: "
892 << int64_t(t2_1 - t2_2) <<
" , length: " << adsGeo.Length()
893 <<
" ,propagation delay: "<<
fPropDelay <<
" , pos z: " << zaxixpos
894 <<
" , center: " << adsGeo.GetCenter() <<
" , zaxis: "<<
geo::Zaxis() <<
" , half length: " << adsGeo.HalfLength()
895 <<
" , actual pos w.rt. z: " << posB <<
'\n';
901 int crossfeb =
std::abs(mac5_1 - mac5_2);
905 if (layer1 && layer2 && region!=
"South" && region!=
"North" ){
906 float avg = 0.5*(posA.Z() + posB.Z());
908 hitpos.SetX(hitpos.X()*1.0/nx);
909 hitpos.SetY(hitpos.Y()*1.0/nx);
913 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"z position in layer 1: "<< posA.Z() <<
" and in layer 2 "<< posB.Z()
914 <<
" average is "<< (posA.Z()+ posB.Z())/2. <<
" ,hitpos z " << hitpos[2] <<
'\n';
917 }
else if ((
int)informationA.size()==1
and (
int)informationB.size()==1
918 and (crossfeb == 7 or crossfeb == 5)
and
919 region!="South" && region!="North"){
921 int z_pos = int64_t(t0_1 - t0_2)/(uint64_t(2*
fPropDelay));
924 hitpos.SetZ(crossfebpos.Z());
925 hitpos.SetX(hitpos.X()*1.0/nx);
926 hitpos.SetY(hitpos.Y()*1.0/nx);
930 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"hello hi namaskar, hitpos z " << hitpos[2] <<
'\n';
932 }
else if (layer1 && region!=
"South" && region!=
"North"){
933 hitpos.SetZ(posA.Z());
934 hitpos.SetX(hitpos.X()*1.0/nx);
935 hitpos.SetY(hitpos.Y()*1.0/nx);
939 mf::LogInfo(
"CRTHitRecoAlg: ") <<
" same layer coincidence: z position in layer 1: "<< posA.Z() <<
" ,hitpos z " << hitpos[2] <<
'\n';
942 }
else if (layer2 && region!=
"South" && region!=
"North" ){
943 hitpos.SetZ(posB.Z());
944 hitpos.SetX(hitpos.X()*1.0/nx);
945 hitpos.SetY(hitpos.Y()*1.0/nx);
949 mf::LogInfo(
"CRTHitRecoAlg: ") <<
" same layer coincidence: z position in layer 2 "<< posB.Z() <<
" ,hitpos z " << hitpos[2] <<
'\n';
951 }
else if (region!=
"South" && region!=
"North" ){
956 if (
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ") <<
" In side CRTs [E/W] x: \t"<< hitpos[0] <<
" ,y: \t" << hitpos[1] <<
" ,z: \t" << hitpos[2]<<
'\n';
965 if(region==
"South") {
971 hitpos.SetX(hitpos.X()/nx);
972 hitpos.SetY(hitpos.Y()/ny);
973 hitpos.SetZ(hitpos.Z()/nz);
978 }
else if (region==
"North"){
991 hitpoint[0] = hitpos.X();
992 hitpoint[1] = hitpos.Y();
993 hitpoint[2] = hitpos.Z();
995 if (region==
"South" && hitpoint[0] >= 366. && hitpoint[1] > 200. &&
fVerbose)
996 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"I am looking for south wall : macs " << (int)macs.back() <<
" x: \t"<< hitpoint[0]
997 <<
" ,y: \t" << hitpoint[1] <<
" ,z: \t" << hitpoint[2] <<
'\n';
1000 if (region==
"North") mf::LogInfo(
"CRTHitRecoAlg: ") <<
"north wall x: \t"<< hitpoint[0] <<
" ,y: \t" << hitpoint[1] <<
" ,z: \t" << hitpoint[2]<<
'\n';
1002 if (
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ") <<
" nx: \t"<< nx <<
" ,ny: \t" << ny <<
" ,nz: \t" << nz<<
'\n';
1003 if (
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ") <<
" x: \t"<< hitpoint[0] <<
" ,y: \t" << hitpoint[1] <<
" ,z: \t" << hitpoint[2]<<
'\n';
1006 uint64_t thit = 0., t1hit = 0.;
1008 ttrigs.resize(
std::distance(ttrigs.begin(), std::unique(ttrigs.begin(), ttrigs.end())));
1012 t1trigs.resize(
std::distance(t1trigs.begin(), std::unique(t1trigs.begin(), t1trigs.end())));
1016 for(uint64_t
const t : ttrigs){
1022 mf::LogInfo(
"CRTHitRecoAlg: ") <<
" Inside the loop: t: \t"<< t <<
"\t" << uint64_t(200.*
fPropDelay)
1023 <<
" t0hit : "<< thit <<
" size ttrig: \t" << ttrigs.size()<<
'\n';
1026 thit=thit/uint64_t(ttrigs.size());
1028 for(
double const t1 : t1trigs)
1032 t1hit=t1hit/uint64_t(t1trigs.size());
1034 if (region==
"South" &&
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ") <<
"..................... Hello ....Welcome to Beam............"<<
'\n';
1035 if (
fVerbose) mf::LogInfo(
"CRTHitRecoAlg: ") <<
" <time>: T0: \t"<< thit <<
" T1 : "<< t1hit <<
" size ttrig: \t" << ttrigs.size()<<
'\n';
1038 mf::LogInfo(
"CRTHitRecoAlg: ") <<
"southt0_h: " << southt0_h <<
" ,southt0_v : " << southt0_v
1039 <<
" ,deltaT: \t" << int64_t(southt0_h - southt0_v) <<
'\n';
1040 if (
foutCSVFile)
filecsv << southt0_h <<
"\t" << southt0_v <<
"\t"<< int64_t(southt0_h - southt0_v) <<
"\n";
1043 auto const& adsGeo = adGeo.SensitiveVolume(adsid_max);
1044 if(region!=
"North" && region!=
"South"){
1045 hitpointerr[0] = (xmax-
xmin)/sqrt(12);
1046 hitpointerr[1] = (ymax-ymin)/sqrt(12);
1047 hitpointerr[2] = (zmax-zmin)/sqrt(12);
1051 if(region==
"North"){
1052 hitpointerr[0] = (xmax-
xmin)/sqrt(12);
1053 hitpointerr[1] = (ymax-ymin)/sqrt(12);
1054 hitpointerr[2] = (zmax-zmin)/sqrt(12);
1057 if(region==
"South"){
1058 hitpointerr[0] = adsGeo.HalfWidth1()*2/sqrt(12);
1059 hitpointerr[1] = adsGeo.HalfWidth1()*2/sqrt(12);
1060 hitpointerr[2] = (zmax-zmin)/sqrt(12);
1063 Long64_t thit1=(Long64_t)(thit-GlobalTrigger[(
int)macs.at(0)]);
1067 hitpoint[1],hitpointerr[1],hitpoint[2],hitpointerr[2],region);
int GetMINOSLayerID(size_t adid)
then if[["$THISISATEST"==1]]
CRTHit FillCRTHit(vector< uint8_t > tfeb_id, map< uint8_t, vector< pair< int, float >>> tpesmap, float peshit, uint64_t time0, Long64_t time1, int plane, double x, double ex, double y, double ey, double z, double ez, string tagger)
bool foutCSVFile
FCL input: Write a CSV File?
std::size_t size(FixedBins< T, C > const &) noexcept
double fPEThresh
threshold[PE] above which charge amplitudes used in hit reco
size_t MacToAuxDetID(uint8_t mac, int chan)
uint64_t fSiPMtoFEBdelay
SiPM to FEB cable induced delay: 11.6 [ns].
int AuxDetRegionNameToNum(string reg)
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
CRTCommonUtils * fCrtutils
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
virtual std::pair< double, double > getSideCRTCalibrationMap(int mac5, int chan) const =0
for($it=0;$it< $RaceTrack_number;$it++)
geo::GeometryCore const * fGeometryService
const icarusDB::IICARUSChannelMap * fChannelMap
pair< uint8_t, uint8_t > ADToMac(size_t adid)
return match has_match and(match.match_pdg==11 or match.match_pdg==-11)
double fPropDelay
propegation time [ns/cm]
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
constexpr Vector Zaxis()
Returns a z axis vector of the specified type.
string GetAuxDetRegion(size_t adid)
TVector3 ChanToWorldCoords(const uint8_t mac, const int chan)
int ChannelToAuxDetSensitiveID(uint8_t mac, int chan)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
BEGIN_PROLOG could also be cout