Loading...
Searching...
No Matches
MdcCalib Class Referenceabstract
#include <MdcCalib.h>
Inheritance diagram for MdcCalib:
Public Member Functions | |
| MdcCalib () | |
| virtual | ~MdcCalib () |
| virtual void | initialize (TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)=0 |
| virtual void | setParam (MdcCalParams ¶m)=0 |
| virtual int | fillHist (MdcCalEvent *event)=0 |
| virtual int | updateConst (MdcCalibConst *calconst)=0 |
| virtual void | printCut () const =0 |
| virtual void | clear ()=0 |
Detailed Description
Definition at line 35 of file MdcCalib.h.
Constructor & Destructor Documentation
◆ MdcCalib()
| MdcCalib::MdcCalib | ( | ) |
Definition at line 36 of file MdcCalib.cxx.
36 {
37 m_nEvt = 0;
38 m_cut1 = 0;
39 m_nTrkAfTrkCut = 0;
40 m_cut2 = 0;
41 m_cut3 = 0;
42 m_cut4 = 0;
43 m_cut5 = 0;
44 m_cut6 = 0;
45 m_cut7 = 0;
46 m_nTrkCal = 0;
47 m_nGrPoint = 0;
48 fgReadWireEff = false;
49
50 int lay;
51 for ( lay = 0; lay < 43; lay++ )
52 {
53 if ( lay < 15 ) m_nEntr[lay] = 1;
54 else m_nEntr[lay] = 2;
55 // m_nEntr[lay] = 3;
56 }
57 m_dwid = 0.5; // mm
58 m_fgIni = false;
59
60 m_phiWid = PI2 / (double)NPhiBin;
61 m_theWid = 2.0 / (double)NThetaBin;
62
63 m_nEvtNtuple = 0;
64
66 {
67 if ( lay < 8 ) m_nBin[lay] = 12;
68 else m_nBin[lay] = 16;
69 }
70
71 // setting boundary layer flags
73 {
74 if ( ( 0 == lay ) || ( 7 == lay ) || ( 8 == lay ) || ( 19 == lay ) || ( 20 == lay ) ||
75 ( 35 == lay ) || ( 36 == lay ) || ( 42 == lay ) )
76 m_layBound[lay] = true;
77 else m_layBound[lay] = false;
78 }
79}
◆ ~MdcCalib()
|
virtual |
Definition at line 81 of file MdcCalib.cxx.
81{}
Member Function Documentation
◆ clear()
|
pure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 83 of file MdcCalib.cxx.
83 {
84 int lay;
85 for ( lay = 0; lay < m_nlayer; lay++ )
86 {
87 delete m_htraw[lay];
88 delete m_htdr[lay];
89 delete m_hresInc[lay];
90 delete m_hresExc[lay];
91 delete m_hresAve[lay];
92 delete m_hadc[lay];
93 for ( int lr = 0; lr < 2; lr++ )
94 {
95 delete m_htdrlr[lay][lr];
96 delete m_hreslrInc[lay][lr];
97 delete m_hreslrExc[lay][lr];
98 delete m_hreslrAve[lay][lr];
99 }
100 }
101
102 delete m_effNtrk;
103 delete m_effNtrkRecHit;
104 delete m_hresAllInc;
105 delete m_hresAllExc;
106 delete m_hresAllAve;
107 for ( int i = 0; i < 14; i++ )
108 {
109 delete m_hresAveAllQ[i];
110 delete m_hresAveOutQ[i];
111 }
112 for ( lay = 0; lay < 43; lay++ )
113 {
114 for ( int i = 0; i < 14; i++ ) delete m_hresAveLayQ[lay][i];
115 }
116 delete m_hresInnInc;
117 delete m_hresInnExc;
118 delete m_hresStpInc;
119 delete m_hresStpExc;
120 delete m_hresOutInc;
121 delete m_hresOutExc;
122 delete m_hresAllExcUp;
123 delete m_hresAllExcDw;
124 delete m_hresInnExcUp;
125 delete m_hresInnExcDw;
126 delete m_hresStpExcUp;
127 delete m_hresStpExcDw;
128 delete m_hresOutExcUp;
129 delete m_hresOutExcDw;
130 for ( int iEs = 0; iEs < m_param.nEsFlag; iEs++ ) delete m_hTes[iEs];
131 delete m_hbbTrkFlg;
132 delete m_hTesAll;
133 delete m_hTesGood;
134 delete m_hTesAllFlag;
135 delete m_hTesRec;
136 delete m_hTesCalFlag;
137 delete m_hTesCalUse;
138 delete m_hnRawHit;
139 delete m_hpt;
140 delete m_hpMax;
141 delete m_hpMaxCms;
142 delete m_hptPos;
143 delete m_hptNeg;
144 delete m_hp;
145 delete m_hp_cms;
146 delete m_hpPos;
147 delete m_hpNeg;
148 delete m_hpPoscms;
149 delete m_hpNegcms;
150 delete m_hp_cut;
151 delete m_hchisq;
152 delete m_hnTrk;
153 delete m_hnTrkCal;
154 delete m_hnhitsRec;
155 delete m_hnhitsRecInn;
156 delete m_hnhitsRecStp;
157 delete m_hnhitsRecOut;
158 delete m_hnhitsCal;
159 delete m_hnhitsCalInn;
160 delete m_hnhitsCalStp;
161 delete m_hnhitsCalOut;
162 delete m_wirehitmap;
163 delete m_layerhitmap;
164 delete m_hnoisephi;
165 delete m_hnoiselay;
166 delete m_hnoisenhits;
167 delete m_hratio;
168 delete m_hdr;
169 delete m_hphi0;
170 delete m_hkap;
171 delete m_hdz;
172 delete m_htanl;
173 delete m_hcosthe;
174 delete m_hcostheNeg;
175 delete m_hcosthePos;
176 delete m_hx0;
177 delete m_hy0;
178 delete m_hdelZ0;
179 delete m_grX0Y0;
180 delete m_hitEffAll;
181 delete m_hitEffRaw;
182 delete m_hitEffRec;
184 int thbin;
186 {
191 for ( thbin = 0; thbin < NThetaBin; thbin++ )
192 {
197 }
198 }
199 for ( thbin = 0; thbin < NThetaBin; thbin++ )
200 {
201 delete m_ppThe[thbin];
202 delete m_pnThe[thbin];
203 delete m_ppTheCms[thbin];
204 delete m_pnTheCms[thbin];
205 }
206
207 for ( unsigned i = 0; i < m_hr2dInc.size(); i++ )
208 {
209 delete m_hr2dInc[i];
210 delete m_hr2dExc[i];
211 }
212 m_hr2dInc.clear();
213 m_hr2dExc.clear();
214 m_mapr2d.clear();
215
216 delete m_fdTime;
217 delete m_fdAdc;
218 delete m_fdres;
219 delete m_fdresAve;
220 delete m_fdres2d;
221 delete m_fdcom;
222 delete m_fdResQ;
223}
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition FoamA.h:85
Referenced by GrXtMdcCalib::clear(), PreT0MdcCalib::clear(), QtMdcCalib::clear(), T0MdcCalib::clear(), Wr2dMdcCalib::clear(), WrMdcCalib::clear(), XtInteMdcCalib::clear(), and XtMdcCalib::clear().
◆ fillHist()
|
pure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 767 of file MdcCalib.cxx.
767 {
768 IMessageSvc* msgSvc;
771 log << MSG::DEBUG << "MdcCalib::fillHist()" << endmsg;
772
773 int i;
774 int k;
775 int lay;
776 int cel;
777 int wir;
778 int lr;
779 int stat;
780
781 int hid;
783 int iEntr;
785
786 int phiBin;
787 int phiBinCms;
788 int theBin;
789 int theBinCms;
790 double lamda;
791 double theta;
792
793 double qhit;
794 double traw;
795 double tdr;
796 double doca;
797 double resiInc;
798 double resiExc;
799 double entr;
800 double pt;
801 double p;
803 double chisq;
804 double ecm = m_param.ecm;
805 double xboost = m_param.boostPar[0] * ecm;
806 double yboost = m_param.boostPar[1];
807 double zboost = m_param.boostPar[2];
808 HepLorentzVector p4;
809
810 double dr;
811 double phi0;
812 double dz;
813 double kap;
814 double tanl;
815
816 double x0;
817 double y0;
818 double zminus = 9999.0;
819 double zplus = -9999.0;
820
821 double hitphi;
822 double philab;
823 double phicms;
824 double thetacms;
825 double costheCMS;
826
827 double dphi;
828 double wphi;
829 double xx;
830 double yy;
831 double rr;
832
833 int nhitlay;
835 bool fgTrk;
836
837 int ntrk = event->getNTrk();
838 int nhitRec;
839 int nhitRecInn;
840 int nhitRecStp;
841 int nhitRecOut;
842 int nhitCal;
843 int nhitCalInn;
844 int nhitCalStp;
845 int nhitCalOut;
846 MdcCalRecTrk* rectrk;
847 MdcCalRecHit* rechit;
848
849 int fgAdd[43]; // for calculating layer efficiency
850
851 // read dead wire
852 if ( !fgReadWireEff )
853 {
855 {
856 int ncel = m_mdcGeomSvc->Layer( lay )->NCell();
857 for ( cel = 0; cel < ncel; cel++ )
858 {
859 double eff = m_mdcFunSvc->getWireEff( lay, cel );
860 if ( eff > 0.5 ) m_fgGoodWire[lay][cel] = true;
861 else m_fgGoodWire[lay][cel] = false;
862 int wireid = m_mdcGeomSvc->Wire( lay, cel )->Id();
863 if ( eff < 0.9 )
864 cout << "dead channel: " << setw( 5 ) << lay << setw( 5 ) << cel << setw( 8 )
865 << wireid << endl;
866 }
867 }
868 fgReadWireEff = true;
869
872 {
873 m_xe[wir] = m_mdcGeomSvc->Wire( wir )->Backward().x();
874 m_ye[wir] = m_mdcGeomSvc->Wire( wir )->Backward().y();
875 m_ze[wir] = m_mdcGeomSvc->Wire( wir )->Backward().z();
876 m_xw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().x();
877 m_yw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().y();
878 m_zw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().z();
879 fwpc << setw( 6 ) << wir << setw( 15 ) << m_xe[wir] << setw( 15 ) << m_ye[wir]
880 << setw( 15 ) << m_ze[wir] << setw( 15 ) << m_xw[wir] << setw( 15 ) << m_yw[wir]
881 << setw( 15 ) << m_zw[wir] << endl;
882 }
883 fwpc.close();
884 }
885
886 int nRawHit = event->getNRawHitTQ();
887 m_hnRawHit->Fill( nRawHit );
888
889 IDataProviderSvc* eventSvc = NULL;
890 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
891
892 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
893 if ( !eventHeader )
894 {
895 log << MSG::FATAL << "Could not find Event Header" << endmsg;
896 return ( StatusCode::FAILURE );
897 }
898 int iRun = eventHeader->runNumber();
899 int iEvt = eventHeader->eventNumber();
900
901 int esTimeflag = event->getEsFlag();
902 double tes = event->getTes();
903 bool esCutFg = event->getEsCutFlag();
904 int iEs = event->getNesCutFlag();
905 // calculate the efficiency of Bhabha event
906 if ( -1 != esTimeflag )
907 {
908 SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc, "/Event/Recon/RecMdcTrackCol" );
909 if ( !newtrkCol )
910 {
911 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endmsg;
912 return ( StatusCode::FAILURE );
913 }
914 int nGoodTrk = 0;
915 int icharge = 0;
916 Vp4 p4p;
917 Vp4 p4m;
918 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
919 for ( ; it_trk != newtrkCol->end(); it_trk++ )
920 {
922 HepLorentzVector p4 = ( *it_trk )->p4( mass );
923 icharge = ( *it_trk )->charge();
924 if ( icharge > 0 ) p4p.push_back( p4 );
925 else p4m.push_back( p4 );
926 }
927 if ( 1 == p4p.size() * p4m.size() )
928 {
929 double dang = p4p[0].vect().angle( p4m[0].vect() );
930 m_hbbTrkFlg->Fill( 1 );
931 if ( dang > 2.94 ) { m_hbbTrkFlg->Fill( 2 ); }
932 }
933 }
934 m_hTesAll->Fill( tes );
935 // cout << "tes " << tes << endl;
936 // if (-1 != esTimeflag) m_hTesGood->Fill(tes);
937 if ( -1 == esTimeflag ) return -1;
938 m_hTesGood->Fill( tes );
939 m_hTesAllFlag->Fill( esTimeflag );
940 if ( ntrk > 0 ) m_hTesRec->Fill( tes );
941 if ( ( iEs >= 0 ) && ( iEs < m_param.nEsFlag ) ) m_hTes[iEs]->Fill( tes );
942
943 if ( esCutFg ) m_hTesCalFlag->Fill( tes );
944 else return -1;
945
946 if ( !m_fgIni )
947 {
949 {
950 if ( lay < 8 ) m_docaMax[lay] = m_param.maxDocaInner;
951 else m_docaMax[lay] = m_param.maxDocaOuter;
952 }
953 m_fgIni = true;
954 }
955
956 bool trkFlag[200]; // for calculating hit efficiency without impact of track fitting
957 for ( i = 0; i < 200; i++ ) trkFlag[i] = false;
958
959 int ntrkCal = 0;
960 double pTrk[2];
961 double pTrkcms[2];
962 double hitphiplus = 9999.0;
963 double hitthetaplus = 9999.0;
964 double hitphiminus = -9999.0;
965 double hitthetaminus = -9999.0;
966 Vp4 pp;
967 Vp4 pm;
968 pp.clear();
969 pm.clear();
970 int phibinp;
971 int phibinm;
972
973 m_hnTrk->Fill( ntrk );
974 if ( ( ntrk < m_param.nTrkCut[0] ) || ( ntrk > m_param.nTrkCut[1] ) )
975 {
976 m_cut1++;
977 return -1;
978 }
979 // if(ntrk > 2) return -1;
980 for ( i = 0; i < ntrk; i++ )
981 {
982 m_nTrkAfTrkCut++;
983 fgTrk = true;
984 rectrk = event->getRecTrk( i );
985 nhitRec = rectrk->getNHits();
986 m_hnhitsRec->Fill( nhitRec );
987
989
990 nhitRecInn = 0;
991 nhitRecStp = 0;
992 nhitRecOut = 0;
993 for ( k = 0; k < nhitRec; k++ )
994 {
995 rechit = rectrk->getRecHit( k );
996 lay = rechit->getLayid();
997 doca = rechit->getDocaExc();
998 resiExc = rechit->getResiExc();
999 fgHitLay[lay] = true;
1000
1001 if ( lay < 8 ) nhitRecInn++;
1002 else if ( lay < 20 ) nhitRecStp++;
1003 else nhitRecOut++;
1004 }
1005 m_hnhitsRecInn->Fill( nhitRecInn );
1006 m_hnhitsRecStp->Fill( nhitRecStp );
1007 m_hnhitsRecOut->Fill( nhitRecOut );
1008
1009 // get momentum
1010 pt = rectrk->getPt();
1011 p = rectrk->getP();
1012
1013 // boost P to the cms
1014 p4 = rectrk->getP4();
1015 HepLorentzVector psip( xboost, yboost, zboost, ecm );
1016 Hep3Vector boostv = psip.boostVector();
1017 p4.boost( -boostv );
1018 p_cms = p4.rho();
1019 phicms = p4.phi();
1021 thetacms = p4.theta();
1022 costheCMS = cos( thetacms );
1024 p4.boost( boostv );
1025 // cout << setw(15) << p << setw(15) << p_cms << setw(15) << costheCMS << endl;
1026
1027 // cos(theta) cut
1028 if ( ( costheCMS < m_param.costheCut[0] ) || ( costheCMS > m_param.costheCut[1] ) )
1029 {
1030 m_cut2++;
1031 continue;
1032 }
1033
1034 // dr cut
1035 dr = rectrk->getDr();
1036 if ( fabs( dr ) > m_param.drCut )
1037 {
1038 m_cut3++;
1039 continue;
1040 }
1041
1042 // dz cut
1043 dz = rectrk->getDz();
1044 if ( fabs( dz ) > m_param.dzCut )
1045 {
1046 m_cut4++;
1047 continue;
1048 }
1049
1050 // momentum cut
1051 if ( ( fabs( p ) < m_param.pCut[0] ) || ( fabs( p ) > m_param.pCut[1] ) )
1052 {
1053 m_cut5++;
1054 continue;
1055 }
1056
1057 // charge cut, to check +/- asymmetry
1058 int charge = 1;
1059 if ( p < 0 ) charge = -1;
1061
1062 // if(! fgTrk) continue;
1063
1064 // hit layer cut
1065 nhitlay = 0;
1067 {
1068 if ( fgHitLay[lay] ) nhitlay++;
1069 }
1070 if ( nhitlay < m_param.nHitLayCut )
1071 {
1072 m_cut6++;
1073 continue;
1074 }
1075
1076 // nhit cut
1077 if ( nhitRec < m_param.nHitCut )
1078 {
1079 m_cut7++;
1080 continue;
1081 }
1082
1083 // bool fgNoise = rectrk->getFgNoiseRatio();
1084 // if(m_param.noiseCut && (!fgNoise)) continue;
1085 // cout << setw(10) << iRun << setw(15) << iEvt << setw(5) << fgNoise << endl;
1086
1087 // if(! ((fgHitLay[0]||fgHitLay[1]) && (fgHitLay[41]||fgHitLay[42])) ){
1088 // continue;
1089 // }
1090
1091 // calculate cell on the track
1092 int cellTrkPass[43];
1093 bool fgPass = getCellTrkPass( event, i, cellTrkPass );
1094 for ( lay = 0; lay < m_nlayer; lay++ )
1095 {
1096 fgAdd[lay] = 0;
1097 // if((16==lay) || (18==lay) || (19==lay) || (41==lay)){ // hv2200 2009-3
1098 if ( ( 15 == lay ) || ( 16 == lay ) || ( 18 == lay ) || ( 19 == lay ) || ( 40 == lay ) ||
1099 ( 41 == lay ) || ( 42 == lay ) )
1100 {
1101 int iCell = cellTrkPass[lay];
1102 if ( fgPass && ( iCell >= 0 ) && m_fgGoodWire[lay][iCell] ) m_effNtrk->Fill( lay );
1103 else fgAdd[lay] = 1;
1104 }
1105 else { m_effNtrk->Fill( lay ); }
1106 }
1107
1108 m_nTrkCal++;
1109 chisq = rectrk->getChisq();
1110 m_hchisq->Fill( chisq );
1111
1112 if ( pt > 0 )
1113 {
1114 m_hpt->Fill( pt );
1115 m_hptPos->Fill( pt );
1116 m_hp->Fill( p );
1117 m_hp_cms->Fill( p_cms );
1118 m_hpPos->Fill( p );
1119 m_hpPoscms->Fill( p_cms );
1120 }
1121 else
1122 {
1123 m_hpt->Fill( -1.0 * pt );
1124 m_hptNeg->Fill( -1.0 * pt );
1125 m_hp->Fill( -1.0 * p );
1126 m_hp_cms->Fill( -1.0 * p_cms );
1127 m_hpNeg->Fill( -1.0 * p );
1128 m_hpNegcms->Fill( -1.0 * p_cms );
1129 }
1130 if ( 2 == ntrk )
1131 {
1132 pTrk[i] = fabs( p );
1133 pTrkcms[i] = fabs( p_cms );
1134 }
1135
1136 dr = rectrk->getDr();
1137 phi0 = rectrk->getPhi0();
1138 kap = rectrk->getKappa();
1139 dz = rectrk->getDz();
1140 tanl = rectrk->getTanLamda();
1141 lamda = atan( tanl );
1142 theta = HFPI - lamda;
1143
1144 m_hdr->Fill( dr );
1145 m_hphi0->Fill( phi0 );
1146 m_hkap->Fill( kap );
1147 m_hdz->Fill( dz );
1148 m_htanl->Fill( tanl );
1149 m_hcosthe->Fill( cos( theta ) );
1152
1153 philab = phi0 + HFPI;
1155 // cout << setw(15) << phi0 << setw(15) << philab << setw(15) << phicms << endl;
1156
1157 phiBin = (int)( philab / m_phiWid );
1158 phiBinCms = (int)( phicms / m_phiWid );
1159 theBin = (int)( ( cos( theta ) + 1.0 ) / m_theWid );
1160 theBinCms = (int)( ( cos( thetacms ) + 1.0 ) / m_theWid );
1161 if ( phiBin < 0 ) phiBin = 0;
1162 if ( phiBin >= NPhiBin ) phiBin = NPhiBin - 1;
1163 if ( phiBinCms < 0 ) phiBinCms = 0;
1164 if ( phiBinCms >= NPhiBin ) phiBinCms = NPhiBin - 1;
1165 if ( theBin < 0 ) theBin = 0;
1166 if ( theBin >= NThetaBin ) theBin = NThetaBin - 1;
1167 if ( theBinCms < 0 ) theBinCms = 0;
1168 if ( theBinCms >= NThetaBin ) theBinCms = NThetaBin - 1;
1169
1170 if ( pt > 0 )
1171 {
1172 m_ppPhi[phiBin]->Fill( p );
1173 m_ppPhiCms[phiBinCms]->Fill( p_cms );
1174 m_ppThe[theBin]->Fill( p );
1175 m_ppTheCms[theBinCms]->Fill( p_cms );
1176 m_ppThePhi[theBin][phiBin]->Fill( p );
1177 m_ppThePhiCms[theBinCms][phiBinCms]->Fill( p_cms );
1178 }
1179 else
1180 {
1181 m_pnPhi[phiBin]->Fill( -1.0 * p );
1182 m_pnPhiCms[phiBinCms]->Fill( -1.0 * p_cms );
1183 m_pnThe[theBin]->Fill( -1.0 * p );
1184 m_pnTheCms[theBinCms]->Fill( -1.0 * p_cms );
1185 m_pnThePhi[theBin][phiBin]->Fill( -1.0 * p );
1186 m_pnThePhiCms[theBinCms][phiBinCms]->Fill( -1.0 * p_cms );
1187 }
1188
1189 x0 = dr * cos( phi0 );
1190 y0 = dr * sin( phi0 );
1191 m_hx0->Fill( x0 );
1192 m_hy0->Fill( y0 );
1193 if ( m_nGrPoint < 10000 )
1194 {
1195 m_grX0Y0->SetPoint( m_nGrPoint, x0, y0 );
1196 m_nGrPoint++;
1197 }
1198
1199 if ( kap < 0 )
1200 {
1201 zminus = dz;
1202 pm.push_back( p4 );
1203 phibinm = phiBinCms;
1204 }
1205 else
1206 {
1207 zplus = dz;
1208 pp.push_back( p4 );
1209 phibinp = phiBinCms;
1210 }
1211
1212 // cout << "phi = " << setw(15) << philab << setw(15) << philab*180./3.14159 << setw(15)
1213 // << p << endl;
1214 ntrkCal++;
1215 trkFlag[i] = true;
1216 nhitCal = 0;
1217 nhitCalInn = 0;
1218 nhitCalStp = 0;
1219 nhitCalOut = 0;
1220 for ( k = 0; k < nhitRec; k++ )
1221 {
1222 rechit = rectrk->getRecHit( k );
1223
1224 lay = rechit->getLayid();
1225 cel = rechit->getCellid();
1226 lr = rechit->getLR();
1227 stat = rechit->getStat();
1228 doca = rechit->getDocaExc();
1229 resiInc = rechit->getResiIncLR();
1230 resiExc = rechit->getResiExcLR();
1231 entr = rechit->getEntra();
1232 tdr = rechit->getTdrift();
1234 wir = m_mdcGeomSvc->Wire( lay, cel )->Id();
1235
1236 m_cel[lay] = (long)cel;
1237 m_lr[lay] = (long)lr;
1238 m_run[lay] = iRun;
1239 m_evt[lay] = iEvt;
1240 m_doca[lay] = doca;
1241 m_dm[lay] = rechit->getDmeas();
1242 m_tdr[lay] = tdr;
1243 m_tdc[lay] = traw;
1244 m_entr[lay] = entr * 180.0 / 3.14;
1245 m_zhit[lay] = rechit->getZhit();
1246 m_qhit[lay] = rechit->getQhit();
1247 m_p[lay] = p_cms;
1248 m_pt[lay] = pt;
1249 m_phi0[lay] = phi0;
1250 m_tanl[lay] = tanl;
1251 qhit = rechit->getQhit();
1252 // cout << "layer " << setw(5) << lay << setw(5) << cel << setw(5) << lr <<
1253 // endl; calculating hitphi
1254 xx =
1255 ( m_zhit[lay] - m_zw[wir] ) * ( m_xe[wir] - m_xw[wir] ) / ( m_ze[wir] - m_zw[wir] ) +
1256 m_xw[wir];
1257 yy =
1258 ( m_zhit[lay] - m_zw[wir] ) * ( m_ye[wir] - m_yw[wir] ) / ( m_ze[wir] - m_zw[wir] ) +
1259 m_yw[wir];
1260 rr = sqrt( ( xx * xx ) + ( yy * yy ) );
1261 dphi = fabs( doca ) / m_radii[lay];
1262
1263 if ( yy >= 0 ) wphi = acos( xx / rr );
1265 if ( 1 == lr ) hitphi = wphi + dphi; // mention
1266 else hitphi = wphi - dphi;
1269
1270 m_hitphi[lay] = hitphi;
1271
1272 if ( ( fabs( doca ) > m_docaMax[lay] ) || ( fabs( resiExc ) > m_param.resiCut[lay] ) )
1273 { continue; }
1274
1275 if ( m_param.fgAdjacLayerCut )
1276 {
1277 if ( 0 == lay )
1278 {
1279 if ( !fgHitLay[1] ) continue;
1280 }
1281 else if ( 42 == lay )
1282 {
1283 if ( !fgHitLay[41] ) continue;
1284 }
1285 else
1286 {
1287 if ( ( !fgHitLay[lay - 1] ) && ( !fgHitLay[lay + 1] ) ) continue;
1288
1289 // for boundary layers
1290 if ( m_param.fgBoundLayerCut && m_layBound[lay] &&
1291 ( ( !fgHitLay[lay - 1] ) || ( !fgHitLay[lay + 1] ) ) )
1292 continue;
1293 }
1294 }
1295
1296 if ( ( 1 == m_param.hitStatCut ) && ( 1 != stat ) ) continue;
1297
1298 // fill xtplot tree
1299 if ( ( 1 == m_param.fillNtuple ) && ( m_nEvtNtuple < m_param.nEvtNtuple ) )
1300 { m_xtTuple[lay]->write(); }
1301
1302 if ( 1 == m_param.hitStatCut )
1303 {
1304 if ( ( 0 == fgAdd[lay] ) && ( 1 == stat ) )
1305 {
1306 m_effNtrkRecHit->Fill( lay );
1307 fgAdd[lay] = 1;
1308 }
1309 }
1310 else
1311 {
1312 if ( 0 == fgAdd[lay] )
1313 {
1314 m_effNtrkRecHit->Fill( lay );
1315 fgAdd[lay] = 1;
1316 }
1317 }
1318
1319 nhitCal++;
1320 if ( lay < 8 ) nhitCalInn++;
1321 else if ( lay < 20 ) nhitCalStp++;
1322 else nhitCalOut++;
1323
1324 m_wirehitmap->Fill( wir );
1325 m_layerhitmap->Fill( lay );
1326
1327 m_htraw[lay]->Fill( traw );
1328 m_htdr[lay]->Fill( tdr );
1329 m_htdrlr[lay][lr]->Fill( tdr );
1330 m_hadc[lay]->Fill( qhit );
1331
1332 m_hresAllInc->Fill( resiInc );
1333 m_hresAllExc->Fill( resiExc );
1334 double resiAve = 0.5 * ( resiInc + resiExc );
1335 m_hresAllAve->Fill( resiAve );
1336 if ( yy > 0 ) m_hresAllExcUp->Fill( resiExc );
1337 else m_hresAllExcDw->Fill( resiExc );
1338
1339 if ( lay < 8 )
1340 {
1341 m_hresInnInc->Fill( resiInc );
1342 m_hresInnExc->Fill( resiExc );
1343 if ( yy > 0 ) m_hresInnExcUp->Fill( resiExc );
1344 else m_hresInnExcDw->Fill( resiExc );
1345 }
1346 else if ( lay < 20 )
1347 {
1348 m_hresStpInc->Fill( resiInc );
1349 m_hresStpExc->Fill( resiExc );
1350 if ( yy > 0 ) m_hresStpExcUp->Fill( resiExc );
1351 else m_hresStpExcDw->Fill( resiExc );
1352 }
1353 else
1354 {
1355 m_hresOutInc->Fill( resiInc );
1356 m_hresOutExc->Fill( resiExc );
1357 if ( yy > 0 ) m_hresOutExcUp->Fill( resiExc );
1358 else m_hresOutExcDw->Fill( resiExc );
1359 }
1360
1361 int qbin = (int)( ( qhit - 100.0 ) / 100.0 );
1362 if ( qbin >= 0 && qbin < 14 )
1363 {
1364 m_hresAveAllQ[qbin]->Fill( resiAve );
1365 m_hresAveLayQ[lay][qbin]->Fill( resiAve );
1366 if ( lay > 7 ) m_hresAveOutQ[qbin]->Fill( resiAve );
1367 }
1368
1369 if ( ( !m_param.cosmicDwTrk ) || ( m_param.cosmicDwTrk && ( yy < 0 ) ) )
1370 {
1371 m_hresInc[lay]->Fill( resiInc );
1372 m_hreslrInc[lay][lr]->Fill( resiInc );
1373 m_hresExc[lay]->Fill( resiExc );
1374 m_hreslrExc[lay][lr]->Fill( resiExc );
1375 m_hresAve[lay]->Fill( resiAve );
1376 m_hreslrAve[lay][lr]->Fill( resiAve );
1377 }
1378
1380 if ( iPhi >= 20 ) iPhi = 19;
1381 m_hresphi[lay][iPhi]->Fill( ( resiInc + resiExc ) * 0.5 );
1382
1383 // bin = (int)(fabs(doca) / m_dwid);
1385 iEntr = m_mdcFunSvc->getSdEntrIndex( entr );
1386 if ( 1 == m_nEntr[lay] ) { iEntr = 0; }
1387 else if ( 2 == m_nEntr[lay] )
1388 {
1389 if ( entr > 0.0 ) iEntr = 1;
1390 else iEntr = 0;
1391 }
1392 else if ( 3 == m_nEntr[lay] )
1393 {
1394 double entrBinWid = 0.174532925; // 10 degree
1395 iEntr = (int)( fabs( entr ) / entrBinWid );
1396 if ( iEntr > 2 ) iEntr = 2;
1397 }
1399 {
1402 {
1403 hid = m_mapr2d[key];
1404 m_hr2dInc[hid]->Fill( resiInc );
1405 m_hr2dExc[hid]->Fill( resiExc );
1406 }
1407 }
1408
1409 if ( ( tdr > 0 ) && ( tdr < 750 ) )
1410 {
1413 m_hr2t[lay][iEntr][lr][bin]->Fill( resiExc );
1414 }
1415 } // loop of nhits
1416 m_nEvtNtuple++;
1417 m_hnhitsCal->Fill( nhitCal );
1418 m_hnhitsCalInn->Fill( nhitCalInn );
1419 m_hnhitsCalStp->Fill( nhitCalStp );
1420 m_hnhitsCalOut->Fill( nhitCalOut );
1421 } // end of track loop
1422 m_hnTrkCal->Fill( ntrkCal );
1423 if ( 2 == ntrkCal )
1424 {
1425 if ( pTrk[0] > pTrk[1] ) m_hpMax->Fill( pTrk[0] );
1426 else m_hpMax->Fill( pTrk[1] );
1427
1428 if ( pTrkcms[0] > pTrkcms[1] ) m_hpMaxCms->Fill( pTrkcms[0] );
1429 else m_hpMaxCms->Fill( pTrkcms[1] );
1430 }
1431 if ( ntrkCal > 0 ) m_hTesCalUse->Fill( tes );
1432
1433 double delZ0;
1434 if ( ( fabs( zminus ) < 9000.0 ) && ( fabs( zplus ) < 9000.0 ) ) delZ0 = zplus - zminus;
1435 m_hdelZ0->Fill( delZ0 );
1436
1437 if ( 1 == pp.size() * pm.size() )
1438 {
1439 HepLorentzVector ptot = pp[0] + pm[0];
1440 bool fourmomcut = false;
1441 // fourmomcut = (ptot.x()>0.02 && ptot.x()<0.06) && (fabs(ptot.y()) < 0.02)
1442 // && (ptot.z()>-0.01 && ptot.z()<0.03) && (ptot.e()>3.4 && ptot.e()<4.0);
1443 fourmomcut = ( fabs( ptot.x() - 0.04 ) < 0.026 ) && ( fabs( ptot.y() ) < 0.026 ) &&
1444 ( fabs( ptot.z() - 0.005 ) < 0.036 ) && ( fabs( ptot.e() - ecm ) < 0.058 );
1445 // cout << "x = " << ptot.x() << ", y = " << ptot.y() << ", z = " << ptot.z() << ", e = "
1446 // << ptot.e() << endl;
1447 if ( fourmomcut )
1448 {
1449 HepLorentzVector psip( xboost, yboost, zboost, ecm );
1450 Hep3Vector boostv = psip.boostVector();
1451 pp[0].boost( -boostv );
1452 pm[0].boost( -boostv );
1453 m_hp_cut->Fill( pp[0].rho() );
1454 m_hp_cut->Fill( pm[0].rho() );
1455 }
1456 }
1457
1458 if ( 2 == ntrk )
1460 if ( ( 5 == m_param.particle ) && ( 2 == ntrk ) && ( fabs( pTrk[0] ) < 5 ) &&
1461 ( fabs( pTrk[1] ) < 5 ) )
1462 {
1463 // if(1==ntrk) p = (event->getRecTrk(0)) -> getP();
1464 // if((5==m_param.particle) && (1==ntrk) && (fabs(p)<5)){
1465 m_tescos = tes;
1466 m_tesFlagcos = esTimeflag;
1467 for ( i = 0; i < ntrk; i++ )
1468 {
1469 rectrk = event->getRecTrk( i );
1470 phi0 = rectrk->getPhi0();
1472
1473 tanl = rectrk->getTanLamda();
1474 lamda = atan( tanl );
1475 theta = HFPI - lamda;
1476
1478 {
1479 m_nhitUpcos = rectrk->getNHits();
1480 m_pUpcos = rectrk->getP();
1481 m_ptUpcos = rectrk->getPt();
1482 m_phiUpcos = phi0;
1483 m_drUpcos = rectrk->getDr();
1484 m_dzUpcos = rectrk->getDz();
1485 m_ctheUpcos = cos( theta );
1486 }
1487 else
1488 {
1489 m_nhitDwcos = rectrk->getNHits();
1490 m_pDwcos = rectrk->getP();
1491 m_ptDwcos = rectrk->getPt();
1492 m_phiDwcos = phi0;
1493 m_drDwcos = rectrk->getDr();
1494 m_dzDwcos = rectrk->getDz();
1495 m_ctheDwcos = cos( theta );
1496
1497 if ( m_pDwcos > 0 ) m_chargecos = 1;
1498 else m_chargecos = 0;
1499 }
1500 }
1501 m_cosmic->write();
1502 }
1503
1504 if ( 1 == m_param.fgCalDetEffi ) calDetEffi();
1505
1506 return 1;
1507}
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
double sin(const BesAngle a)
double cos(const BesAngle a)
*************DOUBLE PRECISION m_pi *DOUBLE PRECISION m_HvecTau2 DOUBLE PRECISION m_HvClone2 DOUBLE PRECISION m_gamma1 DOUBLE PRECISION m_gamma2 DOUBLE PRECISION m_thet1 DOUBLE PRECISION m_thet2 INTEGER m_IFPHOT *COMMON c_Taupair $ !Spin Polarimeter vector first Tau $ !Spin Polarimeter vector second Tau $ !Clone Spin Polarimeter vector first Tau $ !Clone Spin Polarimeter vector second Tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !phi of HvecTau1 $ !theta of HvecTau1 $ !phi of HvecTau2 $ !theta of HvecTau2 $ !super key
Definition Taupair.h:42
std::ofstream ofstream
Referenced by GrXtMdcCalib::fillHist(), QtMdcCalib::fillHist(), T0MdcCalib::fillHist(), Wr2dMdcCalib::fillHist(), WrMdcCalib::fillHist(), XtInteMdcCalib::fillHist(), and XtMdcCalib::fillHist().
◆ initialize()
|
pure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 225 of file MdcCalib.cxx.
226 {
227 IMessageSvc* msgSvc;
230 log << MSG::INFO << "MdcCalib::initialize()" << endmsg;
231
232 m_hlist = hlist;
233 m_mdcGeomSvc = mdcGeomSvc;
234 m_mdcFunSvc = mdcFunSvc;
235 m_mdcUtilitySvc = mdcUtilitySvc;
236
237 int lay;
238 int iEntr;
239 int lr;
241 char hname[200];
242
243 int nbinMom = 2400;
244
245 m_nlayer = m_mdcGeomSvc->getLayerSize();
246
247 for ( lay = 0; lay < m_nlayer; lay++ )
248 { m_radii[lay] = m_mdcGeomSvc->Layer( lay )->Radius(); }
251 {
252 m_xe[wir] = m_mdcGeomSvc->Wire( wir )->Backward().x();
253 m_ye[wir] = m_mdcGeomSvc->Wire( wir )->Backward().y();
254 m_ze[wir] = m_mdcGeomSvc->Wire( wir )->Backward().z();
255 m_xw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().x();
256 m_yw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().y();
257 m_zw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().z();
258 fwpc << setw( 6 ) << wir << setw( 15 ) << m_xe[wir] << setw( 15 ) << m_ye[wir]
259 << setw( 15 ) << m_ze[wir] << setw( 15 ) << m_xw[wir] << setw( 15 ) << m_yw[wir]
260 << setw( 15 ) << m_zw[wir] << endl;
261 }
262 fwpc.close();
263
264 m_fdcom = new TFolder( "common", "common" );
265 m_hlist->Add( m_fdcom );
266
267 m_effNtrk = new TH1F( "effNtrk", "", 43, -0.5, 42.5 );
268 m_fdcom->Add( m_effNtrk );
269
270 m_effNtrkRecHit = new TH1F( "effNtrkRecHit", "", 43, -0.5, 42.5 );
271 m_fdcom->Add( m_effNtrkRecHit );
272
273 m_hresAllInc = new TH1F( "HResAllInc", "", 200, -1.0, 1.0 );
274 m_fdcom->Add( m_hresAllInc );
275
276 m_hresAllExc = new TH1F( "HResAllExc", "", 200, -1.0, 1.0 );
277 m_fdcom->Add( m_hresAllExc );
278
279 m_hresAllAve = new TH1F( "HResAllAve", "", 200, -1.0, 1.0 );
280 m_fdcom->Add( m_hresAllAve );
281
282 m_hresInnInc = new TH1F( "HResInnInc", "", 200, -1.0, 1.0 );
283 m_fdcom->Add( m_hresInnInc );
284
285 m_hresInnExc = new TH1F( "HResInnExc", "", 200, -1.0, 1.0 );
286 m_fdcom->Add( m_hresInnExc );
287
288 m_hresStpInc = new TH1F( "HResStpInc", "", 200, -1.0, 1.0 );
289 m_fdcom->Add( m_hresStpInc );
290
291 m_hresStpExc = new TH1F( "HResStpExc", "", 200, -1.0, 1.0 );
292 m_fdcom->Add( m_hresStpExc );
293
294 m_hresOutInc = new TH1F( "HResOutInc", "", 200, -1.0, 1.0 );
295 m_fdcom->Add( m_hresOutInc );
296
297 m_hresOutExc = new TH1F( "HResOutExc", "", 200, -1.0, 1.0 );
298 m_fdcom->Add( m_hresOutExc );
299
300 m_hresAllExcUp = new TH1F( "HResAllExcUp", "", 200, -1.0, 1.0 );
301 m_fdcom->Add( m_hresAllExcUp );
302
303 m_hresAllExcDw = new TH1F( "HResAllExcDw", "", 200, -1.0, 1.0 );
304 m_fdcom->Add( m_hresAllExcDw );
305
306 m_hresInnExcUp = new TH1F( "HResInnExcUp", "", 200, -1.0, 1.0 );
307 m_fdcom->Add( m_hresInnExcUp );
308
309 m_hresInnExcDw = new TH1F( "HResInnExcDw", "", 200, -1.0, 1.0 );
310 m_fdcom->Add( m_hresInnExcDw );
311
312 m_hresStpExcUp = new TH1F( "HResStpExcUp", "", 200, -1.0, 1.0 );
313 m_fdcom->Add( m_hresStpExcUp );
314
315 m_hresStpExcDw = new TH1F( "HResStpExcDw", "", 200, -1.0, 1.0 );
316 m_fdcom->Add( m_hresStpExcDw );
317
318 m_hresOutExcUp = new TH1F( "HResOutExcUp", "", 200, -1.0, 1.0 );
319 m_fdcom->Add( m_hresOutExcUp );
320
321 m_hresOutExcDw = new TH1F( "HResOutExcDw", "", 200, -1.0, 1.0 );
322 m_fdcom->Add( m_hresOutExcDw );
323
324 m_fdResQ = new TFolder( "ResQ", "ResQ" );
325 m_hlist->Add( m_fdResQ );
326 for ( int i = 0; i < 14; i++ )
327 {
329 m_hresAveAllQ[i] = new TH1F( hname, "", 200, -1, 1 );
330 m_fdResQ->Add( m_hresAveAllQ[i] );
331
333 m_hresAveOutQ[i] = new TH1F( hname, "", 200, -1, 1 );
334 m_fdResQ->Add( m_hresAveOutQ[i] );
335 }
336 for ( lay = 0; lay < 43; lay++ )
337 {
338 for ( int i = 0; i < 14; i++ )
339 {
341 m_hresAveLayQ[lay][i] = new TH1F( hname, "", 200, -1, 1 );
342 m_fdResQ->Add( m_hresAveLayQ[lay][i] );
343 }
344 }
345
346 for ( int iEs = 0; iEs < m_param.nEsFlag; iEs++ )
347 {
349 m_hTes[iEs] = new TH1F( hname, "", 750, 0, 1500 );
350 m_fdcom->Add( m_hTes[iEs] );
351 }
352
353 m_hbbTrkFlg = new TH1F( "BbTrkFlg", "", 100, 0, 6 );
354 m_fdcom->Add( m_hbbTrkFlg );
355
356 m_hTesAll = new TH1F( "TesAll", "", 2000, 0, 2000 );
357 m_fdcom->Add( m_hTesAll );
358
359 m_hTesGood = new TH1F( "TesGood", "", 2000, 0, 2000 );
360 m_fdcom->Add( m_hTesGood );
361
362 m_hTesAllFlag = new TH1F( "TesAllFlag", "", 300, -0.5, 299.5 );
363 m_fdcom->Add( m_hTesAllFlag );
364
365 m_hTesRec = new TH1F( "TesRec", "", 2000, 0, 2000 );
366 m_fdcom->Add( m_hTesRec );
367
368 m_hTesCalFlag = new TH1F( "TesCalFlag", "", 2000, 0, 2000 );
369 m_fdcom->Add( m_hTesCalFlag );
370
371 m_hTesCalUse = new TH1F( "TesCalUse", "", 2000, 0, 2000 );
372 m_fdcom->Add( m_hTesCalUse );
373
374 m_hnRawHit = new TH1F( "nRawHit", "", 6797, -0.5, 6796.5 );
375 m_fdcom->Add( m_hnRawHit );
376
377 m_hpt = new TH1F( "HPt", "", nbinMom, 0, 3 );
378 m_fdcom->Add( m_hpt );
379
380 m_hptPos = new TH1F( "HPtPos", "", nbinMom, 0, 3 );
381 m_fdcom->Add( m_hptPos );
382
383 m_hptNeg = new TH1F( "HPtNeg", "", nbinMom, 0, 3 );
384 m_fdcom->Add( m_hptNeg );
385
386 m_hp = new TH1F( "HP", "", nbinMom, 0, 3 );
387 m_fdcom->Add( m_hp );
388
389 m_hp_cms = new TH1F( "HPCMS", "", nbinMom, 0, 3 );
390 m_fdcom->Add( m_hp_cms );
391
392 m_hpMax = new TH1F( "HPMax", "", nbinMom, 0, 3 );
393 m_fdcom->Add( m_hpMax );
394
395 m_hpMaxCms = new TH1F( "HPMax_Cms", "", nbinMom, 0, 3 );
396 m_fdcom->Add( m_hpMaxCms );
397
398 m_hpPos = new TH1F( "HP_Pos", "", nbinMom, 0, 3 );
399 m_fdcom->Add( m_hpPos );
400
401 m_hpNeg = new TH1F( "HP_Neg", "", nbinMom, 0, 3 );
402 m_fdcom->Add( m_hpNeg );
403
404 m_hpPoscms = new TH1F( "HP_Pos_cms", "", nbinMom, 0, 3 );
405 m_fdcom->Add( m_hpPoscms );
406
407 m_hpNegcms = new TH1F( "HP_Neg_cms", "", nbinMom, 0, 3 );
408 m_fdcom->Add( m_hpNegcms );
409
410 m_hp_cut = new TH1F( "HPCut", "", nbinMom, 0, 3 );
411 m_fdcom->Add( m_hp_cut );
412
413 m_hchisq = new TH1F( "Chisq", "", 10, 0, 100 );
414 m_fdcom->Add( m_hchisq );
415
416 m_hnTrk = new TH1F( "HNtrack", "HNtrack", 10, -0.5, 9.5 );
417 m_fdcom->Add( m_hnTrk );
418
419 m_hnTrkCal = new TH1F( "HNtrackCal", "HNtrackCal", 10, -0.5, 9.5 );
420 m_fdcom->Add( m_hnTrkCal );
421
422 m_hnhitsRec = new TH1F( "HNhitsRec", "", 100, -0.5, 99.5 );
423 m_fdcom->Add( m_hnhitsRec );
424
425 m_hnhitsRecInn = new TH1F( "HNhitsInnRec", "", 60, 0.5, 60.5 );
426 m_fdcom->Add( m_hnhitsRecInn );
427
428 m_hnhitsRecStp = new TH1F( "HNhitsStpRec", "", 60, 0.5, 60.5 );
429 m_fdcom->Add( m_hnhitsRecStp );
430
431 m_hnhitsRecOut = new TH1F( "HNhitsOutRec", "", 60, 0.5, 60.5 );
432 m_fdcom->Add( m_hnhitsRecOut );
433
434 m_hnhitsCal = new TH1F( "HNhitsCal", "", 100, -0.5, 99.5 );
435 m_fdcom->Add( m_hnhitsCal );
436
437 m_hnhitsCalInn = new TH1F( "HNhitsCalInn", "", 60, 0.5, 60.5 );
438 m_fdcom->Add( m_hnhitsCalInn );
439
440 m_hnhitsCalStp = new TH1F( "HNhitsCalStp", "", 60, 0.5, 60.5 );
441 m_fdcom->Add( m_hnhitsCalStp );
442
443 m_hnhitsCalOut = new TH1F( "HNhitsCalOut", "", 60, 0.5, 60.5 );
444 m_fdcom->Add( m_hnhitsCalOut );
445
446 m_wirehitmap = new TH1F( "Wire_HitMap", "Wire_HitMap", 6796, -0.5, 6795.5 );
447 m_fdcom->Add( m_wirehitmap );
448
449 m_layerhitmap = new TH1F( "Layer_HitMap", "Layer_HitMap", 43, -0.5, 42.5 );
450 m_fdcom->Add( m_layerhitmap );
451
452 m_hnoisephi = new TH1F( "phi_noise", "", 100, 0, 6.284 );
453 m_fdcom->Add( m_hnoisephi );
454
455 m_hnoiselay = new TH1F( "Layer_noise", "Layer_noise", 43, -0.5, 42.5 );
456 m_fdcom->Add( m_hnoiselay );
457
458 m_hnoisenhits = new TH1F( "nhits_noise", "nhits_noise", 6796, -0.5, 6795.5 );
459 m_fdcom->Add( m_hnoisenhits );
460
461 m_hratio = new TH1F( "ratio", "", 100, 0, 1 );
462 m_fdcom->Add( m_hratio );
463
464 m_hdr = new TH1F( "dr", "", 2000, -100, 100 );
465 m_fdcom->Add( m_hdr );
466
467 m_hphi0 = new TH1F( "phi0", "", 100, 0, 6.284 );
468 m_fdcom->Add( m_hphi0 );
469
470 m_hkap = new TH1F( "kappa", "", 400, -50, 50 );
471 m_fdcom->Add( m_hkap );
472
473 m_hdz = new TH1F( "dz", "", 2000, -500, 500 );
474 m_fdcom->Add( m_hdz );
475
476 m_htanl = new TH1F( "tanl", "", 200, -5, 5 );
477 m_fdcom->Add( m_htanl );
478
479 m_hcosthe = new TH1F( "costheta", "", 200, -1, 1 );
480 m_fdcom->Add( m_hcosthe );
481
482 m_hcostheNeg = new TH1F( "costhetaNeg", "", 200, -1, 1 );
483 m_fdcom->Add( m_hcostheNeg );
484
485 m_hcosthePos = new TH1F( "costhetaPos", "", 200, -1, 1 );
486 m_fdcom->Add( m_hcosthePos );
487
488 m_hx0 = new TH1F( "x0", "", 100, -10, 10 );
489 m_fdcom->Add( m_hx0 );
490
491 m_hy0 = new TH1F( "y0", "", 100, -10, 10 );
492 m_fdcom->Add( m_hy0 );
493
494 m_hdelZ0 = new TH1F( "delta_z0", "", 100, -50, 50 );
495 m_fdcom->Add( m_hdelZ0 );
496
497 m_grX0Y0 = new TGraph();
498 m_grX0Y0->SetName( "x0y0" );
499 m_fdcom->Add( m_grX0Y0 );
500
501 m_hitEffAll = new TH1F( "hitEffAll", "", 6800, -0.5, 6799.5 );
502 m_fdcom->Add( m_hitEffAll );
503
504 m_hitEffRaw = new TH1F( "hitEffRaw", "", 6800, -0.5, 6799.5 );
505 m_fdcom->Add( m_hitEffRaw );
506
507 m_hitEffRec = new TH1F( "hitEffRec", "", 6800, -0.5, 6799.5 );
508 m_fdcom->Add( m_hitEffRec );
509
510 // histograms for drift time
511 m_fdTime = new TFolder( "time", "time" );
512 m_hlist->Add( m_fdTime );
513
514 for ( lay = 0; lay < m_nlayer; lay++ )
515 {
517 m_htraw[lay] = new TH1F( hname, "", 1000, 0, 2000 );
518 m_fdTime->Add( m_htraw[lay] );
519
521 m_htdr[lay] = new TH1F( hname, "", 510, -10, 500 );
522 m_fdTime->Add( m_htdr[lay] );
523
524 for ( lr = 0; lr < 2; lr++ )
525 {
527 m_htdrlr[lay][lr] = new TH1F( hname, "", 510, -10, 500 );
528 m_fdTime->Add( m_htdrlr[lay][lr] );
529 }
530 }
531
532 // histograms of adc
533 m_fdAdc = new TFolder( "adc", "adc" );
534 m_hlist->Add( m_fdAdc );
535
536 for ( lay = 0; lay < m_nlayer; lay++ )
537 {
539 m_hadc[lay] = new TH1F( hname, "", 1500, 0, 3000 );
540 m_fdAdc->Add( m_hadc[lay] );
541 }
542 // histograms for resolution
543 m_fdres = new TFolder( "resolution", "resolution" );
544 m_hlist->Add( m_fdres );
545
546 m_fdresAve = new TFolder( "resAve", "resAve" );
547 m_hlist->Add( m_fdresAve );
548
549 for ( lay = 0; lay < m_nlayer; lay++ )
550 {
552 m_hresInc[lay] = new TH1F( hname, "", 1000, -5, 5 );
553 m_fdres->Add( m_hresInc[lay] );
554
556 m_hresExc[lay] = new TH1F( hname, "", 1000, -5, 5 );
557 m_fdres->Add( m_hresExc[lay] );
558
560 m_hresAve[lay] = new TH1F( hname, "", 1000, -5, 5 );
561 m_fdresAve->Add( m_hresAve[lay] );
562
563 for ( lr = 0; lr < 2; lr++ )
564 {
566 // m_hreslrInc[lay][lr] = new TH1F(hname, "", 200, -1, 1);
567 m_hreslrInc[lay][lr] = new TH1F( hname, "", 1000, -5, 5 );
568 m_fdres->Add( m_hreslrInc[lay][lr] );
569
571 // m_hreslrExc[lay][lr] = new TH1F(hname, "", 200, -1, 1);
572 m_hreslrExc[lay][lr] = new TH1F( hname, "", 1000, -5, 5 );
573 m_fdres->Add( m_hreslrExc[lay][lr] );
574
576 // m_hreslrAve[lay][lr] = new TH1F(hname, "", 200, -1, 1);
577 m_hreslrAve[lay][lr] = new TH1F( hname, "", 1000, -5, 5 );
578 m_fdresAve->Add( m_hreslrAve[lay][lr] );
579 }
580 for ( int phi = 0; phi < 20; phi++ )
581 {
583 m_hresphi[lay][phi] = new TH1F( hname, "", 200, -1, 1 );
584 m_fdres->Add( m_hresphi[lay][phi] );
585 }
586 }
587
588 /* histograms for momentum vs phi */
589 m_fdmomPhi = new TFolder( "momPhi", "momPhi" );
590 m_hlist->Add( m_fdmomPhi );
591
592 int thbin;
594 {
597 m_fdmomPhi->Add( m_ppPhi[bin] );
598
601 m_fdmomPhi->Add( m_pnPhi[bin] );
602
605 m_fdmomPhi->Add( m_ppPhiCms[bin] );
606
609 m_fdmomPhi->Add( m_pnPhiCms[bin] );
610
611 for ( thbin = 0; thbin < NThetaBin; thbin++ )
612 {
615 m_fdmomPhi->Add( m_ppThePhi[thbin][bin] );
616
619 m_fdmomPhi->Add( m_pnThePhi[thbin][bin] );
620
623 m_fdmomPhi->Add( m_ppThePhiCms[thbin][bin] );
624
627 m_fdmomPhi->Add( m_pnThePhiCms[thbin][bin] );
628 }
629 }
630 for ( thbin = 0; thbin < NThetaBin; thbin++ )
631 {
633 m_ppThe[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
634 m_fdmomPhi->Add( m_ppThe[thbin] );
635
637 m_pnThe[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
638 m_fdmomPhi->Add( m_pnThe[thbin] );
639
641 m_ppTheCms[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
642 m_fdmomPhi->Add( m_ppTheCms[thbin] );
643
645 m_pnTheCms[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
646 m_fdmomPhi->Add( m_pnTheCms[thbin] );
647 }
648
649 // histograms for resolution vs distance
650 m_fdres2d = new TFolder( "res2d", "res2d" );
651 m_hlist->Add( m_fdres2d );
652
653 int hid = 0;
655 TH1F* hist;
656 for ( lay = 0; lay < m_nlayer; lay++ )
657 {
659 {
660 for ( lr = 0; lr < 2; lr++ )
661 {
663 {
665 hist = new TH1F( hname, "", 200, -1, 1 );
666 m_hr2dInc.push_back( hist );
667 m_fdres2d->Add( hist );
668
670 hist = new TH1F( hname, "", 200, -1, 1 );
671 m_hr2dExc.push_back( hist );
672 m_fdres2d->Add( hist );
673
676 hid++;
677 }
678 }
679 }
680 } // end of layer loop
681
682 m_fdres2t = new TFolder( "res2t", "res2t" );
683 // m_hlist -> Add(m_fdres2t);
684
685 for ( lay = 0; lay < m_nlayer; lay++ )
686 {
688 {
689 for ( lr = 0; lr < 2; lr++ )
690 {
692 {
695 m_fdres2t->Add( m_hr2t[lay][iEntr][lr][bin] );
696 }
697 }
698 }
699 }
700
701 INTupleSvc* ntupleSvc;
703 for ( lay = 0; lay < m_nlayer; lay++ )
704 {
706 NTuplePtr nt( ntupleSvc, hname );
707 if ( nt ) m_xtTuple[lay] = nt;
708 else
709 {
711 if ( m_xtTuple[lay] )
712 {
713 m_xtTuple[lay]->addItem( "cel", m_cel[lay] );
714 m_xtTuple[lay]->addItem( "lr", m_lr[lay] );
715 m_xtTuple[lay]->addItem( "run", m_run[lay] );
716 m_xtTuple[lay]->addItem( "evt", m_evt[lay] );
717 m_xtTuple[lay]->addItem( "doca", m_doca[lay] );
718 m_xtTuple[lay]->addItem( "dm", m_dm[lay] );
719 m_xtTuple[lay]->addItem( "tdr", m_tdr[lay] );
720 m_xtTuple[lay]->addItem( "tdc", m_tdc[lay] );
721 m_xtTuple[lay]->addItem( "entr", m_entr[lay] );
722 m_xtTuple[lay]->addItem( "zhit", m_zhit[lay] );
723 m_xtTuple[lay]->addItem( "qhit", m_qhit[lay] );
724 m_xtTuple[lay]->addItem( "p", m_p[lay] );
725 m_xtTuple[lay]->addItem( "pt", m_pt[lay] );
726 m_xtTuple[lay]->addItem( "phi0", m_phi0[lay] );
727 m_xtTuple[lay]->addItem( "tanl", m_tanl[lay] );
728 m_xtTuple[lay]->addItem( "hitphi", m_hitphi[lay] );
729 }
730 else
731 { log << MSG::FATAL << "Cannot book Xt N-tuple:" << long( m_xtTuple[lay] ) << endmsg; }
732 }
733 }
734
735 if ( 5 == m_param.particle )
736 {
738 NTuplePtr nt( ntupleSvc, hname );
739 if ( nt ) m_cosmic = nt;
740 else
741 {
743 if ( m_cosmic )
744 {
745 m_cosmic->addItem( "pUp", m_pUpcos );
746 m_cosmic->addItem( "pDw", m_pDwcos );
747 m_cosmic->addItem( "ptUp", m_ptUpcos );
748 m_cosmic->addItem( "ptDw", m_ptDwcos );
749 m_cosmic->addItem( "phiUp", m_phiUpcos );
750 m_cosmic->addItem( "phiDw", m_phiDwcos );
751 m_cosmic->addItem( "drUp", m_drUpcos );
752 m_cosmic->addItem( "drDw", m_drDwcos );
753 m_cosmic->addItem( "dzUp", m_dzUpcos );
754 m_cosmic->addItem( "dzDw", m_dzDwcos );
755 m_cosmic->addItem( "ctheUp", m_ctheUpcos );
756 m_cosmic->addItem( "ctheDw", m_ctheDwcos );
757 m_cosmic->addItem( "nhitUp", m_nhitUpcos );
758 m_cosmic->addItem( "nhitDw", m_nhitDwcos );
759 m_cosmic->addItem( "char", m_chargecos );
760 m_cosmic->addItem( "tesfg", m_tesFlagcos );
761 m_cosmic->addItem( "tes", m_tescos );
762 }
763 }
764 }
765}
Referenced by GrXtMdcCalib::initialize(), QtMdcCalib::initialize(), T0MdcCalib::initialize(), Wr2dMdcCalib::initialize(), WrMdcCalib::initialize(), XtInteMdcCalib::initialize(), and XtMdcCalib::initialize().
◆ printCut()
|
pure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 1509 of file MdcCalib.cxx.
1509 {
1510 cout << "Events " << m_hTesAll->GetEntries() << ", Tes cut " << m_hTesCalFlag->GetEntries()
1511 << ", nTrkCut " << m_cut1 << endl;
1512 cout << "Tot Tracks " << m_nTrkAfTrkCut << ", cos(theta)_cut " << m_cut2 << ", drCut "
1513 << m_cut3 << ", dzCut " << m_cut4 << ", momCut" << m_cut5 << ", nHitLayer_cut "
1514 << m_cut6 << ", nHit_cut " << m_cut7 << endl;
1515 cout << "nTrack after all cuts: " << m_nTrkCal << endl;
1516}
Referenced by GrXtMdcCalib::printCut(), QtMdcCalib::printCut(), T0MdcCalib::printCut(), Wr2dMdcCalib::printCut(), WrMdcCalib::printCut(), XtInteMdcCalib::printCut(), and XtMdcCalib::printCut().
◆ setParam()
|
inlinepure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 305 of file MdcCalib.h.
305{ m_param = param; }
Referenced by GrXtMdcCalib::setParam(), IniMdcCalib::setParam(), PreT0MdcCalib::setParam(), PreXtMdcCalib::setParam(), QtMdcCalib::setParam(), T0MdcCalib::setParam(), Wr2dMdcCalib::setParam(), WrMdcCalib::setParam(), XtInteMdcCalib::setParam(), and XtMdcCalib::setParam().
◆ updateConst()
|
pure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 1518 of file MdcCalib.cxx.
1518 {
1519 IMessageSvc* msgSvc;
1522 log << MSG::DEBUG << "MdcCalib::updateConst()" << endmsg;
1523
1524 int lay;
1525 double effi;
1526 double effErr;
1527
1528 int nGoodAll = 0;
1529 int nGoodInn = 0;
1530 int nGoodStp = 0;
1531 int nGoodOut = 0;
1532 int nTotAll = 0;
1533 int nTotInn = 0;
1534 int nTotStp = 0;
1535 int nTotOut = 0;
1537 for ( lay = 0; lay < m_nlayer; lay++ )
1538 {
1539 double effNtrk = m_effNtrk->GetBinContent( lay + 1 );
1540 double effGoodHit = m_effNtrkRecHit->GetBinContent( lay + 1 );
1541 nGoodAll += effGoodHit;
1542 if ( lay < 8 ) nGoodInn += effGoodHit;
1543 else if ( lay < 20 ) nGoodStp += effGoodHit;
1544 else nGoodOut += effGoodHit;
1545
1546 nTotAll += effNtrk;
1547 if ( lay < 8 ) nTotInn += effNtrk;
1548 else if ( lay < 20 ) nTotStp += effNtrk;
1549 else nTotOut += effNtrk;
1550
1551 effi = (double)effGoodHit / (double)effNtrk;
1552 effErr = sqrt( effi * ( 1 - effi ) / (double)effNtrk );
1553 feffi << setw( 5 ) << lay << setw( 15 ) << effi << setw( 15 ) << effErr << setw( 15 )
1554 << effGoodHit << setw( 15 ) << effNtrk << endl;
1555 }
1556 double effiAll = (double)nGoodAll / (double)( nTotAll );
1557 double errAll = sqrt( effiAll * ( 1 - effiAll ) / (double)( nTotAll ) );
1558 double effiInn = (double)nGoodInn / (double)( nTotInn );
1559 double errInn = sqrt( effiInn * ( 1 - effiInn ) / (double)( nTotInn ) );
1560 double effiStp = (double)nGoodStp / (double)( nTotStp );
1561 double errStp = sqrt( effiStp * ( 1 - effiStp ) / (double)( nTotStp ) );
1562 double effiOut = (double)nGoodOut / (double)( nTotOut );
1563 double errOut = sqrt( effiOut * ( 1 - effiOut ) / (double)( nTotOut ) );
1564 feffi << endl
1565 << "EffiAll: " << setw( 15 ) << effiAll << setw( 15 ) << errAll << setw( 15 )
1566 << nGoodAll << setw( 15 ) << nTotAll << endl;
1567 feffi << endl
1568 << "EffiInn: " << setw( 15 ) << effiInn << setw( 15 ) << errInn << setw( 15 )
1569 << nGoodInn << setw( 15 ) << nTotInn << endl;
1570 feffi << endl
1571 << "EffiStp: " << setw( 15 ) << effiStp << setw( 15 ) << errStp << setw( 15 )
1572 << nGoodStp << setw( 15 ) << nTotStp << endl;
1573 feffi << endl
1574 << "EffiOut: " << setw( 15 ) << effiOut << setw( 15 ) << errOut << setw( 15 )
1575 << nGoodOut << setw( 15 ) << nTotOut << endl;
1576 feffi.close();
1577
1578 // calculate efficiency without the impact of track fitting
1579 if ( 0 != m_param.fgCalDetEffi )
1580 {
1581 int nHitAll[] = { 0, 0 };
1582 int nHitInn[] = { 0, 0 };
1583 int nHitStp[] = { 0, 0 };
1584 int nHitOut[] = { 0, 0 };
1586 for ( lay = 0; lay < m_nlayer; lay++ )
1587 {
1588 nHitAll[0] += m_hitNum[lay][0];
1589 nHitAll[1] += m_hitNum[lay][1];
1590 if ( lay < 8 )
1591 {
1592 nHitInn[0] += m_hitNum[lay][0];
1593 nHitInn[1] += m_hitNum[lay][1];
1594 }
1595 else if ( lay < 20 )
1596 {
1597 nHitStp[0] += m_hitNum[lay][0];
1598 nHitStp[1] += m_hitNum[lay][1];
1599 }
1600 else
1601 {
1602 nHitOut[0] += m_hitNum[lay][0];
1603 nHitOut[1] += m_hitNum[lay][1];
1604 }
1605
1606 effi = (double)m_hitNum[lay][1] / (double)m_hitNum[lay][0];
1607 effErr = sqrt( effi * ( 1 - effi ) / (double)m_hitNum[lay][0] );
1608 feff2 << setw( 5 ) << lay << setw( 15 ) << effi << setw( 15 ) << effErr << setw( 15 )
1609 << m_hitNum[lay][1] << setw( 15 ) << m_hitNum[lay][0] << endl;
1610 }
1611 effiAll = (double)nHitAll[1] / (double)nHitAll[0];
1612 errAll = sqrt( effiAll * ( 1 - effiAll ) ) / (double)nHitAll[0];
1613 effiInn = (double)nHitInn[1] / (double)nHitInn[0];
1614 errInn = sqrt( effiInn * ( 1 - effiInn ) ) / (double)nHitInn[0];
1615 effiStp = (double)nHitStp[1] / (double)nHitStp[0];
1616 errStp = sqrt( effiStp * ( 1 - effiStp ) ) / (double)nHitStp[0];
1617 effiOut = (double)nHitOut[1] / (double)nHitOut[0];
1618 errOut = sqrt( effiOut * ( 1 - effiOut ) ) / (double)nHitOut[0];
1619 feff2 << endl
1620 << "EffiAll: " << setw( 15 ) << effiAll << setw( 15 ) << errAll << setw( 15 )
1621 << nHitAll[1] << setw( 15 ) << nHitAll[0] << endl;
1622 feff2 << endl
1623 << "EffiInn: " << setw( 15 ) << effiInn << setw( 15 ) << errInn << setw( 15 )
1624 << nHitInn[1] << setw( 15 ) << nHitInn[0] << endl;
1625 feff2 << endl
1626 << "EffiStp: " << setw( 15 ) << effiStp << setw( 15 ) << errStp << setw( 15 )
1627 << nHitStp[1] << setw( 15 ) << nHitStp[0] << endl;
1628 feff2 << endl
1629 << "EffiOut: " << setw( 15 ) << effiOut << setw( 15 ) << errOut << setw( 15 )
1630 << nHitOut[1] << setw( 15 ) << nHitOut[0] << endl;
1631 feff2.close();
1632 }
1633
1634 // get resolution
1635 int i;
1636 int iEntr;
1637 int lr;
1640 int hid;
1641
1642 Stat_t entry;
1644 if ( m_param.calSigma )
1645 {
1647 for ( lay = 0; lay < m_nlayer; lay++ )
1648 {
1649 for ( iEntr = 0; iEntr < m_nEntr[lay]; iEntr++ )
1650 {
1651 for ( lr = 0; lr < 2; lr++ )
1652 {
1653 fr2d << setw( 3 ) << lay << setw( 3 ) << iEntr << setw( 3 ) << lr << endl;
1655 {
1657 hid = m_mapr2d[key];
1658
1659 if ( 1 == m_param.resiType )
1660 {
1661 entry = m_hr2dExc[hid]->GetEntries();
1662 if ( entry > 500 )
1663 {
1664 m_hr2dExc[hid]->Fit( "gaus", "Q" );
1666 }
1669 }
1670 else
1671 {
1672 entry = m_hr2dInc[hid]->GetEntries();
1673 if ( entry > 500 )
1674 {
1675 m_hr2dInc[hid]->Fit( "gaus", "Q" );
1677 }
1680 }
1682 } // end of bin loop
1683
1685 { sigm[bin] = sigm[m_nBin[lay] - 1]; }
1686
1688 {
1689 if ( 1 == m_param.fgCalib[lay] )
1690 {
1691 // calconst -> resetSdpar(lay, iEntr, lr, bin,
1692 // sigm[bin]);
1693 if ( 1 == m_nEntr[lay] )
1694 {
1696 }
1697 else if ( 2 == m_nEntr[lay] )
1698 {
1699 if ( 0 == iEntr )
1700 {
1701 for ( i = 0; i < 3; i++ )
1702 { // entr<0
1704 }
1705 }
1706 else
1707 {
1708 for ( i = 3; i < 6; i++ )
1709 { // entr>0
1711 }
1712 }
1713 }
1714 }
1717 } // end of bin loop
1718 }
1719 } // end of entr loop
1720 }
1721 fr2d.close();
1722 }
1723
1724 return 1;
1725}
void resetSdpar(int lay, int entr, int lr, int bin, double val)
Definition MdcCalibConst.cpp:84
double getSdpar(int lay, int entr, int lr, int bin)
Definition MdcCalibConst.cpp:89
Referenced by GrXtMdcCalib::updateConst(), QtMdcCalib::updateConst(), T0MdcCalib::updateConst(), WrMdcCalib::updateConst(), XtInteMdcCalib::updateConst(), and XtMdcCalib::updateConst().
The documentation for this class was generated from the following files:
