00001 #include <memory>
00002 #include <string>
00003 #include <iostream>
00004 #include <fstream>
00005 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
00006 #include "CalibTracker/SiStripLorentzAngle/interface/SiStripCalibLorentzAngle.h"
00007 #include "FWCore/Framework/interface/ESHandle.h"
00008 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
00009 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00010 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00011 #include "FWCore/ServiceRegistry/interface/Service.h"
00012 #include "CLHEP/Random/RandGauss.h"
00013 #include "CondFormats/DataRecord/interface/SiStripLorentzAngleRcd.h"
00014 #include "DQM/SiStripCommon/interface/SiStripHistoId.h"
00015 #include "DQMServices/Core/interface/MonitorElement.h"
00016 #include "DQMServices/Core/interface/DQMStore.h"
00017 #include "Geometry/CommonTopologies/interface/StripTopology.h"
00018 #include "DQM/SiStripCommon/interface/ExtractTObject.h"
00019 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00020 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00021 #include <TF1.h>
00022 #include <TProfile.h>
00023
00024
00025 SiStripCalibLorentzAngle::SiStripCalibLorentzAngle(edm::ParameterSet const& conf) : ConditionDBWriter<SiStripLorentzAngle>(conf) , conf_(conf){}
00026
00027
00028 void SiStripCalibLorentzAngle::algoBeginJob(const edm::EventSetup& c){
00029
00030
00031 c.get<TrackerDigiGeometryRecord>().get(estracker);
00032 const TrackerGeometry *tracker=&(*estracker);
00033
00034 c.get<IdealMagneticFieldRecord>().get(magfield_);
00035
00036 edm::ESHandle<SiStripLorentzAngle> SiStripLorentzAngle_;
00037 c.get<SiStripLorentzAngleRcd>().get(SiStripLorentzAngle_);
00038 detid_la= SiStripLorentzAngle_->getLorentzAngles();
00039
00040 DQMStore* dbe_ = edm::Service<DQMStore>().operator->();
00041 std::string inputFile_ =conf_.getUntrackedParameter<std::string>("fileName", "LorentzAngle.root");
00042 dbe_->open(inputFile_);
00043
00044 SiStripHistoId hidmanager;
00045
00046 TF1 *fitfunc=0;
00047 double ModuleRangeMin=conf_.getParameter<double>("ModuleFitXMin");
00048 double ModuleRangeMax=conf_.getParameter<double>("ModuleFitXMax");
00049
00050 fitfunc= new TF1("fitfunc","([4]/[3])*[1]*(TMath::Abs(x-[0]))+[2]",-1,1);
00051
00052 std::vector<MonitorElement*> histolist= dbe_->getAllContents("/");
00053 std::vector<MonitorElement*>::iterator histo;
00054 LocalPoint p =LocalPoint(0,0,0);
00055
00056 dbe_->setCurrentFolder("LorentzAngle_Plots");
00057
00058 MonitorElement * LA_plot=dbe_->book1D("TanLAPerTesla","TanLAPerTesla",2000,-0.5,0.5);
00059 MonitorElement * LA_err_plot=dbe_->book1D("TanLAPerTesla Error","TanLAPerTesla Error",1000,0,0.2);
00060 MonitorElement * LA_chi2norm_plot=dbe_->book1D("TanLAPerTesla Chi2norm","TanLAPerTesla Chi2norm",2000,0,20);
00061 MonitorElement * LA_phi_plot=dbe_->bookProfile("TanLAPerTesla_vs_Phi","TanLAPerTesla_vs_Phi",200,-3.5,3.5,1000,-0.5,0.5,"");
00062 MonitorElement * LA_eta_plot=dbe_->bookProfile("TanLAPerTesla_vs_Eta","TanLAPerTesla_vs_Eta",200,-2.6,2.6,1000,-0.5,0.5,"");
00063
00064 MonitorElement * LA_plot_TIB=dbe_->book1D("TanLAPerTesla_TIB","TanLAPerTesla_TIB",2000,-0.5,0.5);
00065 MonitorElement * LA_plot_TIB_1_mono=dbe_->book1D("TanLAPerTesla_TIB_1_MONO","TanLAPerTesla_TIB_1_MONO",2000,-0.5,0.5);
00066 MonitorElement * LA_plot_TIB_1_stereo=dbe_->book1D("TanLAPerTesla_TIB_1_STEREO","TanLAPerTesla_TIB_1_STEREO",2000,-0.5,0.5);
00067 MonitorElement * LA_plot_TIB_2_mono=dbe_->book1D("TanLAPerTesla_TIB_2_MONO","TanLAPerTesla_TIB_2_MONO",2000,-0.5,0.5);
00068 MonitorElement * LA_plot_TIB_2_stereo=dbe_->book1D("TanLAPerTesla_TIB_2_STEREO","TanLAPerTesla_TIB_2_STEREO",2000,-0.5,0.5);
00069 MonitorElement * LA_plot_TIB_3=dbe_->book1D("TanLAPerTesla_TIB_3","TanLAPerTesla_TIB_3",2000,-0.5,0.5);
00070 MonitorElement * LA_plot_TIB_4=dbe_->book1D("TanLAPerTesla_TIB_4","TanLAPerTesla_TIB_4",2000,-0.5,0.5);
00071
00072 MonitorElement * LA_plot_TOB=dbe_->book1D("TanLAPerTesla_TOB","TanLAPerTesla_TOB",2000,-0.5,0.5);
00073 MonitorElement * LA_plot_TOB_1_mono=dbe_->book1D("TanLAPerTesla_TOB_1_MONO","TanLAPerTesla_TOB_1_MONO",2000,-0.5,0.5);
00074 MonitorElement * LA_plot_TOB_1_stereo=dbe_->book1D("TanLAPerTesla_TOB_1_STEREO","TanLAPerTesla_TOB_1_STEREO",2000,-0.5,0.5);
00075 MonitorElement * LA_plot_TOB_2_mono=dbe_->book1D("TanLAPerTesla_TOB_2_MONO","TanLAPerTesla_TOB_2_MONO",2000,-0.5,0.5);
00076 MonitorElement * LA_plot_TOB_2_stereo=dbe_->book1D("TanLAPerTesla_TOB_2_STEREO","TanLAPerTesla_TOB_2_STEREO",2000,-0.5,0.5);
00077 MonitorElement * LA_plot_TOB_3=dbe_->book1D("TanLAPerTesla_TOB_3","TanLAPerTesla_TOB_3",2000,-0.5,0.5);
00078 MonitorElement * LA_plot_TOB_4=dbe_->book1D("TanLAPerTesla_TOB_4","TanLAPerTesla_TOB_4",2000,-0.5,0.5);
00079 MonitorElement * LA_plot_TOB_5=dbe_->book1D("TanLAPerTesla_TOB_5","TanLAPerTesla_TOB_5",2000,-0.5,0.5);
00080 MonitorElement * LA_plot_TOB_6=dbe_->book1D("TanLAPerTesla_TOB_6","TanLAPerTesla_TOB_6",2000,-0.5,0.5);
00081
00082
00083 histos[1] = LA_plot;
00084 histos[2] = LA_err_plot;
00085 histos[3] = LA_chi2norm_plot;
00086 histos[4] = LA_phi_plot;
00087 histos[5] = LA_eta_plot;
00088
00089 histos[6] = LA_plot_TIB;
00090 histos[7] = LA_plot_TIB_1_mono;
00091 histos[8] = LA_plot_TIB_1_stereo;
00092 histos[9] = LA_plot_TIB_2_mono;
00093 histos[10] = LA_plot_TIB_2_stereo;
00094 histos[11] = LA_plot_TIB_3;
00095 histos[12] = LA_plot_TIB_4;
00096
00097 histos[13] = LA_plot_TOB;
00098 histos[14] = LA_plot_TOB_1_mono;
00099 histos[15] = LA_plot_TOB_1_stereo;
00100 histos[16] = LA_plot_TOB_2_mono;
00101 histos[17] = LA_plot_TOB_2_stereo;
00102 histos[18] = LA_plot_TOB_3;
00103 histos[19] = LA_plot_TOB_4;
00104 histos[20] = LA_plot_TOB_5;
00105 histos[21] = LA_plot_TOB_6;
00106
00107 gphi=-99;
00108 geta=-99;
00109
00110 int fitcounter = 0;
00111 int histocounter = 0;
00112 int badfit = 0;
00113 int no_mod_histo = 0;
00114 int Layer = 0;
00115
00116 ofstream LA_neg;
00117 LA_neg.open("Negative_LA.txt");
00118
00119 for(histo=histolist.begin();histo!=histolist.end();++histo){
00120
00121 if((*histo)->getEntries()==0){
00122 uint32_t id0=hidmanager.getComponentId((*histo)->getName());
00123 DetId detid0(id0);
00124 StripSubdetector subid0(id0);
00125 edm::LogInfo("SiStripCalibLorentzAngle")<<"### HISTOGRAM WITH 0 ENTRIES => TYPE:"<<subid0.subdetId();}
00126
00127 if((*histo)->getEntries()>100){
00128 histocounter++;
00129
00130 uint32_t id=hidmanager.getComponentId((*histo)->getName());
00131 DetId detid(id);
00132 StripSubdetector subid(id);
00133
00134 const GeomDetUnit * stripdet;
00135 if(!(stripdet=tracker->idToDetUnit(subid))){
00136 no_mod_histo++;
00137 edm::LogInfo("SiStripCalibLorentzAngle")<<"### NO MODULE HISTOGRAM";}
00138
00139 if(subid.subdetId() == int (StripSubdetector::TIB)){
00140 TIBDetId TIBid=TIBDetId(subid);
00141 Layer = TIBid.layer();
00142 edm::LogInfo("SiStripCalibLorentzAngle")<<"TIB layer = "<<Layer;}
00143
00144 if(subid.subdetId() == int (StripSubdetector::TOB)){
00145 TOBDetId TOBid=TIBDetId(subid);
00146 Layer = TOBid.layer();
00147 edm::LogInfo("SiStripCalibLorentzAngle")<<"TOB layer = "<<Layer;}
00148
00149 edm::LogInfo("SiStripCalibLorentzAngle")<<"id: "<<id;
00150
00151 if(stripdet!=0){
00152 const GlobalPoint gposition = (stripdet->surface()).toGlobal(p);
00153 gphi = gposition.phi();
00154 geta = gposition.eta();}
00155 if(stripdet==0)continue;
00156
00157 float thickness=stripdet->specificSurface().bounds().thickness();
00158 const StripTopology& topol=(StripTopology&)stripdet->topology();
00159 float pitch = topol.localPitch(p);
00160
00161 TProfile* theProfile=ExtractTObject<TProfile>().extract(*histo);
00162
00163 fitfunc->SetParameter(0, 0);
00164 fitfunc->SetParameter(1, 0);
00165 fitfunc->SetParameter(2, 1);
00166 fitfunc->FixParameter(3, pitch);
00167 fitfunc->FixParameter(4, thickness);
00168 int fitresult=theProfile->Fit(fitfunc,"N","",ModuleRangeMin, ModuleRangeMax);
00169
00170 if(fitfunc->GetParameter(0)>0){
00171 LA_neg<<"Type = "<<subid.subdetId()<<" Layer = "<<Layer;
00172 if(subid.subdetId() == int (StripSubdetector::TIB)){
00173 TIBDetId TIBid=TIBDetId(subid);
00174 if(TIBid.string()[0]==1){LA_neg<<" Backward, ";}else{LA_neg<<" Forward, ";}
00175 if(TIBid.string()[1]==1){LA_neg<<" Int String, ";}else{LA_neg<<" Ext String, ";}
00176 LA_neg<<" Nr. String = "<<TIBid.string()[2];}
00177 if(subid.subdetId() == int (StripSubdetector::TOB)){
00178 TOBDetId TOBid=TOBDetId(subid);
00179 if(TOBid.rod()[0]==1){LA_neg<<" Backward, ";}else{LA_neg<<" Forward, ";}
00180 LA_neg<<" Nr. Rod = "<<TOBid.rod()[1];}
00181 LA_neg<<" MonoStereo = "<<subid.stereo()<<" Id = "<<id<<" => Fit_Par0 = "<<fitfunc->GetParameter(0)<<" Fit_Par1 = "<<fitfunc->GetParameter(1)<<" Fit_Par2 = "<<fitfunc->GetParameter(2)<<" Chi2/NDF = "<<fitfunc->GetChisquare()/fitfunc->GetNDF()<<std::endl<<std::endl;
00182 }
00183
00184
00185 if(fitfunc->GetParameter(0)<0&&fitfunc->GetParameter(1)>0&&fitfunc->GetParameter(2)>0){
00186 edm::LogInfo("SiStripCalibLorentzAngle")<<fitfunc->GetParameter(0);
00187 fitcounter++;
00188
00189
00190 const StripGeomDetUnit* det = dynamic_cast<const StripGeomDetUnit*>(estracker->idToDetUnit(detid));
00191 if (det==0){
00192 edm::LogError("SiStripCalibLorentzAngle") << "[SiStripCalibLorentzAngle::getNewObject] the detID " << id << " doesn't seem to belong to Tracker" << std::endl;
00193 continue;
00194 }
00195 LocalVector lbfield=(det->surface()).toLocal(magfield_->inTesla(det->surface().position()));
00196 float theBfield = lbfield.mag();
00197 theBfield = (theBfield > 0) ? theBfield : 0.00001;
00198
00199 float muH = -(fitfunc->GetParameter(0))/theBfield;
00200
00201 detid_la[id]= muH;
00202
00203 histos[1]->Fill(muH);
00204 histos[2]->Fill(fitfunc->GetParError(0)/theBfield);
00205 histos[3]->Fill(fitfunc->GetChisquare()/fitfunc->GetNDF());
00206 histos[4]->Fill(gphi,muH);
00207 histos[5]->Fill(geta,muH);
00208
00209 if(subid.subdetId() == int (StripSubdetector::TIB)){
00210 histos[6]->Fill(muH);
00211 if((Layer==1)&&(subid.stereo()==0))histos[7]->Fill(muH);
00212 if((Layer==1)&&(subid.stereo()==1))histos[8]->Fill(muH);
00213 if((Layer==2)&&(subid.stereo()==0))histos[9]->Fill(muH);
00214 if((Layer==2)&&(subid.stereo()==1))histos[10]->Fill(muH);
00215 if(Layer==3)histos[11]->Fill(muH);
00216 if(Layer==4)histos[12]->Fill(muH);
00217 }
00218
00219 if(subid.subdetId() == int (StripSubdetector::TOB)){
00220 histos[13]->Fill(muH);
00221 if((Layer==1)&&(subid.stereo()==0))histos[14]->Fill(muH);
00222 if((Layer==1)&&(subid.stereo()==1))histos[15]->Fill(muH);
00223 if((Layer==2)&&(subid.stereo()==0))histos[16]->Fill(muH);
00224 if((Layer==2)&&(subid.stereo()==1))histos[17]->Fill(muH);
00225 if(Layer==3)histos[18]->Fill(muH);
00226 if(Layer==4)histos[19]->Fill(muH);
00227 if(Layer==5)histos[20]->Fill(muH);
00228 if(Layer==6)histos[21]->Fill(muH);
00229 }
00230
00231 }else{
00232 badfit++;}
00233 }
00234 }
00235
00236 edm::LogInfo("SiStripCalibLorentzAngle")<<"### NR.OF TIB AND TOB MODULES = 7932";
00237 edm::LogInfo("SiStripCalibLorentzAngle")<<"### NR.OF HISTOS WITH ENTRIES > 100 = "<<histocounter;
00238 edm::LogInfo("SiStripCalibLorentzAngle")<<"### NO MODULE HISTOS = "<<no_mod_histo;
00239 edm::LogInfo("SiStripCalibLorentzAngle")<<"### NR.OF GOOD FIT = "<<fitcounter;
00240 edm::LogInfo("SiStripCalibLorentzAngle")<<"### NR.OF BAD FIT = "<<badfit;
00241
00242 LA_neg.close();
00243
00244 dbe_->save("LA_plots.root","LorentzAngle_Plots");
00245 }
00246
00247
00248 SiStripCalibLorentzAngle::~SiStripCalibLorentzAngle() {
00249 }
00250
00251
00252
00253
00254 SiStripLorentzAngle* SiStripCalibLorentzAngle::getNewObject(){
00255
00256 SiStripLorentzAngle* LorentzAngle = new SiStripLorentzAngle();
00257
00258 for(std::map<uint32_t, float>::iterator it = detid_la.begin(); it != detid_la.end(); it++){
00259
00260 float langle=it->second;
00261 if ( ! LorentzAngle->putLorentzAngle(it->first,langle) )
00262 edm::LogError("SiStripCalibLorentzAngle")<<"[SiStripCalibLorentzAngle::analyze] detid already exists"<<std::endl;
00263 }
00264
00265 return LorentzAngle;
00266 }
1.5.4