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 | clear ()=0 |
Detailed Description
Definition at line 36 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_cut2=0;
40 m_cut3=0;
41 m_cut4=0;
42 m_cut5=0;
43 m_cut6=0;
44 m_nGrPoint = 0;
45 fgReadWireEff = false;
46
47 int lay;
48 for(lay=0; lay<43; lay++){
49 if(lay < 15) m_nEntr[lay] = 1;
50 else m_nEntr[lay] = 2;
51 }
52 m_dwid = 0.5; // mm
53 m_fgIni = false;
54
55 m_phiWid = PI2 / (double)NPhiBin;
56 m_theWid = 2.0 / (double)NThetaBin;
57
58 m_nEvtNtuple = 0;
59
61 if(lay < 8) m_nBin[lay] = 12;
62 else m_nBin[lay] = 16;
63 }
64
65 // setting boundary layer flags
67 if((0==lay) || (7==lay) || (8==lay) || (19==lay) || (20==lay) ||
68 (35==lay) || (36==lay) || (42==lay) ) m_layBound[lay] = true;
69 else m_layBound[lay] = false;
70 }
71}
◆ ~MdcCalib()
|
virtual |
Definition at line 73 of file MdcCalib.cxx.
73 {
74}
Member Function Documentation
◆ clear()
|
pure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 76 of file MdcCalib.cxx.
76 {
77 int lay;
78 for(lay=0; lay<m_nlayer; lay++){
79 delete m_htraw[lay];
80 delete m_htdr[lay];
81 delete m_hresInc[lay];
82 delete m_hresExc[lay];
83 delete m_hresAve[lay];
84 delete m_hadc[lay];
85 for (int lr=0; lr<2; lr++){
86 delete m_htdrlr[lay][lr];
87 delete m_hreslrInc[lay][lr];
88 delete m_hreslrExc[lay][lr];
89 delete m_hreslrAve[lay][lr];
90 }
91 }
92
93 delete m_effNtrk;
94 delete m_effNtrkRecHit;
95 delete m_hresAllInc;
96 delete m_hresAllExc;
97 delete m_hresAllAve;
98 for(int i=0; i<14; i++){
99 delete m_hresAveAllQ[i];
100 delete m_hresAveOutQ[i];
101 }
102 for(lay=0; lay<43; lay++){
103 for(int i=0; i<14; i++) delete m_hresAveLayQ[lay][i];
104 }
105 delete m_hresInnInc;
106 delete m_hresInnExc;
107 delete m_hresStpInc;
108 delete m_hresStpExc;
109 delete m_hresOutInc;
110 delete m_hresOutExc;
112 delete m_hbbTrkFlg;
113 delete m_hTesAll;
114 delete m_hTesGood;
115 delete m_hTesAllFlag;
116 delete m_hTesRec;
117 delete m_hTesCalFlag;
118 delete m_hTesCalUse;
119 delete m_hnRawHit;
120 delete m_hpt;
121 delete m_hpMax;
122 delete m_hpMaxCms;
123 delete m_hptPos;
124 delete m_hptNeg;
125 delete m_hp;
126 delete m_hp_cms;
127 delete m_hpPos;
128 delete m_hpNeg;
129 delete m_hpPoscms;
130 delete m_hpNegcms;
131 delete m_hp_cut;
132 delete m_hchisq;
133 delete m_hnTrk;
134 delete m_hnTrkCal;
135 delete m_hnhitsRec;
136 delete m_hnhitsRecInn;
137 delete m_hnhitsRecStp;
138 delete m_hnhitsRecOut;
139 delete m_hnhitsCal;
140 delete m_hnhitsCalInn;
141 delete m_hnhitsCalStp;
142 delete m_hnhitsCalOut;
143 delete m_wirehitmap;
144 delete m_layerhitmap;
145 delete m_hnoisephi;
146 delete m_hnoiselay;
147 delete m_hnoisenhits;
148 delete m_hratio;
149 delete m_hdr;
150 delete m_hphi0;
151 delete m_hkap;
152 delete m_hdz;
153 delete m_htanl;
154 delete m_hcosthe;
155 delete m_hcostheNeg;
156 delete m_hcosthePos;
157 delete m_hx0;
158 delete m_hy0;
159 delete m_hdelZ0;
160 delete m_grX0Y0;
161 delete m_hitEffAll;
162 delete m_hitEffRaw;
163 delete m_hitEffRec;
165 int thbin;
171 for(thbin=0; thbin<NThetaBin; thbin++){
176 }
177 }
178 for(thbin=0; thbin<NThetaBin; thbin++){
179 delete m_ppThe[thbin];
180 delete m_pnThe[thbin];
181 delete m_ppTheCms[thbin];
182 delete m_pnTheCms[thbin];
183 }
184
185 for(unsigned i=0; i<m_hr2dInc.size(); i++){
186 delete m_hr2dInc[i];
187 delete m_hr2dExc[i];
188 }
189 m_hr2dInc.clear();
190 m_hr2dExc.clear();
191 m_mapr2d.clear();
192
193 delete m_fdTime;
194 delete m_fdAdc;
195 delete m_fdres;
196 delete m_fdresAve;
197 delete m_fdres2d;
198 delete m_fdcom;
199 delete m_fdResQ;
200}
*******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(), XtMdcCalib::clear(), and MdcCalibAlg::finalize().
◆ fillHist()
|
pure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 692 of file MdcCalib.cxx.
692 {
693 IMessageSvc* msgSvc;
696 log << MSG::DEBUG << "MdcCalib::fillHist()" << endreq;
697
698 int i;
699 int k;
700 int lay;
701 int cel;
702 int wir;
703 int lr;
704 int stat;
705
706 int hid;
708 int iEntr;
710
711 int phiBin;
712 int phiBinCms;
713 int theBin;
714 int theBinCms;
715 double lamda;
716 double theta;
717
718 double qhit;
719 double traw;
720 double tdr;
721 double doca;
722 double resiInc;
723 double resiExc;
724 double entr;
725 double pt;
726 double p;
728 double chisq;
733 HepLorentzVector p4;
734
735 double dr;
736 double phi0;
737 double dz;
738 double kap;
739 double tanl;
740
741 double x0;
742 double y0;
743 double zminus = 9999.0;
744 double zplus = -9999.0;
745
746 double hitphi;
747 double philab;
748 double phicms;
749 double thetacms;
750 double costheCMS;
751
752 double dphi;
753 double wphi;
754 double xx;
755 double yy;
756 double rr;
757
758 int nhitlay;
760 bool fgTrk;
761
762 int ntrk = event -> getNTrk();
763 int nhitRec;
764 int nhitRecInn;
765 int nhitRecStp;
766 int nhitRecOut;
767 int nhitCal;
768 int nhitCalInn;
769 int nhitCalStp;
770 int nhitCalOut;
771 MdcCalRecTrk* rectrk;
772 MdcCalRecHit* rechit;
773
774 int fgAdd[43]; // for calculating layer efficiency
775
776 // read dead wire
777 if(!fgReadWireEff){
780 for(cel=0; cel<ncel; cel++){
782 if(eff > 0.5) m_fgGoodWire[lay][cel] = true;
783 else m_fgGoodWire[lay][cel] = false;
784 if(eff<0.9) cout << "dead channel: " << setw(5) << lay << setw(5) << cel << endl;
785 }
786 }
787 fgReadWireEff = true;
788 }
789
790 int nRawHit = event->getNRawHitTQ();
791 m_hnRawHit->Fill(nRawHit);
792
793 IDataProviderSvc* eventSvc = NULL;
794 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
795
796 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
797 if (!eventHeader) {
798 log << MSG::FATAL << "Could not find Event Header" << endreq;
799 return( StatusCode::FAILURE);
800 }
801 int iRun = eventHeader->runNumber();
802 int iEvt = eventHeader->eventNumber();
803
804 int esTimeflag = event->getEsFlag();
805 double tes = event->getTes();
806 bool esCutFg = event->getEsCutFlag();
807 int iEs = event->getNesCutFlag();
808 //calculate the efficiency of Bhabha event
809 if (-1 != esTimeflag) {
810 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc, "/Event/Recon/RecMdcTrackCol");
811 if(!newtrkCol){
812 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endreq;
813 return ( StatusCode::FAILURE );
814 }
815 int nGoodTrk = 0;
816 int icharge = 0;
817 Vp4 p4p;
818 Vp4 p4m;
819 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
820 for(; it_trk != newtrkCol->end(); it_trk++){
822 HepLorentzVector p4 = (*it_trk)->p4(mass);
823 icharge = (*it_trk)->charge();
824 if (icharge > 0) p4p.push_back(p4);
825 else p4m.push_back(p4);
826 }
827 if (1 == p4p.size() * p4m.size()){
828 double dang = p4p[0].vect().angle(p4m[0].vect());
829 m_hbbTrkFlg->Fill(1);
830 if (dang > 2.94) {
831 m_hbbTrkFlg->Fill(2);
832 }
833 }
834
835 }
836 m_hTesAll->Fill(tes);
837// cout << "tes " << tes << endl;
838 if (-1 != esTimeflag) m_hTesGood->Fill(tes);
839 m_hTesAllFlag->Fill(esTimeflag);
840 if(ntrk > 0) m_hTesRec->Fill(tes);
842 if( esCutFg ) m_hTesCalFlag->Fill(tes);
843 else return -1;
844
845 if(! m_fgIni){
849 }
850 m_fgIni = true;
851 }
852
853 bool trkFlag[200]; // for calculating hit efficiency without impact of track fitting
854 for(i=0; i<200; i++) trkFlag[i] = false;
855
856 int ntrkCal = 0;
857 double pTrk[2];
858 double pTrkcms[2];
859 double hitphiplus = 9999.0;
860 double hitthetaplus = 9999.0;
861 double hitphiminus = -9999.0;
862 double hitthetaminus = -9999.0;
863 Vp4 pp;
864 Vp4 pm;
865 pp.clear();
866 pm.clear();
867 int phibinp;
868 int phibinm;
869
870 m_hnTrk->Fill(ntrk);
872 m_cut1++;
873 return -1;
874 }
875// if(ntrk > 2) return -1;
876 for(i=0; i<ntrk; i++){
877 fgTrk = true;
878 rectrk = event -> getRecTrk(i);
879 nhitRec = rectrk -> getNHits();
880 m_hnhitsRec -> Fill( nhitRec );
881
883 fgHitLay[lay] = false;
884 }
885
886 nhitRecInn = 0;
887 nhitRecStp = 0;
888 nhitRecOut = 0;
889 for(k=0; k<nhitRec; k++){
890 rechit = rectrk -> getRecHit(k);
891 lay = rechit -> getLayid();
892 doca = rechit -> getDocaExc();
893 resiExc = rechit -> getResiExc();
894 fgHitLay[lay] = true;
895
896 if(lay < 8) nhitRecInn++;
897 else if(lay < 20) nhitRecStp++;
898 else nhitRecOut++;
899 }
900 m_hnhitsRecInn->Fill(nhitRecInn);
901 m_hnhitsRecStp->Fill(nhitRecStp);
902 m_hnhitsRecOut->Fill(nhitRecOut);
903
904 // get momentum
905 pt = rectrk -> getPt();
906 p = rectrk -> getP();
907
908 // boost P to the cms
909 p4 = rectrk->getP4();
910 HepLorentzVector psip(xboost, yboost, zboost, ecm);
911 Hep3Vector boostv = psip.boostVector();
912 p4.boost(- boostv);
913 p_cms = p4.rho();
914 phicms = p4.phi();
916 thetacms = p4.theta();
917 costheCMS = cos(thetacms);
919 p4.boost(boostv);
920// cout << setw(15) << p << setw(15) << p_cms << setw(15) << costheCMS << endl;
921
922 // cos(theta) cut
924 m_cut2++;
925 continue;
926 }
927
928 // dr cut
929 dr = rectrk->getDr();
931 m_cut3++;
932 continue;
933 }
934
935 // dz cut
936 dz = rectrk->getDz();
938 m_cut4++;
939 continue;
940 }
941
942// if(! fgTrk) continue;
943
944 // hit layer cut
945 nhitlay = 0;
947 if(fgHitLay[lay]) nhitlay++;
948 }
950 m_cut5++;
951 continue;
952 }
953
954 // nhit cut
956 m_cut6++;
957 continue;
958 }
959
960// bool fgNoise = rectrk->getFgNoiseRatio();
961// if(m_param.noiseCut && (!fgNoise)) continue;
962// cout << setw(10) << iRun << setw(15) << iEvt << setw(5) << fgNoise << endl;
963
964// if(! ((fgHitLay[0]||fgHitLay[1]) && (fgHitLay[41]||fgHitLay[42])) ){
965// continue;
966// }
967
968 // calculate cell on the track
969 int cellTrkPass[43];
970 bool fgPass = getCellTrkPass(event, i, cellTrkPass);
971 for(lay=0; lay<m_nlayer; lay++){
972 fgAdd[lay] = 0;
973// if((16==lay) || (18==lay) || (19==lay) || (41==lay)){ // hv2200 2009-3
974 if((15==lay) || (16==lay) || (18==lay) || (19==lay) || (40==lay) || (41==lay) || (42==lay)){
975 int iCell = cellTrkPass[lay];
976 if(fgPass && (iCell >= 0) && m_fgGoodWire[lay][iCell]) m_effNtrk->Fill(lay);
977 else fgAdd[lay] = 1;
978 } else{
979 m_effNtrk->Fill(lay);
980 }
981 }
982
983 chisq = rectrk -> getChisq();
984 m_hchisq -> Fill( chisq );
985
986 if(pt > 0){
987 m_hpt -> Fill(pt);
988 m_hptPos -> Fill(pt);
989 m_hp -> Fill(p);
990 m_hp_cms -> Fill(p_cms);
991 m_hpPos -> Fill(p);
992 m_hpPoscms -> Fill(p_cms);
993 } else{
994 m_hpt -> Fill(-1.0*pt);
995 m_hptNeg -> Fill(-1.0*pt);
996 m_hp -> Fill(-1.0*p);
997 m_hp_cms -> Fill(-1.0*p_cms);
998 m_hpNeg -> Fill(-1.0*p);
999 m_hpNegcms -> Fill(-1.0*p_cms);
1000 }
1001 if(2 == ntrk){
1002 pTrk[i] = fabs(p);
1003 pTrkcms[i] = fabs(p_cms);
1004 }
1005
1006 dr = rectrk -> getDr();
1007 phi0 = rectrk -> getPhi0();
1008 kap = rectrk -> getKappa();
1009 dz = rectrk -> getDz();
1010 tanl = rectrk -> getTanLamda();
1011 lamda = atan(tanl);
1012 theta = HFPI - lamda;
1013
1014 m_hdr -> Fill(dr);
1015 m_hphi0 -> Fill(phi0);
1016 m_hkap -> Fill(kap);
1017 m_hdz -> Fill(dz);
1018 m_htanl -> Fill(tanl);
1019 m_hcosthe -> Fill(cos(theta));
1022
1023 philab = phi0 + HFPI;
1025// cout << setw(15) << phi0 << setw(15) << philab << setw(15) << phicms << endl;
1026
1027 phiBin = (int)(philab / m_phiWid);
1028 phiBinCms = (int)(phicms / m_phiWid);
1029 theBin = (int)((cos(theta) + 1.0) / m_theWid);
1030 theBinCms = (int)((cos(thetacms) + 1.0) / m_theWid);
1031 if(phiBin < 0) phiBin = 0;
1032 if(phiBin >= NPhiBin) phiBin = NPhiBin-1;
1033 if(phiBinCms < 0) phiBinCms = 0;
1034 if(phiBinCms >= NPhiBin) phiBinCms = NPhiBin-1;
1035 if(theBin < 0) theBin = 0;
1036 if(theBin >= NThetaBin) theBin = NThetaBin-1;
1037 if(theBinCms < 0) theBinCms = 0;
1038 if(theBinCms >= NThetaBin) theBinCms = NThetaBin-1;
1039
1040 if(pt > 0){
1041 m_ppPhi[phiBin]->Fill(p);
1042 m_ppPhiCms[phiBinCms]->Fill(p_cms);
1043 m_ppThe[theBin]->Fill(p);
1044 m_ppTheCms[theBinCms]->Fill(p_cms);
1045 m_ppThePhi[theBin][phiBin]->Fill(p);
1046 m_ppThePhiCms[theBinCms][phiBinCms]->Fill(p_cms);
1047 } else{
1048 m_pnPhi[phiBin]->Fill(-1.0*p);
1049 m_pnPhiCms[phiBinCms]->Fill(-1.0*p_cms);
1050 m_pnThe[theBin]->Fill(-1.0*p);
1051 m_pnTheCms[theBinCms]->Fill(-1.0*p_cms);
1052 m_pnThePhi[theBin][phiBin]->Fill(-1.0*p);
1053 m_pnThePhiCms[theBinCms][phiBinCms]->Fill(-1.0*p_cms);
1054 }
1055
1056 x0 = dr * cos(phi0);
1057 y0 = dr * sin(phi0);
1058 m_hx0 -> Fill(x0);
1059 m_hy0 -> Fill(y0);
1060 if(m_nGrPoint < 10000){
1061 m_grX0Y0->SetPoint(m_nGrPoint, x0, y0);
1062 m_nGrPoint++;
1063 }
1064
1065 if(kap < 0) {
1066 zminus = dz;
1067 pm.push_back(p4);
1068 phibinm = phiBinCms;
1069 } else {
1070 zplus = dz;
1071 pp.push_back(p4);
1072 phibinp = phiBinCms;
1073 }
1074
1075// cout << "phi = " << setw(15) << philab << setw(15) << philab*180./3.14159 << setw(15) << p << endl;
1076 ntrkCal++;
1077 trkFlag[i] = true;
1078 nhitCal = 0;
1079 nhitCalInn = 0;
1080 nhitCalStp = 0;
1081 nhitCalOut = 0;
1082 for(k=0; k<nhitRec; k++){
1083 rechit = rectrk -> getRecHit(k);
1084
1085 lay = rechit -> getLayid();
1086 cel = rechit -> getCellid();
1087 lr = rechit -> getLR();
1088 stat = rechit -> getStat();
1089 doca = rechit -> getDocaExc();
1090 resiInc = rechit -> getResiIncLR();
1091 resiExc = rechit -> getResiExcLR();
1092 entr = rechit -> getEntra();
1093 tdr = rechit -> getTdrift();
1094 traw = (rechit -> getTdc()) * MdcCalTdcCnv;
1095 wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
1096
1097 m_cel[lay] = (long)cel;
1098 m_lr[lay] = (long)lr;
1099 m_run[lay] = iRun;
1100 m_evt[lay] = iEvt;
1101 m_doca[lay] = doca;
1102 m_dm[lay] = rechit -> getDmeas();
1103 m_tdr[lay] = tdr;
1104 m_tdc[lay] = traw;
1105 m_entr[lay] = entr*180.0/3.14;
1106 m_zhit[lay] = rechit -> getZhit();
1107 m_qhit[lay] = rechit -> getQhit();
1108 m_p[lay] = p;
1109 m_pt[lay] = pt;
1110 m_phi0[lay] = phi0;
1111 m_tanl[lay] = tanl;
1112 qhit = rechit -> getQhit();
1113
1114 // calculating hitphi
1115 xx = (m_zhit[lay] - m_zw[wir]) * (m_xe[wir] - m_xw[wir]) /
1116 (m_ze[wir] - m_zw[wir]) + m_xw[wir];
1117 yy = (m_zhit[lay] - m_zw[wir]) * (m_ye[wir] - m_yw[wir]) /
1118 (m_ze[wir] - m_zw[wir]) + m_yw[wir];
1119 rr = sqrt( (xx * xx) + (yy * yy) );
1120 dphi = fabs(doca) / m_radii[lay];
1121
1122 if( yy >= 0 ) wphi = acos(xx / rr);
1124 if(1 == lr) hitphi = wphi + dphi; // mention
1125 else hitphi = wphi - dphi;
1128
1129 m_hitphi[lay] = hitphi;
1130
1131 if( (fabs(doca) > m_docaMax[lay]) ||
1132 (fabs(resiExc) > m_param.resiCut[lay]) ){
1133 continue;
1134 }
1135
1137 if(0 == lay){
1138 if( ! fgHitLay[1] ) continue;
1139 } else if(42 == lay){
1140 if( ! fgHitLay[41] ) continue;
1141 } else{
1142 if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ) continue;
1143
1144 // for boundary layers
1146 ((!fgHitLay[lay-1]) || (!fgHitLay[lay+1])) ) continue;
1147 }
1148 }
1149
1151
1152 // fill xtplot tree
1154 m_xtTuple[lay] -> write();
1155 }
1156
1158 if( (0 == fgAdd[lay]) && (1 == stat) ){
1159 m_effNtrkRecHit->Fill(lay);
1160 fgAdd[lay] = 1;
1161 }
1162 } else{
1163 if(0 == fgAdd[lay]){
1164 m_effNtrkRecHit->Fill(lay);
1165 fgAdd[lay] = 1;
1166 }
1167 }
1168
1169 nhitCal++;
1170 if(lay < 8) nhitCalInn++;
1171 else if(lay < 20) nhitCalStp++;
1172 else nhitCalOut++;
1173
1174 m_wirehitmap -> Fill(wir);
1175 m_layerhitmap -> Fill(lay);
1176
1177 m_htraw[lay] -> Fill(traw);
1178 m_htdr[lay] -> Fill(tdr);
1179 m_htdrlr[lay][lr]->Fill(tdr);
1180 m_hadc[lay] -> Fill(qhit);
1181
1182 m_hresAllInc -> Fill(resiInc);
1183 m_hresAllExc -> Fill(resiExc);
1184 double resiAve = 0.5 * (resiInc + resiExc);
1185 m_hresAllAve -> Fill(resiAve);
1186
1187 if(lay < 8){
1188 m_hresInnInc -> Fill(resiInc);
1189 m_hresInnExc -> Fill(resiExc);
1190 } else if(lay < 20){
1191 m_hresStpInc -> Fill(resiInc);
1192 m_hresStpExc -> Fill(resiExc);
1193 } else{
1194 m_hresOutInc -> Fill(resiInc);
1195 m_hresOutExc -> Fill(resiExc);
1196 }
1197
1198 int qbin = (int)((qhit-100.0)/100.0);
1199 if(qbin>=0 && qbin<14){
1200 m_hresAveAllQ[qbin]->Fill(resiAve);
1201 m_hresAveLayQ[lay][qbin]->Fill(resiAve);
1202 if(lay > 7) m_hresAveOutQ[qbin]->Fill(resiAve);
1203 }
1204
1205 m_hresInc[lay] -> Fill(resiInc);
1206 m_hreslrInc[lay][lr]->Fill(resiInc);
1207 m_hresExc[lay] -> Fill(resiExc);
1208 m_hreslrExc[lay][lr]->Fill(resiExc);
1209 m_hresAve[lay] -> Fill(resiAve);
1210 m_hreslrAve[lay][lr]->Fill(resiAve);
1211
1213 if(iPhi>=20) iPhi = 19;
1214 m_hresphi[lay][iPhi]->Fill((resiInc+resiExc)*0.5);
1215
1216// bin = (int)(fabs(doca) / m_dwid);
1218 iEntr = m_mdcFunSvc -> getSdEntrIndex(entr);
1219 if(1 == m_nEntr[lay]){
1220 iEntr = 0;
1221 } else if(2 == m_nEntr[lay]){
1222 if(entr > 0.0) iEntr = 1;
1223 else iEntr = 0;
1224 }
1228 hid = m_mapr2d[key];
1229 m_hr2dInc[hid] -> Fill(resiInc);
1230 m_hr2dExc[hid] -> Fill(resiExc);
1231 }
1232 }
1233
1234 if((tdr>0) && (tdr<750)){
1237 m_hr2t[lay][iEntr][lr][bin]->Fill(resiExc);
1238 }
1239 } // loop of nhits
1240 m_nEvtNtuple++;
1241 m_hnhitsCal->Fill(nhitCal);
1242 m_hnhitsCalInn->Fill(nhitCalInn);
1243 m_hnhitsCalStp->Fill(nhitCalStp);
1244 m_hnhitsCalOut->Fill(nhitCalOut);
1245 } // end of track loop
1246 m_hnTrkCal->Fill(ntrkCal);
1247 if(2 == ntrkCal){
1248 if(pTrk[0] > pTrk[1]) m_hpMax->Fill(pTrk[0]);
1249 else m_hpMax->Fill(pTrk[1]);
1250
1251 if(pTrkcms[0] > pTrkcms[1]) m_hpMaxCms->Fill(pTrkcms[0]);
1252 else m_hpMaxCms->Fill(pTrkcms[1]);
1253 }
1254 if(ntrkCal > 0) m_hTesCalUse->Fill(tes);
1255
1256 double delZ0;
1257 if((fabs(zminus) < 9000.0) && (fabs(zplus) < 9000.0)) delZ0 = zplus - zminus;
1258 m_hdelZ0 -> Fill(delZ0);
1259
1260 if (1 == pp.size() * pm.size()){
1261 HepLorentzVector ptot = pp[0] + pm[0];
1262 bool fourmomcut = false;
1263// fourmomcut = (ptot.x()>0.02 && ptot.x()<0.06) && (fabs(ptot.y()) < 0.02)
1264// && (ptot.z()>-0.01 && ptot.z()<0.03) && (ptot.e()>3.4 && ptot.e()<4.0);
1265 fourmomcut = (fabs(ptot.x()-0.04)<0.026) && (fabs(ptot.y()) < 0.026)
1266 && (fabs(ptot.z()-0.005)<0.036) && (fabs(ptot.e()-ecm)<0.058);
1267 //cout << "x = " << ptot.x() << ", y = " << ptot.y() << ", z = " << ptot.z() << ", e = " << ptot.e() << endl;
1268 if (fourmomcut) {
1269 HepLorentzVector psip(xboost, yboost, zboost, ecm);
1270 Hep3Vector boostv = psip.boostVector();
1271 pp[0].boost(- boostv);
1272 pm[0].boost(- boostv);
1273 m_hp_cut->Fill(pp[0].rho());
1274 m_hp_cut->Fill(pm[0].rho());
1275 }
1276 }
1277
1278 if(2==ntrk) for(i=0; i<ntrk; i++) pTrk[i] = (event -> getRecTrk(i)) -> getP();
1280// if(1==ntrk) p = (event->getRecTrk(0)) -> getP();
1281// if((5==m_param.particle) && (1==ntrk) && (fabs(p)<5)){
1282 m_tescos = tes;
1283 m_tesFlagcos = esTimeflag;
1284 for(i=0; i<ntrk; i++){
1285 rectrk = event -> getRecTrk(i);
1286 phi0 = rectrk -> getPhi0();
1288
1289 tanl = rectrk -> getTanLamda();
1290 lamda = atan(tanl);
1291 theta = HFPI - lamda;
1292
1294 m_nhitUpcos = rectrk -> getNHits();
1295 m_pUpcos = rectrk -> getP();
1296 m_ptUpcos = rectrk -> getPt();
1297 m_phiUpcos = phi0;
1298 m_drUpcos = rectrk->getDr();
1299 m_dzUpcos = rectrk->getDz();
1300 m_ctheUpcos = cos(theta);
1301 } else{
1302 m_nhitDwcos = rectrk -> getNHits();
1303 m_pDwcos = rectrk -> getP();
1304 m_ptDwcos = rectrk -> getPt();
1305 m_phiDwcos = phi0;
1306 m_drDwcos = rectrk->getDr();
1307 m_dzDwcos = rectrk->getDz();
1308 m_ctheDwcos = cos(theta);
1309
1310 if(m_pDwcos > 0) m_chargecos = 1;
1311 else m_chargecos = 0;
1312 }
1313 }
1314 m_cosmic->write();
1315 }
1316
1318
1319 return 1;
1320}
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
*************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
virtual double getWireEff(int layid, int cellid) const =0
virtual const MdcGeoLayer *const Layer(unsigned id)=0
Definition MdcCalRecHit.h:18
Definition MdcCalRecTrk.h:19
int NCell(void) const
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 202 of file MdcCalib.cxx.
203 {
204 IMessageSvc* msgSvc;
207 log << MSG::INFO << "MdcCalib::initialize()" << endreq;
208
209 m_hlist = hlist;
210 m_mdcGeomSvc = mdcGeomSvc;
211 m_mdcFunSvc = mdcFunSvc;
212 m_mdcUtilitySvc = mdcUtilitySvc;
213
214 int lay;
215 int iEntr;
216 int lr;
218 char hname[200];
219
220 m_nlayer = m_mdcGeomSvc -> getLayerSize();
221
222 for(lay=0; lay<m_nlayer; lay++){
224 }
225 ofstream fwpc("wirelog.txt");
233 fwpc << setw(6) << wir << setw(15) << m_xe[wir] << setw(15) << m_ye[wir]
234 << setw(15) << m_ze[wir] << setw(15) << m_xw[wir]
235 << setw(15) << m_yw[wir] << setw(15) << m_zw[wir] << endl;
236 }
237 fwpc.close();
238
239 m_fdcom = new TFolder("common", "common");
240 m_hlist -> Add(m_fdcom);
241
242 m_effNtrk = new TH1F("effNtrk", "", 43, -0.5, 42.5);
243 m_fdcom->Add(m_effNtrk);
244
245 m_effNtrkRecHit = new TH1F("effNtrkRecHit", "", 43, -0.5, 42.5);
246 m_fdcom->Add(m_effNtrkRecHit);
247
248 m_hresAllInc = new TH1F("HResAllInc", "", 200, -1.0, 1.0);
249 m_fdcom -> Add(m_hresAllInc);
250
251 m_hresAllExc = new TH1F("HResAllExc", "", 200, -1.0, 1.0);
252 m_fdcom -> Add(m_hresAllExc);
253
254 m_hresAllAve = new TH1F("HResAllAve", "", 200, -1.0, 1.0);
255 m_fdcom -> Add(m_hresAllAve);
256
257 m_hresInnInc = new TH1F("HResInnInc", "", 200, -1.0, 1.0);
258 m_fdcom -> Add(m_hresInnInc);
259
260 m_hresInnExc = new TH1F("HResInnExc", "", 200, -1.0, 1.0);
261 m_fdcom -> Add(m_hresInnExc);
262
263 m_hresStpInc = new TH1F("HResStpInc", "", 200, -1.0, 1.0);
264 m_fdcom -> Add(m_hresStpInc);
265
266 m_hresStpExc = new TH1F("HResStpExc", "", 200, -1.0, 1.0);
267 m_fdcom -> Add(m_hresStpExc);
268
269 m_hresOutInc = new TH1F("HResOutInc", "", 200, -1.0, 1.0);
270 m_fdcom -> Add(m_hresOutInc);
271
272 m_hresOutExc = new TH1F("HResOutExc", "", 200, -1.0, 1.0);
273 m_fdcom -> Add(m_hresOutExc);
274
275 m_fdResQ = new TFolder("ResQ", "ResQ");
276 m_hlist->Add(m_fdResQ);
277 for(int i=0; i<14; i++){
278 sprintf(hname, "resoAll_qbin%02d", i);
279 m_hresAveAllQ[i] = new TH1F(hname, "", 200, -1, 1);
280 m_fdResQ->Add(m_hresAveAllQ[i]);
281
282 sprintf(hname, "resoOut_qbin%02d", i);
283 m_hresAveOutQ[i] = new TH1F(hname, "", 200, -1, 1);
284 m_fdResQ->Add(m_hresAveOutQ[i]);
285 }
286 for(lay=0; lay<43; lay++){
287 for(int i=0; i<14; i++){
288 sprintf(hname, "resoLay%02d_qbin%02d", lay, i);
289 m_hresAveLayQ[lay][i] = new TH1F(hname, "", 200, -1, 1);
290 m_fdResQ->Add(m_hresAveLayQ[lay][i]);
291 }
292 }
293
296 m_hTes[iEs] = new TH1F(hname, "", 750, 0, 1500);
297 m_fdcom->Add(m_hTes[iEs]);
298 }
299
300 m_hbbTrkFlg = new TH1F("BbTrkFlg", "", 100, 0, 6);
301 m_fdcom -> Add(m_hbbTrkFlg);
302
303 m_hTesAll = new TH1F("TesAll", "", 1000, 0, 2000);
304 m_fdcom -> Add(m_hTesAll);
305
306 m_hTesGood = new TH1F("TesGood", "", 1000, 0, 2000);
307 m_fdcom -> Add(m_hTesGood);
308
309 m_hTesAllFlag = new TH1F("TesAllFlag", "", 300, -0.5, 299.5);
310 m_fdcom -> Add(m_hTesAllFlag);
311
312 m_hTesRec = new TH1F("TesRec", "", 1000, 0, 2000);
313 m_fdcom -> Add(m_hTesRec);
314
315 m_hTesCalFlag = new TH1F("TesCalFlag", "", 1000, 0, 2000);
316 m_fdcom -> Add(m_hTesCalFlag);
317
318 m_hTesCalUse = new TH1F("TesCalUse", "", 1000, 0, 2000);
319 m_fdcom -> Add(m_hTesCalUse);
320
321 m_hnRawHit = new TH1F("nRawHit", "", 6797, -0.5, 6796.5);
322 m_fdcom -> Add(m_hnRawHit);
323
324 m_hpt = new TH1F("HPt", "", 800, 0, 3);
325 m_fdcom -> Add(m_hpt);
326
327 m_hptPos = new TH1F("HPtPos", "", 800, 0, 3);
328 m_fdcom -> Add(m_hptPos);
329
330 m_hptNeg = new TH1F("HPtNeg", "", 800, 0, 3);
331 m_fdcom -> Add(m_hptNeg);
332
333 m_hp = new TH1F("HP", "", 800, 0, 3);
334 m_fdcom -> Add(m_hp);
335
336 m_hp_cms = new TH1F("HPCMS", "", 800, 0, 3);
337 m_fdcom -> Add(m_hp_cms);
338
339 m_hpMax = new TH1F("HPMax", "", 800, 0, 3);
340 m_fdcom -> Add(m_hpMax);
341
342 m_hpMaxCms = new TH1F("HPMax_Cms", "", 800, 0, 3);
343 m_fdcom -> Add(m_hpMaxCms);
344
345 m_hpPos = new TH1F("HP_Pos", "", 800, 0, 3);
346 m_fdcom -> Add(m_hpPos);
347
348 m_hpNeg = new TH1F("HP_Neg", "", 800, 0, 3);
349 m_fdcom -> Add(m_hpNeg);
350
351 m_hpPoscms = new TH1F("HP_Pos_cms", "", 800, 0, 3);
352 m_fdcom -> Add(m_hpPoscms);
353
354 m_hpNegcms = new TH1F("HP_Neg_cms", "", 800, 0, 3);
355 m_fdcom -> Add(m_hpNegcms);
356
357 m_hp_cut = new TH1F("HPCut", "", 800, 0, 3);
358 m_fdcom -> Add(m_hp_cut);
359
360 m_hchisq = new TH1F("Chisq", "", 10, 0, 100);
361 m_fdcom -> Add(m_hchisq);
362
363 m_hnTrk = new TH1F("HNtrack", "HNtrack", 10, -0.5, 9.5);
364 m_fdcom -> Add(m_hnTrk);
365
366 m_hnTrkCal = new TH1F("HNtrackCal", "HNtrackCal", 10, -0.5, 9.5);
367 m_fdcom -> Add(m_hnTrkCal);
368
369 m_hnhitsRec = new TH1F("HNhitsRec", "", 100, -0.5, 99.5);
370 m_fdcom -> Add(m_hnhitsRec);
371
372 m_hnhitsRecInn = new TH1F("HNhitsInnRec", "", 60, 0.5, 60.5);
373 m_fdcom -> Add(m_hnhitsRecInn);
374
375 m_hnhitsRecStp = new TH1F("HNhitsStpRec", "", 60, 0.5, 60.5);
376 m_fdcom -> Add(m_hnhitsRecStp);
377
378 m_hnhitsRecOut = new TH1F("HNhitsOutRec", "", 60, 0.5, 60.5);
379 m_fdcom -> Add(m_hnhitsRecOut);
380
381 m_hnhitsCal = new TH1F("HNhitsCal", "", 100, -0.5, 99.5);
382 m_fdcom -> Add(m_hnhitsCal);
383
384 m_hnhitsCalInn = new TH1F("HNhitsCalInn", "", 60, 0.5, 60.5);
385 m_fdcom -> Add(m_hnhitsCalInn);
386
387 m_hnhitsCalStp = new TH1F("HNhitsCalStp", "", 60, 0.5, 60.5);
388 m_fdcom -> Add(m_hnhitsCalStp);
389
390 m_hnhitsCalOut = new TH1F("HNhitsCalOut", "", 60, 0.5, 60.5);
391 m_fdcom -> Add(m_hnhitsCalOut);
392
393 m_wirehitmap = new TH1F("Wire_HitMap", "Wire_HitMap", 6796, -0.5, 6795.5);
394 m_fdcom -> Add(m_wirehitmap);
395
396 m_layerhitmap = new TH1F("Layer_HitMap", "Layer_HitMap", 43, -0.5, 42.5);
397 m_fdcom -> Add(m_layerhitmap);
398
399 m_hnoisephi = new TH1F("phi_noise", "", 100, 0, 6.284);
400 m_fdcom -> Add(m_hnoisephi);
401
402 m_hnoiselay = new TH1F("Layer_noise", "Layer_noise", 43, -0.5, 42.5);
403 m_fdcom -> Add(m_hnoiselay);
404
405 m_hnoisenhits = new TH1F("nhits_noise", "nhits_noise", 6796, -0.5, 6795.5);
406 m_fdcom -> Add(m_hnoisenhits);
407
408 m_hratio = new TH1F("ratio", "", 100, 0, 1);
409 m_fdcom -> Add(m_hratio);
410
411 m_hdr = new TH1F("dr", "", 500, -500, 500);
412 m_fdcom -> Add(m_hdr);
413
414 m_hphi0 = new TH1F("phi0", "", 100, 0, 6.284);
415 m_fdcom -> Add(m_hphi0);
416
417 m_hkap = new TH1F("kappa", "", 400, -50, 50);
418 m_fdcom -> Add(m_hkap);
419
420 m_hdz = new TH1F("dz", "", 500, -1000, 1000);
421 m_fdcom -> Add(m_hdz);
422
423 m_htanl = new TH1F("tanl", "", 200, -5, 5);
424 m_fdcom -> Add(m_htanl);
425
426 m_hcosthe = new TH1F("costheta", "", 200, -1, 1);
427 m_fdcom -> Add(m_hcosthe);
428
429 m_hcostheNeg = new TH1F("costhetaNeg", "", 200, -1, 1);
430 m_fdcom -> Add(m_hcostheNeg);
431
432 m_hcosthePos = new TH1F("costhetaPos", "", 200, -1, 1);
433 m_fdcom -> Add(m_hcosthePos);
434
435 m_hx0 = new TH1F("x0", "", 100, -10, 10);
436 m_fdcom -> Add(m_hx0);
437
438 m_hy0 = new TH1F("y0", "", 100, -10, 10);
439 m_fdcom -> Add(m_hy0);
440
441 m_hdelZ0 = new TH1F("delta_z0", "", 100, -50, 50);
442 m_fdcom -> Add(m_hdelZ0);
443
444 m_grX0Y0 = new TGraph();
445 m_grX0Y0->SetName("x0y0");
446 m_fdcom -> Add(m_grX0Y0);
447
448 m_hitEffAll = new TH1F("hitEffAll", "", 6800, -0.5, 6799.5);
449 m_fdcom -> Add(m_hitEffAll);
450
451 m_hitEffRaw = new TH1F("hitEffRaw", "", 6800, -0.5, 6799.5);
452 m_fdcom -> Add(m_hitEffRaw);
453
454 m_hitEffRec = new TH1F("hitEffRec", "", 6800, -0.5, 6799.5);
455 m_fdcom -> Add(m_hitEffRec);
456
457 // histograms for drift time
458 m_fdTime = new TFolder("time", "time");
459 m_hlist -> Add(m_fdTime);
460
461 for(lay=0; lay<m_nlayer; lay++){
462 sprintf(hname, "Traw%02d", lay);
463 m_htraw[lay] = new TH1F(hname, "", 1000, 0, 2000);
464 m_fdTime -> Add(m_htraw[lay]);
465
466 sprintf(hname, "Tdr%02d", lay);
467 m_htdr[lay] = new TH1F(hname, "", 510, -10, 500);
468 m_fdTime -> Add(m_htdr[lay]);
469
470 for (lr=0; lr<2; lr++){
471 sprintf(hname, "Tdr%02d_lr%01d", lay, lr);
472 m_htdrlr[lay][lr] = new TH1F(hname, "", 510, -10, 500);
473 m_fdTime -> Add(m_htdrlr[lay][lr]);
474 }
475 }
476
477 // histograms of adc
478 m_fdAdc = new TFolder("adc", "adc");
479 m_hlist -> Add(m_fdAdc);
480
481 for(lay=0; lay<m_nlayer; lay++){
482 sprintf(hname, "adc%02d", lay);
483 m_hadc[lay] = new TH1F(hname, "", 1500, 0, 3000);
484 m_fdAdc -> Add(m_hadc[lay]);
485 }
486 // histograms for resolution
487 m_fdres = new TFolder("resolution", "resolution");
488 m_hlist -> Add(m_fdres);
489
490 m_fdresAve = new TFolder("resAve", "resAve");
491 m_hlist -> Add(m_fdresAve);
492
493 for(lay=0; lay<m_nlayer; lay++){
494 sprintf(hname, "Reso%02dInc", lay);
495 m_hresInc[lay] = new TH1F(hname, "", 1000, -5, 5);
496 m_fdres -> Add(m_hresInc[lay]);
497
498 sprintf(hname, "Reso%02dExc", lay);
499 m_hresExc[lay] = new TH1F(hname, "", 1000, -5, 5);
500 m_fdres -> Add(m_hresExc[lay]);
501
502 sprintf(hname, "Reso%02d", lay);
503 m_hresAve[lay] = new TH1F(hname, "", 1000, -5, 5);
504 m_fdresAve -> Add(m_hresAve[lay]);
505
506 for (lr=0; lr<2; lr++){
507 sprintf(hname, "Reso%02dInc_lr%01d", lay, lr);
508// m_hreslrInc[lay][lr] = new TH1F(hname, "", 200, -1, 1);
509 m_hreslrInc[lay][lr] = new TH1F(hname, "", 1000, -5, 5);
510 m_fdres->Add(m_hreslrInc[lay][lr]);
511
512 sprintf(hname, "Reso%02dExc_lr%01d", lay, lr);
513// m_hreslrExc[lay][lr] = new TH1F(hname, "", 200, -1, 1);
514 m_hreslrExc[lay][lr] = new TH1F(hname, "", 1000, -5, 5);
515 m_fdres->Add(m_hreslrExc[lay][lr]);
516
517 sprintf(hname, "Reso%02d_lr%01d", lay, lr);
518// m_hreslrAve[lay][lr] = new TH1F(hname, "", 200, -1, 1);
519 m_hreslrAve[lay][lr] = new TH1F(hname, "", 1000, -5, 5);
520 m_fdresAve->Add(m_hreslrAve[lay][lr]);
521 }
522 for(int phi=0; phi<20; phi++){
523 sprintf(hname, "ResoPhi%02d_phi%02d", lay, phi);
524 m_hresphi[lay][phi] = new TH1F(hname, "", 200, -1, 1);
525 m_fdres->Add(m_hresphi[lay][phi]);
526 }
527 }
528
529 /* histograms for momentum vs phi */
530 m_fdmomPhi = new TFolder("momPhi", "momPhi");
531 m_hlist -> Add(m_fdmomPhi);
532
533 int thbin;
537 m_fdmomPhi->Add(m_ppPhi[bin]);
538
541 m_fdmomPhi->Add(m_pnPhi[bin]);
542
545 m_fdmomPhi->Add(m_ppPhiCms[bin]);
546
549 m_fdmomPhi->Add(m_pnPhiCms[bin]);
550
551 for(thbin=0; thbin<NThetaBin; thbin++){
554 m_fdmomPhi->Add(m_ppThePhi[thbin][bin]);
555
558 m_fdmomPhi->Add(m_pnThePhi[thbin][bin]);
559
562 m_fdmomPhi->Add(m_ppThePhiCms[thbin][bin]);
563
566 m_fdmomPhi->Add(m_pnThePhiCms[thbin][bin]);
567 }
568 }
569 for(thbin=0; thbin<NThetaBin; thbin++){
570 sprintf(hname, "hPpos_the%02d", thbin);
571 m_ppThe[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
572 m_fdmomPhi->Add(m_ppThe[thbin]);
573
574 sprintf(hname, "hPneg_the%02d", thbin);
575 m_pnThe[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
576 m_fdmomPhi->Add(m_pnThe[thbin]);
577
578 sprintf(hname, "hPposCms_the%02d", thbin);
579 m_ppTheCms[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
580 m_fdmomPhi->Add(m_ppTheCms[thbin]);
581
582 sprintf(hname, "hPnegCms_the%02d", thbin);
583 m_pnTheCms[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
584 m_fdmomPhi->Add(m_pnTheCms[thbin]);
585 }
586
587 // histograms for resolution vs distance
588 m_fdres2d = new TFolder("res2d", "res2d");
589 m_hlist -> Add(m_fdres2d);
590
591 int hid = 0;
593 TH1F* hist;
594 for(lay=0; lay<m_nlayer; lay++){
596 for(lr=0; lr<2; lr++){
599 hist = new TH1F(hname, "", 200, -1, 1);
600 m_hr2dInc.push_back(hist);
601 m_fdres2d -> Add(hist);
602
604 hist = new TH1F(hname, "", 200, -1, 1);
605 m_hr2dExc.push_back(hist);
606 m_fdres2d -> Add(hist);
607
610 hid++;
611 }
612 }
613 }
614 } // end of layer loop
615
616 m_fdres2t = new TFolder("res2t", "res2t");
617// m_hlist -> Add(m_fdres2t);
618
619 for(lay=0; lay<m_nlayer; lay++){
621 for(lr=0; lr<2; lr++){
625 m_fdres2t -> Add(m_hr2t[lay][iEntr][lr][bin]);
626 }
627 }
628 }
629 }
630
631 INTupleSvc* ntupleSvc;
633 for(lay=0; lay<m_nlayer; lay++){
634 sprintf(hname, "FILE136/xt%02d", lay);
635 NTuplePtr nt(ntupleSvc, hname);
636 if ( nt ) m_xtTuple[lay] = nt;
637 else{
639 if( m_xtTuple[lay] ){
640 m_xtTuple[lay]->addItem("cel", m_cel[lay]);
641 m_xtTuple[lay]->addItem("lr", m_lr[lay]);
642 m_xtTuple[lay]->addItem("run", m_run[lay]);
643 m_xtTuple[lay]->addItem("evt", m_evt[lay]);
644 m_xtTuple[lay]->addItem("doca", m_doca[lay]);
645 m_xtTuple[lay]->addItem("dm", m_dm[lay]);
646 m_xtTuple[lay]->addItem("tdr", m_tdr[lay]);
647 m_xtTuple[lay]->addItem("tdc", m_tdc[lay]);
648 m_xtTuple[lay]->addItem("entr", m_entr[lay]);
649 m_xtTuple[lay]->addItem("zhit", m_zhit[lay]);
650 m_xtTuple[lay]->addItem("qhit", m_qhit[lay]);
651 m_xtTuple[lay]->addItem("p", m_p[lay]);
652 m_xtTuple[lay]->addItem("pt", m_pt[lay]);
653 m_xtTuple[lay]->addItem("phi0", m_phi0[lay]);
654 m_xtTuple[lay]->addItem("tanl", m_tanl[lay]);
655 m_xtTuple[lay]->addItem("hitphi", m_hitphi[lay]);
656 } else{
657 log << MSG::FATAL << "Cannot book Xt N-tuple:"
658 << long(m_xtTuple[lay]) << endreq;
659 }
660 }
661 }
662
664 sprintf(hname, "FILE136/cosmic");
665 NTuplePtr nt(ntupleSvc, hname);
666 if ( nt ) m_cosmic = nt;
667 else{
669 if( m_cosmic ){
670 m_cosmic->addItem("pUp", m_pUpcos);
671 m_cosmic->addItem("pDw", m_pDwcos);
672 m_cosmic->addItem("ptUp", m_ptUpcos);
673 m_cosmic->addItem("ptDw", m_ptDwcos);
674 m_cosmic->addItem("phiUp", m_phiUpcos);
675 m_cosmic->addItem("phiDw", m_phiDwcos);
676 m_cosmic->addItem("drUp", m_drUpcos);
677 m_cosmic->addItem("drDw", m_drDwcos);
678 m_cosmic->addItem("dzUp", m_dzUpcos);
679 m_cosmic->addItem("dzDw", m_dzDwcos);
680 m_cosmic->addItem("ctheUp", m_ctheUpcos);
681 m_cosmic->addItem("ctheDw", m_ctheDwcos);
682 m_cosmic->addItem("nhitUp", m_nhitUpcos);
683 m_cosmic->addItem("nhitDw", m_nhitDwcos);
684 m_cosmic->addItem("char", m_chargecos);
685 m_cosmic->addItem("tesfg", m_tesFlagcos);
686 m_cosmic->addItem("tes", m_tescos);
687 }
688 }
689 }
690}
virtual const MdcGeoWire *const Wire(unsigned id)=0
double Radius(void) const
HepPoint3D Forward(void) const
HepPoint3D Backward(void) const
Referenced by GrXtMdcCalib::initialize(), QtMdcCalib::initialize(), T0MdcCalib::initialize(), Wr2dMdcCalib::initialize(), WrMdcCalib::initialize(), XtInteMdcCalib::initialize(), and XtMdcCalib::initialize().
◆ setParam()
|
inlinepure virtual |
Implemented in GrXtMdcCalib, IniMdcCalib, PreT0MdcCalib, PreXtMdcCalib, QtMdcCalib, T0MdcCalib, Wr2dMdcCalib, WrMdcCalib, XtInteMdcCalib, and XtMdcCalib.
Definition at line 293 of file MdcCalib.h.
293 {
294 m_param = param;
295}
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 1322 of file MdcCalib.cxx.
1322 {
1323 IMessageSvc* msgSvc;
1326 log << MSG::DEBUG << "MdcCalib::updateConst()" << endreq;
1327
1328 cout << "Tot " << m_hnTrk->GetEntries()
1329 << ", nTrkCut " << m_cut1 << ", cos(theta)_cut " << m_cut2 << ", drCut " << m_cut3
1330 << ", dzCut " << m_cut4 << ", nHitLayer_cut " << m_cut5 << ", nHit_cut " << m_cut6 << endl;
1331
1332 int lay;
1333 double effi;
1334 double effErr;
1335
1336 int nGoodAll = 0;
1337 int nGoodInn = 0;
1338 int nGoodStp = 0;
1339 int nGoodOut = 0;
1340 int nTotAll = 0;
1341 int nTotInn = 0;
1342 int nTotStp = 0;
1343 int nTotOut = 0;
1344 ofstream feffi("MdcLayerEffi.dat");
1345 for(lay=0; lay<m_nlayer; lay++){
1346 double effNtrk = m_effNtrk->GetBinContent(lay+1);
1347 double effGoodHit = m_effNtrkRecHit->GetBinContent(lay+1);
1348 nGoodAll += effGoodHit;
1349 if(lay < 8) nGoodInn += effGoodHit;
1350 else if (lay < 20) nGoodStp += effGoodHit;
1351 else nGoodOut += effGoodHit;
1352
1353 nTotAll += effNtrk;
1354 if(lay < 8) nTotInn += effNtrk;
1355 else if (lay < 20) nTotStp += effNtrk;
1356 else nTotOut += effNtrk;
1357
1358 effi = (double)effGoodHit / (double)effNtrk;
1359 effErr = sqrt(effi * (1-effi) / (double)effNtrk);
1360 feffi << setw(5) << lay << setw(15) << effi << setw(15) << effErr
1361 << setw(15) << effGoodHit << setw(15) << effNtrk << endl;
1362 }
1363 double effiAll = (double)nGoodAll / (double)(nTotAll);
1364 double errAll = sqrt(effiAll * (1-effiAll) / (double)(nTotAll));
1365 double effiInn = (double)nGoodInn / (double)(nTotInn);
1366 double errInn = sqrt(effiInn * (1-effiInn) / (double)(nTotInn));
1367 double effiStp = (double)nGoodStp / (double)(nTotStp);
1368 double errStp = sqrt(effiStp * (1-effiStp) / (double)(nTotStp));
1369 double effiOut = (double)nGoodOut / (double)(nTotOut);
1370 double errOut = sqrt(effiOut * (1-effiOut) / (double)(nTotOut));
1371 feffi << endl << "EffiAll: " << setw(15) << effiAll << setw(15) << errAll
1372 << setw(15) << nGoodAll << setw(15) << nTotAll << endl;
1373 feffi << endl << "EffiInn: " << setw(15) << effiInn << setw(15) << errInn
1374 << setw(15) << nGoodInn << setw(15) << nTotInn << endl;
1375 feffi << endl << "EffiStp: " << setw(15) << effiStp << setw(15) << errStp
1376 << setw(15) << nGoodStp << setw(15) << nTotStp << endl;
1377 feffi << endl << "EffiOut: " << setw(15) << effiOut << setw(15) << errOut
1378 << setw(15) << nGoodOut << setw(15) << nTotOut << endl;
1379 feffi.close();
1380
1381 // calculate efficiency without the impact of track fitting
1383 int nHitAll[] = {0, 0};
1384 int nHitInn[] = {0, 0};
1385 int nHitStp[] = {0, 0};
1386 int nHitOut[] = {0, 0};
1387 ofstream feff2("MdcHitEffi.dat");
1388 for(lay=0; lay<m_nlayer; lay++){
1389 nHitAll[0] += m_hitNum[lay][0];
1390 nHitAll[1] += m_hitNum[lay][1];
1391 if(lay < 8){
1392 nHitInn[0] += m_hitNum[lay][0];
1393 nHitInn[1] += m_hitNum[lay][1];
1394 } else if (lay < 20){
1395 nHitStp[0] += m_hitNum[lay][0];
1396 nHitStp[1] += m_hitNum[lay][1];
1397 } else{
1398 nHitOut[0] += m_hitNum[lay][0];
1399 nHitOut[1] += m_hitNum[lay][1];
1400 }
1401
1402 effi = (double)m_hitNum[lay][1] / (double)m_hitNum[lay][0];
1403 effErr = sqrt(effi * (1-effi) / (double)m_hitNum[lay][0]);
1404 feff2 << setw(5) << lay << setw(15) << effi << setw(15) << effErr
1405 << setw(15) << m_hitNum[lay][1] << setw(15) << m_hitNum[lay][0] << endl;
1406 }
1407 effiAll = (double)nHitAll[1] / (double)nHitAll[0];
1408 errAll = sqrt(effiAll * (1-effiAll)) / (double)nHitAll[0];
1409 effiInn = (double)nHitInn[1] / (double)nHitInn[0];
1410 errInn = sqrt(effiInn * (1-effiInn)) / (double)nHitInn[0];
1411 effiStp = (double)nHitStp[1] / (double)nHitStp[0];
1412 errStp = sqrt(effiStp * (1-effiStp)) / (double)nHitStp[0];
1413 effiOut = (double)nHitOut[1] / (double)nHitOut[0];
1414 errOut = sqrt(effiOut * (1-effiOut)) / (double)nHitOut[0];
1415 feff2 << endl << "EffiAll: " << setw(15) << effiAll << setw(15) << errAll
1416 << setw(15) << nHitAll[1] << setw(15) << nHitAll[0] << endl;
1417 feff2 << endl << "EffiInn: " << setw(15) << effiInn << setw(15) << errInn
1418 << setw(15) << nHitInn[1] << setw(15) << nHitInn[0] << endl;
1419 feff2 << endl << "EffiStp: " << setw(15) << effiStp << setw(15) << errStp
1420 << setw(15) << nHitStp[1] << setw(15) << nHitStp[0] << endl;
1421 feff2 << endl << "EffiOut: " << setw(15) << effiOut << setw(15) << errOut
1422 << setw(15) << nHitOut[1] << setw(15) << nHitOut[0] << endl;
1423 feff2.close();
1424 }
1425
1426 // get resolution
1427 int i;
1428 int iEntr;
1429 int lr;
1432 int hid;
1433
1434 Stat_t entry;
1437 ofstream fr2d("logr2d.dat");
1438 for(lay=0; lay<m_nlayer; lay++){
1439 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
1440 for(lr=0; lr<2; lr++){
1441 fr2d << setw(3) << lay << setw(3) << iEntr << setw(3) << lr << endl;
1444 hid = m_mapr2d[key];
1445
1447 entry = m_hr2dExc[hid] -> GetEntries();
1448 if(entry > 500){
1451 } else if(entry > 100){
1452 sigm[bin] = m_hr2dExc[hid]->GetRMS();
1453 } else{
1454 sigm[bin] = 0.2;
1455 }
1456 } else{
1457 entry = m_hr2dInc[hid] -> GetEntries();
1458 if(entry > 500){
1461 } else if(entry > 100){
1462 sigm[bin] = m_hr2dInc[hid]->GetRMS();
1463 } else{
1464 sigm[bin] = 0.2;
1465 }
1466 }
1468 } // end of bin loop
1469
1471 sigm[bin] = sigm[m_nBin[lay]-1];
1472 }
1473
1476// calconst -> resetSdpar(lay, iEntr, lr, bin, sigm[bin]);
1477 if(1 == m_nEntr[lay]){
1479 } else if(2 == m_nEntr[lay]){
1480 if(0 == iEntr){
1481 for(i=0; i<3; i++){ // entr<0
1483 }
1484 } else{
1485 for(i=3; i<6; i++){ // entr>0
1487 }
1488 }
1489 }
1490 } else{
1492 }
1494 } // end of bin loop
1495 }
1496 } // end of entr loop
1497 }
1498 fr2d.close();
1499 }
1500
1501 return 1;
1502}
gr Fit("g1","R")
double getSdpar(int lay, int entr, int lr, int bin)
Definition MdcCalibConst.cpp:105
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:
