00001 #include "CalDetDST/BetterMCTables.h"
00002
00003 #include <cmath>
00004 #include <iostream>
00005
00006 #include "DatabaseInterface/DbiCascader.h"
00007 #include "DatabaseInterface/DbiWriter.h"
00008 #include "Plex/PlexStripEndId.h"
00009 #include "Calibrator/Calibrator.h"
00010 #include "Calibrator/CalMIPCalibration.h"
00011 #include "Validity/VldTimeStamp.h"
00012 #include "Validity/VldContext.h"
00013 #include "Validity/VldRange.h"
00014 #include "UgliGeometry/UgliGeomHandle.h"
00015 #include "UgliGeometry/UgliStripHandle.h"
00016
00017 using std::cout;
00018 using std::endl;
00019
00020 float calc_acor(float lc, float lg, float c1, float c2, float l1, float l2, float l3){
00021
00022 float a=(1.0/(c1+c2))*(c1*std::exp(-lg/l1) + c2*std::exp(-lg/l2));
00023 float b=std::exp(-lc/l3);
00024 return (a*b);
00025 }
00026
00027 void load_mc_mipcal(float scale){
00028
00029
00030
00031
00032 DbiCascader& cascader=DbiTableProxyRegistry::Instance().GetCascader();
00033
00034
00035 string tableDescr ="(SEQNO int, SEIDKEY int, STRIPENDID int, SCALE float)";
00036
00037 Int_t dbNoTemp = cascader.CreateTemporaryTable("CALMIPCALIBRATION",tableDescr);
00038
00039 if (dbNoTemp<0){
00040 cout << "No database to will accept temporary tables. " << endl;
00041 return;
00042 }
00043 else{
00044 cout<<"Database number "<<dbNoTemp<<endl;
00045 }
00046 VldTimeStamp vts(2001,1,1,0,0,0);
00047 VldTimeStamp vte(2010,1,1,0,0,0);
00048
00049 VldRange vr(Detector::kCalDet,SimFlag::kMC,vts,vte,"by hand");
00050 VldTimeStamp create;
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 DbiWriter<CalMIPCalibration> tcal(vr, -1, 0, create);
00064
00065
00066 for(int plane=0;plane<60;plane++){
00067 for(int strip=0;strip<24;strip++){
00068 PlexStripEndId seidn(Detector::kCalDet,plane,strip,
00069 StripEnd::kNegative);
00070 PlexStripEndId seidp(Detector::kCalDet,plane,strip,
00071 StripEnd::kPositive);
00072 Int_t seidkeyn = seidn.BuildPlnStripEndKey();
00073 Int_t seidkeyp = seidp.BuildPlnStripEndKey();
00074 CalMIPCalibration tcn(seidkeyn, StripEnd::kNegative,scale);
00075 tcal<<tcn;
00076 CalMIPCalibration tcp(seidkeyp, StripEnd::kPositive,scale);
00077 tcal<<tcp;
00078 }
00079 }
00080 if(!tcal.CanOutput()){
00081 cout<<"MakeTimeCal: Writer can't output!"<<endl;
00082 }
00083
00084 cout<<"Done."<<endl;
00085 tcal.Close();
00086
00087 }
00088
00089
00090 void load_better_mipcal(float scale, VldTimeStamp vtm){
00091
00092
00093
00094
00095 DbiCascader& cascader=DbiTableProxyRegistry::Instance().GetCascader();
00096
00097
00098 string tableDescr ="(SEQNO int, SEIDKEY int, STRIPENDID int, SCALE float)";
00099
00100 Int_t dbNoTemp = cascader.CreateTemporaryTable("CALMIPCALIBRATION",tableDescr);
00101
00102 if (dbNoTemp<0){
00103 cout << "No database to will accept temporary tables. " << endl;
00104 return;
00105 }
00106 else{
00107 cout<<"Database number "<<dbNoTemp<<endl;
00108 }
00109
00110 VldTimeStamp vts(2001,1,1,0,0,0);
00111 VldTimeStamp vte(2010,1,1,0,0,0);
00112
00113
00114 VldRange vr(Detector::kCalDet,SimFlag::kMC,vts,vte,"by hand");
00115 VldTimeStamp create;
00116
00117 DbiWriter<CalMIPCalibration> tcal(vr, -1, 0, create);
00118
00119
00120 VldContext vcm(Detector::kCalDet, SimFlag::kMC, vtm);
00121
00122 cout<<"better_mc_tables! \nstart: "<<vts
00123 <<"\nend: "<<vte
00124 <<"\n vtm: "<<vtm
00125 <<"\ncontext: "<<vcm<<endl;
00126
00127
00128 UgliGeomHandle ugh(vcm);
00129
00130
00131
00132 Calibrator::Instance().Reset(vcm);
00133 const float about_1pe=60;
00134
00135
00136 for(int plane=0;plane<60;plane++){
00137 for(int strip=0;strip<24;strip++){
00138 PlexStripEndId seidn(Detector::kCalDet,plane,strip,
00139 StripEnd::kNegative);
00140 PlexStripEndId seidp(Detector::kCalDet,plane,strip,
00141 StripEnd::kPositive);
00142 const Int_t seidkeyn = seidn.BuildPlnStripEndKey();
00143 const Int_t seidkeyp = seidp.BuildPlnStripEndKey();
00144 UgliStripHandle ugs_n = ugh.GetStripHandle(seidn);
00145 UgliStripHandle ugs_p = ugh.GetStripHandle(seidp);
00146
00147 const bool nvalid = ugs_n.IsValid();
00148 const bool pvalid = ugs_p.IsValid();
00149 if(!nvalid) cout<<"neg not valid!"<<endl;
00150 if(!pvalid) cout<<"pos not valid!"<<endl;
00151 if(!(pvalid & nvalid)) continue;
00152
00153 const float clear_n = ugs_n.ClearFiber(seidn.GetEnd());
00154 const float green_n = ugs_n.GetHalfLength()
00155 +ugs_n.WlsPigtail(seidn.GetEnd());
00156
00157 const float clear_p = ugs_p.ClearFiber(seidp.GetEnd());
00158 const float green_p = ugs_p.GetHalfLength()
00159 +ugs_p.WlsPigtail(seidp.GetEnd());
00160
00161
00162 float calc_acor(float lc, float lg, float c1=0.3333, float c2=0.6667,
00163 float l1=1.052, float l2=7.078, float l3=11.4);
00164
00165
00166 const float an = calc_acor(clear_n, green_n,
00167 0.3333, 0.6667, 1.052, 7.078, 11.4);
00168 const float ap = calc_acor(clear_p, green_p,
00169 0.3333, 0.6667, 1.052, 7.078, 11.4);
00170
00171
00172
00173
00174 const float gn = about_1pe/(Calibrator::Instance().GetPhotoElectrons(about_1pe,seidn));
00175 const float gp = about_1pe/(Calibrator::Instance().GetPhotoElectrons(about_1pe,seidp));
00176
00177
00178
00179
00180
00181
00182
00183 const float norm=(calc_acor(0, 4.65)+calc_acor(6,0.65))/2.0;
00184
00185 const float sn=scale*gn*an/norm;
00186 const float sp=scale*gp*ap/norm;
00187
00188 cout<<plane<<" "<<strip<<" "<<sn<<" "<<sp<<endl;
00189
00190
00191 CalMIPCalibration tcn(seidkeyn, StripEnd::kNegative,sn);
00192 tcal<<tcn;
00193 CalMIPCalibration tcp(seidkeyp, StripEnd::kPositive,sp);
00194 tcal<<tcp;
00195 }
00196 }
00197 if(!tcal.CanOutput()){
00198 cout<<"MakeTimeCal: Writer can't output!"<<endl;
00199 }
00200
00201 cout<<"Done."<<endl;
00202 tcal.Close();
00203 cout<<"vtm: "<<vtm<<endl;
00204 }