00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "RecoLocalMuon/CSCEfficiency/interface/CSCEfficiency.h"
00010
00011
00012 #define SQR(x) ((x)*(x))
00013
00014 #define XMIN -70.
00015 #define XMAX 70.
00016 #define YMIN -165.
00017 #define YMAX 165.
00018 #define LAYER_MIN -0.5
00019 #define LAYER_MAX 9.5
00020
00021 template <class T>
00022 inline std::string to_string (const T& t)
00023 {
00024 std::stringstream ss;
00025 ss << t;
00026 return ss.str();
00027 }
00028
00029 void Rotate(double Xinit, double Yinit, double angle, double & Xrot, double & Yrot);
00030
00031 using namespace std;
00032 using namespace edm;
00033
00034
00035
00036
00037
00038 CSCEfficiency::CSCEfficiency(const ParameterSet& pset) : theSimHitMap( pset.getParameter<edm::InputTag>("MuonCSCHits") ){
00039 const float Xmin = XMIN;
00040 const float Xmax = XMAX;
00041 const int nXbins = int(4.*(Xmax - Xmin));
00042 const float Ymin = YMIN;
00043 const float Ymax = YMAX;
00044 const int nYbins = int(2.*(Ymax - Ymin));
00045 const float Layer_min = LAYER_MIN;
00046 const float Layer_max = LAYER_MAX;
00047 const int nLayer_bins = int(Layer_max - Layer_min);
00048
00049
00050
00051 stripDigiTag_ = pset.getParameter<edm::InputTag>("stripDigiTag") ;
00052 wireDigiTag_ = pset.getParameter<edm::InputTag>("wireDigiTag") ;
00053 rootFileName = pset.getUntrackedParameter<string>("rootFileName");
00054 WorkInEndcap = pset.getUntrackedParameter<int>("WorkInEndcap");
00055 ExtrapolateFromStation = pset.getUntrackedParameter<int>("ExtrapolateFromStation");
00056 ExtrapolateToStation = pset.getUntrackedParameter<int>("ExtrapolateToStation");
00057 ExtrapolateToRing = pset.getUntrackedParameter<int>("ExtrapolateToRing");
00058 DATA = pset.getUntrackedParameter<bool>("runOnData");
00059 update = pset.getUntrackedParameter<bool>("update");
00060
00061
00062
00063 if(!DATA){
00064 mycscunpacker = pset.getUntrackedParameter<string>("mycscunpacker");
00065 }
00066
00067 nEventsAnalyzed = 0;
00068 std::string Path = "AllChambers/";
00069 std::string FullName;
00070 if(update){
00071
00072 theFile = new TFile(rootFileName.c_str(), "UPDATE");
00073 }
00074 else{
00075
00076 theFile = new TFile(rootFileName.c_str(), "RECREATE");
00077 }
00078 theFile->cd();
00079
00080 char SpecName[50];
00081
00082 sprintf(SpecName,"DataFlow");
00083 if(!update){
00084 DataFlow =
00085 new TH1F(SpecName,"Data flow;condition number;entries",30,-0.5,29.5);
00086 }
00087 else{
00088 FullName = Path + to_string(SpecName);
00089 strcpy(SpecName, FullName.c_str());
00090 DataFlow = (TH1F*)(theFile)->Get(SpecName);
00091 }
00092
00093 int Chan = 100;
00094 float minChan = 0.;
00095 float maxChan = 30.;
00096 sprintf(SpecName,"chi2_ndf");
00097 if(!update){
00098 Chi2 =
00099 new TH1F(SpecName,"Chi2/ndf;chi2/ndf;entries",Chan,minChan,maxChan);
00100 }
00101 else{
00102 FullName = Path + to_string(SpecName);
00103 strcpy(SpecName, FullName.c_str());
00104 Chi2 = (TH1F*)(theFile)->Get(SpecName);
00105 }
00106
00107 sprintf(SpecName,"chi2_ndf-ME1_a");
00108 if(!update){
00109 Chi2_ME1_a =
00110 new TH1F(SpecName,"Chi2/ndf-ME1_a;chi2/ndf;entries",Chan,minChan,maxChan);
00111 }
00112 else{
00113 FullName = Path + to_string(SpecName);
00114 strcpy(SpecName, FullName.c_str());
00115 Chi2_ME1_a = (TH1F*)(theFile)->Get(SpecName);
00116 }
00117
00118 sprintf(SpecName,"chi2_ndf-ME1_b");
00119 if(!update){
00120 Chi2_ME1_b =
00121 new TH1F(SpecName,"Chi2/ndf-ME1_b;chi2/ndf;entries",Chan,minChan,maxChan);
00122 }
00123 else{
00124 FullName = Path + to_string(SpecName);
00125 strcpy(SpecName, FullName.c_str());
00126 Chi2_ME1_b = (TH1F*)(theFile)->Get(SpecName);
00127 }
00128
00129 sprintf(SpecName,"chi2_ndf-ME1_2");
00130 if(!update){
00131 Chi2_ME1_2 =
00132 new TH1F(SpecName,"Chi2/ndf-ME1_2;chi2/ndf;entries",Chan,minChan,maxChan);
00133 }
00134 else{
00135 FullName = Path + to_string(SpecName);
00136 strcpy(SpecName, FullName.c_str());
00137 Chi2_ME1_2 = (TH1F*)(theFile)->Get(SpecName);
00138 }
00139
00140 sprintf(SpecName,"chi2_ndf-ME1_3");
00141 if(!update){
00142 Chi2_ME1_3 =
00143 new TH1F(SpecName,"Chi2/ndf-ME1_3;chi2/ndf;entries",Chan,minChan,maxChan);
00144 }
00145 else{
00146 FullName = Path + to_string(SpecName);
00147 strcpy(SpecName, FullName.c_str());
00148 Chi2_ME1_3 = (TH1F*)(theFile)->Get(SpecName);
00149 }
00150
00151 sprintf(SpecName,"chi2_ndf-ME2_1");
00152 if(!update){
00153 Chi2_ME2_1 =
00154 new TH1F(SpecName,"Chi2/ndf-ME2_1;chi2/ndf;entries",Chan,minChan,maxChan);
00155 }
00156 else{
00157 FullName = Path + to_string(SpecName);
00158 strcpy(SpecName, FullName.c_str());
00159 Chi2_ME2_1 = (TH1F*)(theFile)->Get(SpecName);
00160 }
00161
00162 sprintf(SpecName,"chi2_ndf-ME2_2");
00163 if(!update){
00164 Chi2_ME2_2 =
00165 new TH1F(SpecName,"Chi2/ndf-ME2_2;chi2/ndf;entries",Chan,minChan,maxChan);
00166 }
00167 else{
00168 FullName = Path + to_string(SpecName);
00169 strcpy(SpecName, FullName.c_str());
00170 Chi2_ME2_2 = (TH1F*)(theFile)->Get(SpecName);
00171 }
00172
00173 sprintf(SpecName,"chi2_ndf-ME3_1");
00174 if(!update){
00175 Chi2_ME3_1 =
00176 new TH1F(SpecName,"Chi2/ndf-ME3_1;chi2/ndf;entries",Chan,minChan,maxChan);
00177 }
00178 else{
00179 FullName = Path + to_string(SpecName);
00180 strcpy(SpecName, FullName.c_str());
00181 Chi2_ME3_1 = (TH1F*)(theFile)->Get(SpecName);
00182 }
00183
00184 sprintf(SpecName,"chi2_ndf-ME3_2");
00185 if(!update){
00186 Chi2_ME3_2 =
00187 new TH1F(SpecName,"Chi2/ndf-ME3_2;chi2/ndf;entries",Chan,minChan,maxChan);
00188 }
00189 else{
00190 FullName = Path + to_string(SpecName);
00191 strcpy(SpecName, FullName.c_str());
00192 Chi2_ME3_2 = (TH1F*)(theFile)->Get(SpecName);
00193 }
00194
00195 sprintf(SpecName,"chi2_ndf-ME4_1");
00196 if(!update){
00197 Chi2_ME4_1 =
00198 new TH1F(SpecName,"Chi2/ndf-ME4_1;chi2/ndf;entries",Chan,minChan,maxChan);
00199 }
00200 else{
00201 FullName = Path + to_string(SpecName);
00202 strcpy(SpecName, FullName.c_str());
00203 Chi2_ME4_1 = (TH1F*)(theFile)->Get(SpecName);
00204 }
00205
00206 sprintf(SpecName,"XY_ALCTmissing");
00207 if(!update){
00208 XY_ALCTmissing =
00209 new TH2F(SpecName,"XY - ALCT missing;cm;cm",nXbins,XMIN,XMAX,nYbins,YMIN,YMAX);
00210 }
00211 else{
00212 FullName = Path + to_string(SpecName);
00213 strcpy(SpecName, FullName.c_str());
00214 XY_ALCTmissing = (TH2F*)(theFile)->Get(SpecName);
00215 }
00216
00217 sprintf(SpecName,"dydz_Eff_ALCT");
00218 if(!update){
00219 dydz_Eff_ALCT =
00220 new TH1F(SpecName,"ALCT efficient events vs. dy/dz of the segment in ref. station;dydz;entries",30,-1.5,1.5);
00221 }
00222 else{
00223 FullName = Path + to_string(SpecName);
00224 strcpy(SpecName, FullName.c_str());
00225 dydz_Eff_ALCT = (TH1F*)(theFile)->Get(SpecName);
00226 }
00227
00228 sprintf(SpecName,"dydz_All_ALCT");
00229 if(!update){
00230 dydz_All_ALCT =
00231 new TH1F(SpecName,"ALCT events vs. dy/dz of the segment in ref. station ;dydz;entries",30,-1.5,1.5);
00232 }
00233 else{
00234 FullName = Path + to_string(SpecName);
00235 strcpy(SpecName, FullName.c_str());
00236 dydz_All_ALCT = (TH1F*)(theFile)->Get(SpecName);
00237 }
00238
00239 sprintf(SpecName,"EfficientSegments");
00240 if(!update){
00241 EfficientSegments =
00242 new TH1F(SpecName,"Efficient segments;chamber number;entries",NumCh, FirstCh-0.5, LastCh+0.5);
00243 }
00244 else{
00245 FullName = Path + to_string(SpecName);
00246 strcpy(SpecName, FullName.c_str());
00247 EfficientSegments = (TH1F*)(theFile)->Get(SpecName);
00248 }
00249
00250 sprintf(SpecName,"AllSegments");
00251 if(!update){
00252 AllSegments =
00253 new TH1F(SpecName,"All segments;chamber number;entries",NumCh, FirstCh-0.5, LastCh+0.5);
00254 }
00255 else{
00256 FullName = Path + to_string(SpecName);
00257 strcpy(SpecName, FullName.c_str());
00258 AllSegments = (TH1F*)(theFile)->Get(SpecName);
00259 }
00260
00261 sprintf(SpecName,"EfficientSegments_theta");
00262 if(!update){
00263 EfficientSegments_theta =
00264 new TH1F(SpecName,"Efficient segments in theta;theta;entries",79, 0., 3.16);
00265 }
00266 else{
00267 FullName = Path + to_string(SpecName);
00268 strcpy(SpecName, FullName.c_str());
00269 EfficientSegments_theta = (TH1F*)(theFile)->Get(SpecName);
00270 }
00271
00272 sprintf(SpecName,"AllSegments_theta");
00273 if(!update){
00274 AllSegments_theta =
00275 new TH1F(SpecName,"All segments in theta;theta;entries",79, 0., 3.16);
00276 }
00277 else{
00278 FullName = Path + to_string(SpecName);
00279 strcpy(SpecName, FullName.c_str());
00280 AllSegments_theta = (TH1F*)(theFile)->Get(SpecName);
00281 }
00282
00283
00284 sprintf(SpecName,"EfficientRechits_inSegment");
00285 if(!update){
00286 EfficientRechits_inSegment =
00287 new TH1F(SpecName,"Existing RecHit given a segment;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00288 }
00289 else{
00290 FullName = Path + to_string(SpecName)+"_AllCh";
00291 strcpy(SpecName, FullName.c_str());
00292 EfficientRechits_inSegment = (TH1F*)((theFile))->Get(SpecName);
00293 }
00294 sprintf(SpecName,"InefficientSingleHits");
00295 if(!update){
00296 InefficientSingleHits =
00297 new TH1F(SpecName,"Single RecHits not in the segment;layers (1-6);entries ",nLayer_bins,Layer_min,Layer_max);
00298 }
00299 else{
00300 FullName = Path + to_string(SpecName)+"_AllCh";
00301 strcpy(SpecName, FullName.c_str());
00302 InefficientSingleHits = (TH1F*)((theFile))->Get(SpecName);
00303 }
00304 sprintf(SpecName,"AllSingleHits");
00305 if(!update){
00306 AllSingleHits =
00307 new TH1F(SpecName,"Single RecHits given a segment; layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00308 }
00309 else{
00310 FullName = Path + to_string(SpecName)+"_AllCh";
00311 strcpy(SpecName, FullName.c_str());
00312 AllSingleHits = (TH1F*)((theFile))->Get(SpecName);
00313 }
00314
00315 sprintf(SpecName,"XvsY_InefficientRecHits");
00316 if(!update){
00317 XvsY_InefficientRecHits =
00318 new TH2F(SpecName,"Rechits if one or more layers have no any (local system);X, cm; Y, cm",
00319 nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00320 }
00321 else{
00322 FullName = Path + to_string(SpecName)+"_AllCh";
00323 strcpy(SpecName, FullName.c_str());
00324 XvsY_InefficientRecHits = (TH2F*)((theFile))->Get(SpecName);
00325 }
00326 sprintf(SpecName,"XvsY_InefficientRecHits_good");
00327 if(!update){
00328 XvsY_InefficientRecHits_good =
00329 new TH2F(SpecName,"Rechits if one or more layers have no any (local system) - sensitive area only;X, cm; Y, cm",
00330 nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00331 }
00332 else{
00333 FullName = Path + to_string(SpecName)+"_AllCh";
00334 strcpy(SpecName, FullName.c_str());
00335 XvsY_InefficientRecHits_good = (TH2F*)((theFile))->Get(SpecName);
00336 }
00337
00338 sprintf(SpecName,"XvsY_InefficientSegments");
00339 if(!update){
00340 XvsY_InefficientSegments =
00341 new TH2F(SpecName,"Segments with less than 6 hits;X, cm; Y, cm",nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00342 }
00343 else{
00344 FullName = Path + to_string(SpecName)+"_AllCh";
00345 strcpy(SpecName, FullName.c_str());
00346 XvsY_InefficientSegments = (TH2F*)((theFile))->Get(SpecName);
00347 }
00348
00349 sprintf(SpecName,"XvsY_InefficientSegments_good");
00350 if(!update){
00351 XvsY_InefficientSegments_good =
00352 new TH2F(SpecName,"Segments with less than 6 hits - sensitive area only;X, cm; Y, cm",
00353 nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00354 }
00355 else{
00356 FullName = Path + to_string(SpecName)+"_AllCh";
00357 strcpy(SpecName, FullName.c_str());
00358 XvsY_InefficientSegments_good = (TH2F*)((theFile))->Get(SpecName);
00359 }
00360
00361 sprintf(SpecName,"EfficientRechits");
00362 if(!update){
00363 EfficientRechits =
00364 new TH1F(SpecName,"Existing RecHit;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00365 }
00366 else{
00367 FullName = Path + to_string(SpecName)+"_AllCh";
00368 strcpy(SpecName, FullName.c_str());
00369 EfficientRechits = (TH1F*)((theFile))->Get(SpecName);
00370 }
00371
00372 sprintf(SpecName,"EfficientRechits_good");
00373 if(!update){
00374 EfficientRechits_good =
00375 new TH1F(SpecName,"Existing RecHit - sensitive area only;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00376 }
00377 else{
00378 FullName = Path + to_string(SpecName)+"_AllCh";
00379 strcpy(SpecName, FullName.c_str());
00380 EfficientRechits_good = (TH1F*)((theFile))->Get(SpecName);
00381 }
00382 sprintf(SpecName,"EfficientLCTs");
00383 if(!update){
00384 EfficientLCTs =
00385 new TH1F(SpecName,"Existing LCTs (1-a, 2-c, 3-corr);3 sets + normalization;entries",30,0.5,30.5);
00386 }
00387 else{
00388 FullName = Path + to_string(SpecName)+"_AllCh";
00389 strcpy(SpecName, FullName.c_str());
00390 EfficientLCTs = (TH1F*)((theFile))->Get(SpecName);
00391 }
00392
00393 sprintf(SpecName,"EfficientStrips");
00394 if(!update){
00395 EfficientStrips =
00396 new TH1F(SpecName,"Existing strip;layer (1-6); entries",nLayer_bins,Layer_min,Layer_max);
00397 }
00398 else{
00399 FullName = Path + to_string(SpecName)+"_AllCh";
00400 strcpy(SpecName, FullName.c_str());
00401 EfficientStrips = (TH1F*)((theFile))->Get(SpecName);
00402 }
00403 sprintf(SpecName,"EfficientWireGroups");
00404 if(!update){
00405 EfficientWireGroups =
00406 new TH1F(SpecName,"Existing WireGroups;layer (1-6); entries ",nLayer_bins,Layer_min,Layer_max);
00407 }
00408 else{
00409 FullName = Path + to_string(SpecName)+"_AllCh";
00410 strcpy(SpecName, FullName.c_str());
00411 EfficientWireGroups = (TH1F*)((theFile))->Get(SpecName);
00412 }
00413
00414 for(int iLayer=0; iLayer<6;iLayer++){
00415 sprintf(SpecName,"XvsY_InefficientRecHits_inSegment_L%d",iLayer);
00416 if(!update){
00417 XvsY_InefficientRecHits_inSegment.push_back
00418 (new TH2F(SpecName,"Missing RecHit/layer in a segment (local system, good region);X, cm; Y, cm",
00419 nXbins,Xmin,Xmax,nYbins,Ymin, Ymax));
00420 }
00421 else{
00422 FullName = Path + to_string(SpecName)+"_AllCh";
00423 strcpy(SpecName, FullName.c_str());
00424 XvsY_InefficientRecHits_inSegment.push_back( (TH2F*)((theFile))->Get(SpecName));
00425 }
00426
00427 sprintf(SpecName,"Y_InefficientRecHits_inSegment_L%d",iLayer);
00428 if(!update){
00429 Y_InefficientRecHits_inSegment.push_back
00430 (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
00431 nYbins,Ymin, Ymax));
00432 }
00433 else{
00434 FullName = Path + to_string(SpecName)+"_AllCh";
00435 strcpy(SpecName, FullName.c_str());
00436 Y_InefficientRecHits_inSegment.push_back( (TH1F*)((theFile))->Get(SpecName));
00437 }
00438
00439 sprintf(SpecName,"Y_AllRecHits_inSegment_L%d",iLayer);
00440 if(!update){
00441 Y_AllRecHits_inSegment.push_back
00442 (new TH1F(SpecName,"All (extrapolated from the segment) RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
00443 nYbins,Ymin, Ymax));
00444 }
00445 else{
00446 FullName = Path + to_string(SpecName)+"_AllCh";
00447 strcpy(SpecName, FullName.c_str());
00448 Y_AllRecHits_inSegment.push_back( (TH1F*)((theFile))->Get(SpecName));
00449 }
00450 }
00451
00452
00453
00454
00455 sprintf(SpecName,"Sim_Rechits");
00456 if(!update){
00457 SimRechits =
00458 new TH1F(SpecName,"Existing RecHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00459 }
00460 else{
00461 FullName = Path + to_string(SpecName)+"_AllCh";
00462 strcpy(SpecName, FullName.c_str());
00463 SimRechits = (TH1F*)((theFile))->Get(SpecName);
00464 }
00465
00466 sprintf(SpecName,"Sim_Simhits");
00467 if(!update){
00468 SimSimhits =
00469 new TH1F(SpecName,"Existing SimHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00470 }
00471 else{
00472 FullName = Path + to_string(SpecName)+"_AllCh";
00473 strcpy(SpecName, FullName.c_str());
00474 SimSimhits = (TH1F*)((theFile))->Get(SpecName);
00475 }
00476
00477 sprintf(SpecName,"Sim_Rechits_each");
00478 if(!update){
00479 SimRechits_each =
00480 new TH1F(SpecName,"Existing RecHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00481 }
00482 else{
00483 FullName = Path + to_string(SpecName)+"_AllCh";
00484 strcpy(SpecName, FullName.c_str());
00485 SimRechits_each = (TH1F*)((theFile))->Get(SpecName);
00486 }
00487
00488 sprintf(SpecName,"Sim_Simhits_each");
00489 if(!update){
00490 SimSimhits_each =
00491 new TH1F(SpecName,"Existing SimHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00492 }
00493 else{
00494 FullName = Path + to_string(SpecName)+"_AllCh";
00495 strcpy(SpecName, FullName.c_str());
00496 SimSimhits_each = (TH1F*)((theFile))->Get(SpecName);
00497 }
00498
00499
00500
00501 for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
00502 sprintf(SpecName,"Chamber_%d",iChamber);
00503 if(!update){
00504 theFile->mkdir(SpecName);
00505 }
00506 theFile->cd(SpecName);
00507 std::string Path = to_string(SpecName)+"/";
00508 sprintf(SpecName,"EfficientRechits_inSegment_Ch%d",iChamber);
00509 if(!update){
00510 ChHist[iChamber-FirstCh].EfficientRechits_inSegment =
00511 new TH1F(SpecName,"Existing RecHit given a segment;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00512 }
00513 else{
00514 FullName = Path + to_string(SpecName);
00515 strcpy(SpecName, FullName.c_str());
00516 ChHist[iChamber-FirstCh].EfficientRechits_inSegment = (TH1F*)((theFile))->Get(SpecName);
00517 }
00518
00519 sprintf(SpecName,"InefficientSingleHits_Ch%d",iChamber);
00520 if(!update){
00521 ChHist[iChamber-FirstCh].InefficientSingleHits =
00522 new TH1F(SpecName,"Single RecHits not in the segment;layers (1-6);entries ",nLayer_bins,Layer_min,Layer_max);
00523 }
00524 else{
00525 FullName = Path + to_string(SpecName);
00526 strcpy(SpecName, FullName.c_str());
00527 ChHist[iChamber-FirstCh].InefficientSingleHits = (TH1F*)((theFile))->Get(SpecName);
00528 }
00529
00530 sprintf(SpecName,"AllSingleHits_Ch%d",iChamber);
00531 if(!update){
00532 ChHist[iChamber-FirstCh].AllSingleHits =
00533 new TH1F(SpecName,"Single RecHits given a segment; layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00534 }
00535 else{
00536 FullName = Path + to_string(SpecName);
00537 strcpy(SpecName, FullName.c_str());
00538 ChHist[iChamber-FirstCh].AllSingleHits = (TH1F*)((theFile))->Get(SpecName);
00539 }
00540
00541 sprintf(SpecName,"XvsY_InefficientRecHits_Ch%d ",iChamber);
00542 if(!update){
00543 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits =
00544 new TH2F(SpecName,"Rechits if one or more layers have no any (local system);X, cm; Y, cm",
00545 nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00546 }
00547 else{
00548 FullName = Path + to_string(SpecName);
00549 strcpy(SpecName, FullName.c_str());
00550 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits = (TH2F*)((theFile))->Get(SpecName);
00551 }
00552
00553 sprintf(SpecName,"XvsY_InefficientRecHits_good_Ch%d ",iChamber);
00554 if(!update){
00555 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_good =
00556 new TH2F(SpecName,"Rechits if one or more layers have no any (local system) - sensitive area only;X, cm; Y, cm",
00557 nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00558 }
00559 else{
00560 FullName = Path + to_string(SpecName);
00561 strcpy(SpecName, FullName.c_str());
00562 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_good = (TH2F*)((theFile))->Get(SpecName);
00563 }
00564
00565 sprintf(SpecName,"XvsY_InefficientSegments_Ch%d",iChamber);
00566 if(!update){
00567 ChHist[iChamber-FirstCh].XvsY_InefficientSegments =
00568 new TH2F(SpecName,"Segments with less than 6 hits;X, cm; Y, cm",nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00569 }
00570 else{
00571 FullName = Path + to_string(SpecName);
00572 strcpy(SpecName, FullName.c_str());
00573 ChHist[iChamber-FirstCh].XvsY_InefficientSegments = (TH2F*)((theFile))->Get(SpecName);
00574 }
00575
00576 sprintf(SpecName,"XvsY_InefficientSegments_good_Ch%d",iChamber);
00577 if(!update){
00578 ChHist[iChamber-FirstCh].XvsY_InefficientSegments_good =
00579 new TH2F(SpecName,"Segments with less than 6 hits - sensitive area only;X, cm; Y, cm",
00580 nXbins,Xmin,Xmax,nYbins,Ymin,Ymax);
00581 }
00582 else{
00583 FullName = Path + to_string(SpecName);
00584 strcpy(SpecName, FullName.c_str());
00585 ChHist[iChamber-FirstCh].XvsY_InefficientSegments_good = (TH2F*)((theFile))->Get(SpecName);
00586 }
00587
00588
00589 sprintf(SpecName,"EfficientRechits_Ch%d",iChamber);
00590 if(!update){
00591 ChHist[iChamber-FirstCh].EfficientRechits =
00592 new TH1F(SpecName,"Existing RecHit;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00593 }
00594 else{
00595 FullName = Path + to_string(SpecName);
00596 strcpy(SpecName, FullName.c_str());
00597 ChHist[iChamber-FirstCh].EfficientRechits = (TH1F*)((theFile))->Get(SpecName);
00598 }
00599
00600 sprintf(SpecName,"EfficientRechits_good_Ch%d",iChamber);
00601 if(!update){
00602 ChHist[iChamber-FirstCh].EfficientRechits_good =
00603 new TH1F(SpecName,"Existing RecHit - sensitive area only;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00604 }
00605 else{
00606 FullName = Path + to_string(SpecName);
00607 strcpy(SpecName, FullName.c_str());
00608 ChHist[iChamber-FirstCh].EfficientRechits_good = (TH1F*)((theFile))->Get(SpecName);
00609 }
00610
00611 sprintf(SpecName,"EfficientLCTs_Ch%d",iChamber);
00612 if(!update){
00613 ChHist[iChamber-FirstCh].EfficientLCTs =
00614 new TH1F(SpecName,"Existing LCTs (1-a, 2-c, 3-corr);3 sets + normalization;entries",30,0.5,30.5);
00615 }
00616 else{
00617 FullName = Path + to_string(SpecName);
00618 strcpy(SpecName, FullName.c_str());
00619 ChHist[iChamber-FirstCh].EfficientLCTs = (TH1F*)((theFile))->Get(SpecName);
00620 }
00621
00622 sprintf(SpecName,"EfficientStrips_Ch%d",iChamber);
00623 if(!update){
00624 ChHist[iChamber-FirstCh].EfficientStrips =
00625 new TH1F(SpecName,"Existing strip;layer (1-6); entries",nLayer_bins,Layer_min,Layer_max);
00626 }
00627 else{
00628 FullName = Path + to_string(SpecName);
00629 strcpy(SpecName, FullName.c_str());
00630 ChHist[iChamber-FirstCh].EfficientStrips = (TH1F*)((theFile))->Get(SpecName);
00631 }
00632
00633 sprintf(SpecName,"EfficientWireGroups_Ch%d",iChamber);
00634 if(!update){
00635 ChHist[iChamber-FirstCh].EfficientWireGroups =
00636 new TH1F(SpecName,"Existing WireGroups;layer (1-6); entries ",nLayer_bins,Layer_min,Layer_max);
00637 }
00638 else{
00639 FullName = Path + to_string(SpecName);
00640 strcpy(SpecName, FullName.c_str());
00641 ChHist[iChamber-FirstCh].EfficientWireGroups = (TH1F*)((theFile))->Get(SpecName);
00642 }
00643
00644 for(int iLayer=0; iLayer<6;iLayer++){
00645 sprintf(SpecName,"XvsY_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
00646 if(!update){
00647 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_inSegment.push_back
00648 (new TH2F(SpecName,"Missing RecHit/layer in a segment (local system, good region);X, cm; Y, cm",
00649 nXbins,Xmin,Xmax,nYbins,Ymin, Ymax));
00650 }
00651 else{
00652 FullName = Path + to_string(SpecName);
00653 strcpy(SpecName, FullName.c_str());
00654 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_inSegment.push_back((TH2F*)((theFile))->Get(SpecName));
00655 }
00656
00657 sprintf(SpecName,"Y_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
00658 if(!update){
00659 ChHist[iChamber-FirstCh].Y_InefficientRecHits_inSegment.push_back
00660 (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
00661 nYbins,Ymin, Ymax));
00662 }
00663 else{
00664 FullName = Path + to_string(SpecName);
00665 strcpy(SpecName, FullName.c_str());
00666 ChHist[iChamber-FirstCh].Y_InefficientRecHits_inSegment.push_back((TH1F*)((theFile))->Get(SpecName));
00667 }
00668
00669 sprintf(SpecName,"Y_AllRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
00670 if(!update){
00671 ChHist[iChamber-FirstCh].Y_AllRecHits_inSegment.push_back
00672 (new TH1F(SpecName,"All (extrapolated from the segment) RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
00673 nYbins,Ymin, Ymax));
00674 }
00675 else{
00676 FullName = Path + to_string(SpecName);
00677 strcpy(SpecName, FullName.c_str());
00678 ChHist[iChamber-FirstCh].Y_AllRecHits_inSegment.push_back((TH1F*)((theFile))->Get(SpecName));
00679 }
00680 }
00681
00682 sprintf(SpecName,"Sim_Rechits_Ch%d",iChamber);
00683 if(!update){
00684 ChHist[iChamber-FirstCh].SimRechits =
00685 new TH1F(SpecName,"Existing RecHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00686 }
00687 else{
00688 FullName = Path + to_string(SpecName)+"_AllCh";
00689 strcpy(SpecName, FullName.c_str());
00690 ChHist[iChamber-FirstCh].SimRechits = (TH1F*)((theFile))->Get(SpecName);
00691 }
00692
00693 sprintf(SpecName,"Sim_Simhits_Ch%d",iChamber);
00694 if(!update){
00695 ChHist[iChamber-FirstCh].SimSimhits =
00696 new TH1F(SpecName,"Existing SimHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00697 }
00698 else{
00699 FullName = Path + to_string(SpecName)+"_AllCh";
00700 strcpy(SpecName, FullName.c_str());
00701 ChHist[iChamber-FirstCh].SimSimhits = (TH1F*)((theFile))->Get(SpecName);
00702 }
00703
00704 sprintf(SpecName,"Sim_Rechits_each_Ch%d",iChamber);
00705 if(!update){
00706 ChHist[iChamber-FirstCh].SimRechits_each =
00707 new TH1F(SpecName,"Existing RecHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00708 }
00709 else{
00710 FullName = Path + to_string(SpecName)+"_AllCh";
00711 strcpy(SpecName, FullName.c_str());
00712 ChHist[iChamber-FirstCh].SimRechits_each = (TH1F*)((theFile))->Get(SpecName);
00713 }
00714
00715 sprintf(SpecName,"Sim_Simhits_each_Ch%d",iChamber);
00716 if(!update){
00717 ChHist[iChamber-FirstCh].SimSimhits_each =
00718 new TH1F(SpecName,"Existing SimHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
00719 }
00720 else{
00721 FullName = Path + to_string(SpecName)+"_AllCh";
00722 strcpy(SpecName, FullName.c_str());
00723 ChHist[iChamber-FirstCh].SimSimhits_each = (TH1F*)((theFile))->Get(SpecName);
00724 }
00725
00726
00727
00728
00729
00730
00731 theFile->cd();
00732 }
00733 }
00734
00735
00736 CSCEfficiency::~CSCEfficiency(){
00737
00738 theFile->cd();
00739
00740 char SpecName[20];
00741 int Nbins;
00742 std::vector<float> bins, Efficiency, EffError;
00743 TH1F * readHisto;
00744 TH1F * writeHisto;
00745 std::vector<float> eff(2);
00746
00747
00748 const float Ymin = YMIN;
00749 const float Ymax = YMAX;
00750 const int nYbins = int(2.*(Ymax - Ymin));
00751 const float Layer_min = LAYER_MIN;
00752 const float Layer_max = LAYER_MAX-2.;
00753 const int nLayer_bins = int(Layer_max - Layer_min);
00754
00755
00756
00757 const int chMin = 1;
00758 const int chMax = 36;
00759 int chRange = chMax - chMin +1;
00760
00761 sprintf(SpecName,"eff_S");
00762 TH1F * h_effStrips =
00763 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00764 sprintf(SpecName,"eff_S2");
00765 TH1F * h_allStrips =
00766 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00767
00768 sprintf(SpecName,"eff_W");
00769 TH1F * h_effWGs =
00770 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00771 TH1F * h_allWGs =
00772 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00773
00774 sprintf(SpecName,"eff_R");
00775 TH1F * h_effRHs =
00776 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00777 TH1F * h_allRHs =
00778 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00779
00780 sprintf(SpecName,"eff_R_g");
00781 TH1F * h_effRHs_good =
00782 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00783 TH1F * h_allRHs_good =
00784 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00785
00786 sprintf(SpecName,"eff_R_in");
00787 TH1F * h_effRHs_inSeg =
00788 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00789 TH1F * h_allRHs_inSeg =
00790 new TH1F(SpecName,"Efficiency",chRange,float(chMin) -0.5,float(chMax) +0.5);
00791
00792
00793 for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
00794 sprintf(SpecName,"Chamber_%d",iChamber);
00795
00796 if(!update){
00797 if(iChamber==FirstCh){
00798 const char *current_title;
00799 const char *changed_title;
00800
00801 AllSingleHits = (TH1F*)ChHist[iChamber-FirstCh].AllSingleHits->Clone();
00802 current_title = AllSingleHits->GetName();
00803 changed_title = ChangeTitle(current_title);
00804 AllSingleHits->SetName(changed_title);
00805
00806 EfficientRechits_inSegment = (TH1F*)ChHist[iChamber-FirstCh].EfficientRechits_inSegment->Clone();
00807 current_title = EfficientRechits_inSegment->GetName();
00808 changed_title = ChangeTitle(current_title);
00809 EfficientRechits_inSegment->SetName(changed_title);
00810
00811 InefficientSingleHits = (TH1F*)ChHist[iChamber-FirstCh].InefficientSingleHits->Clone();
00812 current_title = InefficientSingleHits->GetName();
00813 changed_title = ChangeTitle(current_title);
00814 InefficientSingleHits->SetName(changed_title);
00815
00816 XvsY_InefficientRecHits = (TH2F*)ChHist[iChamber-FirstCh].XvsY_InefficientRecHits->Clone();
00817 current_title = XvsY_InefficientRecHits->GetName();
00818 changed_title = ChangeTitle(current_title);
00819 XvsY_InefficientRecHits->SetName(changed_title);
00820
00821 XvsY_InefficientRecHits_good = (TH2F*)ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_good->Clone();
00822 current_title = XvsY_InefficientRecHits_good->GetName();
00823 changed_title = ChangeTitle(current_title);
00824 XvsY_InefficientRecHits_good->SetName(changed_title);
00825
00826 XvsY_InefficientSegments = (TH2F*)ChHist[iChamber-FirstCh].XvsY_InefficientSegments->Clone();
00827 current_title = XvsY_InefficientSegments->GetName();
00828 changed_title = ChangeTitle(current_title);
00829 XvsY_InefficientSegments->SetName(changed_title);
00830
00831 XvsY_InefficientSegments_good = (TH2F*)ChHist[iChamber-FirstCh].XvsY_InefficientSegments_good->Clone();
00832 current_title = XvsY_InefficientSegments_good->GetName();
00833 changed_title = ChangeTitle(current_title);
00834 XvsY_InefficientSegments_good->SetName(changed_title);
00835
00836 EfficientRechits = (TH1F*)ChHist[iChamber-FirstCh].EfficientRechits->Clone();
00837 current_title = EfficientRechits->GetName();
00838 changed_title = ChangeTitle(current_title);
00839 EfficientRechits->SetName(changed_title);
00840
00841 EfficientRechits_good = (TH1F*)ChHist[iChamber-FirstCh].EfficientRechits_good->Clone();
00842 current_title = EfficientRechits_good->GetName();
00843 changed_title = ChangeTitle(current_title);
00844 EfficientRechits_good->SetName(changed_title);
00845
00846 EfficientLCTs = (TH1F*)ChHist[iChamber-FirstCh].EfficientLCTs->Clone();
00847 current_title = EfficientLCTs->GetName();
00848 changed_title = ChangeTitle(current_title);
00849 EfficientLCTs->SetName(changed_title);
00850
00851 EfficientStrips = (TH1F*)ChHist[iChamber-FirstCh].EfficientStrips->Clone();
00852 current_title = EfficientStrips->GetName();
00853 changed_title = ChangeTitle(current_title);
00854 EfficientStrips->SetName(changed_title);
00855
00856 EfficientWireGroups = (TH1F*)ChHist[iChamber-FirstCh].EfficientWireGroups->Clone();
00857 current_title = EfficientWireGroups->GetName();
00858 changed_title = ChangeTitle(current_title);
00859 EfficientWireGroups->SetName(changed_title);
00860 for(int iLayer=0; iLayer<6;iLayer++){
00861 XvsY_InefficientRecHits_inSegment[iLayer] =
00862 (TH2F*)ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_inSegment[iLayer]->Clone();
00863 current_title = XvsY_InefficientRecHits_inSegment[iLayer]->GetName();
00864 changed_title = ChangeTitle(current_title);
00865 XvsY_InefficientRecHits_inSegment[iLayer]->SetName(changed_title);
00866
00867 Y_InefficientRecHits_inSegment[iLayer] =
00868 (TH1F*)ChHist[iChamber-FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Clone();
00869 current_title = Y_InefficientRecHits_inSegment[iLayer]->GetName();
00870 changed_title = ChangeTitle(current_title);
00871 Y_InefficientRecHits_inSegment[iLayer]->SetName(changed_title);
00872
00873 Y_AllRecHits_inSegment[iLayer] =
00874 (TH1F*)ChHist[iChamber-FirstCh].Y_AllRecHits_inSegment[iLayer]->Clone();
00875 current_title = Y_AllRecHits_inSegment[iLayer]->GetName();
00876 changed_title = ChangeTitle(current_title);
00877 Y_AllRecHits_inSegment[iLayer]->SetName(changed_title);
00878 }
00879
00880 SimRechits=(TH1F*)ChHist[iChamber-FirstCh].SimRechits->Clone();
00881 current_title = SimRechits->GetName();
00882 changed_title = ChangeTitle(current_title);
00883 SimRechits->SetName(changed_title);
00884
00885 SimSimhits=(TH1F*)ChHist[iChamber-FirstCh].SimSimhits->Clone();
00886 current_title = SimSimhits->GetName();
00887 changed_title = ChangeTitle(current_title);
00888 SimSimhits->SetName(changed_title);
00889
00890 SimRechits_each=(TH1F*)ChHist[iChamber-FirstCh].SimRechits_each->Clone();
00891 current_title = SimRechits_each->GetName();
00892 changed_title = ChangeTitle(current_title);
00893 SimRechits_each->SetName(changed_title);
00894
00895 SimSimhits_each=(TH1F*)ChHist[iChamber-FirstCh].SimSimhits_each->Clone();
00896 current_title = SimSimhits_each->GetName();
00897 changed_title = ChangeTitle(current_title);
00898 SimSimhits_each->SetName(changed_title);
00899
00900 }
00901 else{
00902 AllSingleHits->Add(ChHist[iChamber-FirstCh].AllSingleHits);
00903 EfficientRechits_inSegment->Add(ChHist[iChamber-FirstCh].EfficientRechits_inSegment);
00904 InefficientSingleHits->Add(ChHist[iChamber-FirstCh].InefficientSingleHits);
00905 XvsY_InefficientRecHits->Add(ChHist[iChamber-FirstCh].XvsY_InefficientRecHits);
00906 XvsY_InefficientRecHits_good->Add(ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_good);
00907 XvsY_InefficientSegments->Add(ChHist[iChamber-FirstCh].XvsY_InefficientSegments);
00908 XvsY_InefficientSegments_good->Add(ChHist[iChamber-FirstCh].XvsY_InefficientSegments_good);
00909 EfficientRechits->Add(ChHist[iChamber-FirstCh].EfficientRechits);
00910 EfficientRechits_good->Add(ChHist[iChamber-FirstCh].EfficientRechits_good);
00911 EfficientLCTs->Add(ChHist[iChamber-FirstCh].EfficientLCTs);
00912 EfficientStrips->Add(ChHist[iChamber-FirstCh].EfficientStrips);
00913 EfficientWireGroups->Add(ChHist[iChamber-FirstCh].EfficientWireGroups);
00914 for(int iLayer=0; iLayer<6;iLayer++){
00915 XvsY_InefficientRecHits_inSegment[iLayer]->
00916 Add(ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_inSegment[iLayer]);
00917 Y_InefficientRecHits_inSegment[iLayer]->
00918 Add(ChHist[iChamber-FirstCh].Y_InefficientRecHits_inSegment[iLayer]);
00919 Y_AllRecHits_inSegment[iLayer]->Add(ChHist[iChamber-FirstCh].Y_AllRecHits_inSegment[iLayer]);
00920 }
00921 SimRechits->Add(ChHist[iChamber-FirstCh].SimRechits);
00922 SimSimhits->Add(ChHist[iChamber-FirstCh].SimSimhits);
00923 SimRechits_each->Add(ChHist[iChamber-FirstCh].SimRechits_each);
00924 SimSimhits_each->Add(ChHist[iChamber-FirstCh].SimSimhits_each);
00925 }
00926
00927 theFile->cd(SpecName);
00928
00929 ChHist[iChamber-FirstCh].EfficientRechits_inSegment->Write();
00930 ChHist[iChamber-FirstCh].AllSingleHits->Write();
00931 ChHist[iChamber-FirstCh].InefficientSingleHits->Write();
00932 ChHist[iChamber-FirstCh].XvsY_InefficientSegments->Write();
00933 ChHist[iChamber-FirstCh].XvsY_InefficientSegments_good->Write();
00934 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_good->Write();
00935 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits->Write();
00936 ChHist[iChamber-FirstCh].EfficientRechits->Write();
00937 ChHist[iChamber-FirstCh].EfficientRechits_good->Write();
00938 ChHist[iChamber-FirstCh].EfficientLCTs->Write();
00939 ChHist[iChamber-FirstCh].EfficientStrips->Write();
00940 ChHist[iChamber-FirstCh].EfficientWireGroups->Write();
00941 for(unsigned int iLayer = 0; iLayer< 6; iLayer++){
00942 ChHist[iChamber-FirstCh].XvsY_InefficientRecHits_inSegment[iLayer]->Write();
00943 ChHist[iChamber-FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Write();
00944 ChHist[iChamber-FirstCh].Y_AllRecHits_inSegment[iLayer]->Write();
00945 }
00946 ChHist[iChamber-FirstCh].SimRechits->Write();
00947 ChHist[iChamber-FirstCh].SimSimhits->Write();
00948 ChHist[iChamber-FirstCh].SimRechits_each->Write();
00949 ChHist[iChamber-FirstCh].SimSimhits_each->Write();
00950
00951 }
00952
00953 const int min = 2;
00954 const int max = 7;
00955 int reference = 9;
00956 float effStrips = ChHist[iChamber-FirstCh].EfficientStrips->Integral(min,max);
00957 float allStrips = 6.*ChHist[iChamber-FirstCh].EfficientStrips->GetBinContent(reference);
00958 h_effStrips->SetBinContent(iChamber,effStrips);
00959 h_allStrips->SetBinContent(iChamber,allStrips);
00960
00961 float effWGs = ChHist[iChamber-FirstCh].EfficientWireGroups->Integral(min,max);
00962 float allWGs = 6.*ChHist[iChamber-FirstCh].EfficientWireGroups->GetBinContent(reference);
00963 h_effWGs->SetBinContent(iChamber,effWGs);
00964 h_allWGs->SetBinContent(iChamber,allWGs);
00965
00966 float effRHs = ChHist[iChamber-FirstCh].EfficientRechits->Integral(min,max);
00967 float allRHs = 6.*ChHist[iChamber-FirstCh].EfficientRechits->GetBinContent(reference);
00968 h_effRHs->SetBinContent(iChamber,effRHs);
00969 h_allRHs->SetBinContent(iChamber,allRHs);
00970
00971 float effRHs_good = ChHist[iChamber-FirstCh].EfficientRechits_good->Integral(min,max);
00972 float allRHs_good = 6.*ChHist[iChamber-FirstCh].EfficientRechits_good->GetBinContent(reference);
00973 h_effRHs_good->SetBinContent(iChamber,effRHs_good);
00974 h_allRHs_good->SetBinContent(iChamber,allRHs_good);
00975
00976 reference = 10;
00977 float effRHs_inSeg = ChHist[iChamber-FirstCh].EfficientRechits_inSegment->Integral(min,max);
00978 float allRHs_inSeg = 6.*ChHist[iChamber-FirstCh].EfficientRechits_inSegment->GetBinContent(reference);
00979 h_effRHs_inSeg->SetBinContent(iChamber,effRHs_inSeg);
00980 h_allRHs_inSeg->SetBinContent(iChamber,allRHs_inSeg);
00981
00982 reference = 9;
00983
00984
00985 theFile->cd(SpecName);
00986
00987 sprintf(SpecName,"FINAL_Rechit_inSegment_Efficiency_Ch%d",iChamber);
00988 ChHist[iChamber-FirstCh].FINAL_Rechit_inSegment_Efficiency =
00989 new TH1F(SpecName,"Rechit in segment Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
00990 readHisto = ChHist[iChamber-FirstCh].EfficientRechits_inSegment;
00991 writeHisto = ChHist[iChamber-FirstCh].FINAL_Rechit_inSegment_Efficiency;
00992 histoEfficiency(readHisto, writeHisto,10);
00993 ChHist[iChamber-FirstCh].FINAL_Rechit_inSegment_Efficiency->Write("",TObject::kOverwrite);
00994
00995
00996 sprintf(SpecName,"FINAL_Attachment_Efficiency_Ch%d",iChamber);
00997 ChHist[iChamber-FirstCh].FINAL_Attachment_Efficiency =
00998 new TH1F(SpecName,"Attachment Efficiency (rechit to segment);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
00999 ChHist[iChamber-FirstCh].FINAL_Attachment_Efficiency->Sumw2();
01000 sprintf(SpecName,"efficientSegments_Ch%d",iChamber);
01001 TH1F * efficientSegments = new TH1F(SpecName,"Attachment Efficiency (rechit to segment);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01002 efficientSegments = (TH1F*)ChHist[iChamber-FirstCh].AllSingleHits->Clone();
01003 efficientSegments->Add(ChHist[iChamber-FirstCh].InefficientSingleHits,-1.);
01004 ChHist[iChamber-FirstCh].FINAL_Attachment_Efficiency->
01005 Divide(efficientSegments,
01006 ChHist[iChamber-FirstCh].AllSingleHits,
01007 1.,1.,"B");
01008 delete efficientSegments;
01009 ChHist[iChamber-FirstCh].FINAL_Attachment_Efficiency->Write();
01010
01011
01012 sprintf(SpecName,"FINAL_Rechit_Efficiency_Ch%d",iChamber);
01013 ChHist[iChamber-FirstCh].FINAL_Rechit_Efficiency =
01014 new TH1F(SpecName,"Rechit Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01015 readHisto = ChHist[iChamber-FirstCh].EfficientRechits;
01016 writeHisto = ChHist[iChamber-FirstCh].FINAL_Rechit_Efficiency;
01017 histoEfficiency(readHisto, writeHisto,9);
01018 ChHist[iChamber-FirstCh].FINAL_Rechit_Efficiency->Write();
01019
01020
01021 sprintf(SpecName,"FINAL_Rechit_Efficiency_good_Ch%d",iChamber);
01022 ChHist[iChamber-FirstCh].FINAL_Rechit_Efficiency_good =
01023 new TH1F(SpecName,"Rechit Efficiency - sensitive area only;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01024 readHisto = ChHist[iChamber-FirstCh].EfficientRechits_good;
01025 writeHisto = ChHist[iChamber-FirstCh].FINAL_Rechit_Efficiency_good;
01026 histoEfficiency(readHisto, writeHisto,9);
01027 ChHist[iChamber-FirstCh].FINAL_Rechit_Efficiency_good->Write();
01028
01029
01030 sprintf(SpecName,"FINAL_LCTs_Efficiency_Ch%d",iChamber);
01031 ChHist[iChamber-FirstCh].FINAL_LCTs_Efficiency = new TH1F(SpecName,"LCTs Efficiency;1-a, 2-c, 3-corr (3 sets);efficiency",30,0.5,30.5);
01032 Nbins = ChHist[iChamber-FirstCh].EfficientLCTs->GetSize()-2;
01033 bins.clear();
01034 bins.resize(Nbins);
01035 Efficiency.clear();
01036 Efficiency.resize(Nbins);
01037 EffError.clear();
01038 EffError.resize(Nbins);
01039 bins[Nbins-1] = ChHist[iChamber-FirstCh].EfficientLCTs->GetBinContent(Nbins);
01040 bins[Nbins-2] = ChHist[iChamber-FirstCh].EfficientLCTs->GetBinContent(Nbins-1);
01041 bins[Nbins-3] = ChHist[iChamber-FirstCh].EfficientLCTs->GetBinContent(Nbins-2);
01042 for (int i=0;i<Nbins;i++){
01043 bins[i] = ChHist[iChamber-FirstCh].EfficientLCTs->GetBinContent(i+1);
01044 float Norm = bins[Nbins-1];
01045
01046 if(i>19){
01047 Norm = bins[Nbins-3];
01048 }
01049 getEfficiency(bins[i], Norm, eff);
01050 Efficiency[i] = eff[0];
01051 EffError[i] = eff[1];
01052 ChHist[iChamber-FirstCh].FINAL_LCTs_Efficiency->SetBinContent(i+1, Efficiency[i]);
01053 ChHist[iChamber-FirstCh].FINAL_LCTs_Efficiency->SetBinError(i+1, EffError[i]);
01054 }
01055 ChHist[iChamber-FirstCh].FINAL_LCTs_Efficiency->Write();
01056
01057
01058 sprintf(SpecName,"FINAL_Strip_Efficiency_Ch%d",iChamber);
01059 ChHist[iChamber-FirstCh].FINAL_Strip_Efficiency =
01060 new TH1F(SpecName,"Strip Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01061 readHisto = ChHist[iChamber-FirstCh].EfficientStrips;
01062 writeHisto = ChHist[iChamber-FirstCh].FINAL_Strip_Efficiency;
01063 histoEfficiency(readHisto, writeHisto,9);
01064 ChHist[iChamber-FirstCh].FINAL_Strip_Efficiency->Write();
01065
01066
01067 sprintf(SpecName,"FINAL_WireGroup_Efficiency_Ch%d",iChamber);
01068 ChHist[iChamber-FirstCh].FINAL_WireGroup_Efficiency =
01069 new TH1F(SpecName,"WireGroup Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01070 readHisto = ChHist[iChamber-FirstCh].EfficientWireGroups;
01071 writeHisto = ChHist[iChamber-FirstCh].FINAL_WireGroup_Efficiency;
01072 histoEfficiency(readHisto, writeHisto,9);
01073 ChHist[iChamber-FirstCh].FINAL_WireGroup_Efficiency->Write();
01074
01075 for(int iLayer=0; iLayer<6;iLayer++){
01076 sprintf(SpecName,"FINAL_Y_RecHit_InSegment_Efficiency_Ch%d_L%d",iChamber,iLayer);
01077 ChHist[iChamber-FirstCh].FINAL_Y_RecHit_InSegment_Efficiency.push_back
01078 (new TH1F(SpecName,"RecHit/layer in a segment efficiency (local system, whole chamber);Y, cm;entries",
01079 nYbins,Ymin, Ymax));
01080 ChHist[iChamber-FirstCh].FINAL_Y_RecHit_InSegment_Efficiency.back()->Sumw2();
01081 sprintf(SpecName,"efficientRecHits_Ch%d_L%d",iChamber,iLayer);
01082 TH1F *efficientRecHits_Y = new TH1F(SpecName,"RecHit/layer in a segment efficiency (local system, whole chamber);Y, cm;entries",
01083 nYbins,Ymin, Ymax);
01084 efficientRecHits_Y = (TH1F*)ChHist[iChamber-FirstCh].Y_AllRecHits_inSegment.back()->Clone();
01085 efficientRecHits_Y->Add(ChHist[iChamber-FirstCh].Y_InefficientRecHits_inSegment.back(),-1.);
01086 ChHist[iChamber-FirstCh].FINAL_Y_RecHit_InSegment_Efficiency.back()->
01087 Divide(efficientRecHits_Y,
01088 ChHist[iChamber-FirstCh].Y_AllRecHits_inSegment[iLayer],
01089 1.,1.,"B");
01090 delete efficientRecHits_Y;
01091 ChHist[iChamber-FirstCh].FINAL_Y_RecHit_InSegment_Efficiency.back()->Write();
01092 }
01093
01094 sprintf(SpecName,"FINAL_SimRechit_Efficiency_Ch%d",iChamber);
01095 ChHist[iChamber-FirstCh].FINAL_SimRechit_Efficiency =
01096 new TH1F(SpecName,"Rechit Efficiency (Nrechit layers / Nsimhit layers);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01097 ChHist[iChamber-FirstCh].FINAL_SimRechit_Efficiency->Sumw2();
01098 sprintf(SpecName,"SimRechits_Ch%d",iChamber);
01099 TH1F * SimRechits = new TH1F(SpecName,"Rechit Efficiency (Nrechit layers / Nsimhit layers);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01100 SimRechits = (TH1F*)ChHist[iChamber-FirstCh].SimRechits->Clone();
01101 ChHist[iChamber-FirstCh].FINAL_SimRechit_Efficiency->
01102 Divide(SimRechits,
01103 ChHist[iChamber-FirstCh].SimSimhits,
01104 1.,1.,"B");
01105 delete SimRechits;
01106 ChHist[iChamber-FirstCh].FINAL_SimRechit_Efficiency->Write();
01107
01108 sprintf(SpecName,"FINAL_SimRechit_each_Efficiency_Ch%d",iChamber);
01109 ChHist[iChamber-FirstCh].FINAL_SimRechit_each_Efficiency =
01110 new TH1F(SpecName,"Rechit Efficiency (Nrechits/Nsimhits);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01111 ChHist[iChamber-FirstCh].FINAL_SimRechit_each_Efficiency->Sumw2();
01112 sprintf(SpecName,"SimRechits_each_Ch%d",iChamber);
01113 TH1F * SimRechits_each = new TH1F(SpecName,"Rechit Efficiency (Nrechits/ Nsimhits);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01114 SimRechits_each = (TH1F*)ChHist[iChamber-FirstCh].SimRechits_each->Clone();
01115 ChHist[iChamber-FirstCh].FINAL_SimRechit_each_Efficiency->
01116 Divide(SimRechits_each,
01117 ChHist[iChamber-FirstCh].SimSimhits_each,
01118 1.,1.,"B");
01119 delete SimRechits_each;
01120 ChHist[iChamber-FirstCh].FINAL_SimRechit_each_Efficiency->Write();
01121
01122 theFile->cd();
01123
01124 }
01125 sprintf(SpecName,"AllChambers");
01126 if(!update){
01127 theFile->mkdir(SpecName);
01128 theFile->cd(SpecName);
01129 DataFlow->Write();
01130
01131 Chi2->Write();
01132 Chi2_ME1_a->Write();
01133 Chi2_ME1_b->Write();
01134 Chi2_ME1_2->Write();
01135 Chi2_ME1_3->Write();
01136 Chi2_ME2_1->Write();
01137 Chi2_ME2_2->Write();
01138 Chi2_ME3_1->Write();
01139 Chi2_ME3_2->Write();
01140 Chi2_ME4_1->Write();
01141
01142 XY_ALCTmissing->Write();
01143 EfficientSegments->Write();
01144 AllSegments->Write();
01145 EfficientSegments_theta->Write();
01146 AllSegments_theta->Write();
01147
01148 EfficientRechits_inSegment->Write();
01149 AllSingleHits->Write();
01150 InefficientSingleHits->Write();
01151 XvsY_InefficientRecHits->Write();
01152 XvsY_InefficientRecHits_good->Write();
01153 XvsY_InefficientSegments->Write();
01154 XvsY_InefficientSegments_good->Write();
01155 EfficientRechits->Write();
01156 EfficientRechits_good->Write();
01157 EfficientLCTs->Write();
01158 EfficientStrips->Write();
01159 EfficientWireGroups->Write();
01160 for(unsigned int iLayer = 0; iLayer< 6; iLayer++){
01161 XvsY_InefficientRecHits_inSegment[iLayer]->Write();
01162 Y_InefficientRecHits_inSegment[iLayer]->Write();
01163 Y_AllRecHits_inSegment[iLayer]->Write();
01164 }
01165 SimRechits->Write();
01166 SimSimhits->Write();
01167 SimRechits_each->Write();
01168 SimSimhits_each->Write();
01169
01170 }
01171 theFile->cd(SpecName);
01172
01173 TGraphAsymmErrors* g_effStrips = new TGraphAsymmErrors(chRange);
01174 g_effStrips->BayesDivide(h_effStrips,h_allStrips);
01175 g_effStrips->SetTitle("Strip efficiency per chamber ;Chamber # ; efficiency");
01176 g_effStrips->SetName("FINAL_Strip_Efficiency_perChamber");
01177
01178 g_effStrips->Write();
01179 delete h_effStrips;
01180 delete h_allStrips;
01181
01182 TGraphAsymmErrors* g_effWGs = new TGraphAsymmErrors(chRange);
01183 g_effWGs->BayesDivide(h_effWGs,h_allWGs);
01184 g_effWGs->SetTitle("Wire group efficiency per chamber ;Chamber # ; efficiency");
01185 g_effWGs->SetName("FINAL_WireGroup_Efficiency_perChamber");
01186 g_effWGs->Write();
01187 delete h_effWGs;
01188 delete h_allWGs;
01189
01190 TGraphAsymmErrors* g_effRHs = new TGraphAsymmErrors(chRange);
01191 g_effRHs->BayesDivide(h_effRHs,h_allRHs);
01192 g_effRHs->SetTitle("Rechit efficiency per chamber ;Chamber # ; efficiency");
01193 g_effRHs->SetName("FINAL_Rechit_Efficiency_perChamber");
01194 g_effRHs->Write();
01195
01196 delete h_effRHs;
01197 delete h_allRHs;
01198
01199 TGraphAsymmErrors* g_effRHs_good = new TGraphAsymmErrors(chRange);
01200 g_effRHs_good->BayesDivide(h_effRHs_good,h_allRHs_good);
01201 g_effRHs_good->SetTitle("Rechit (sensitive region) efficiency per chamber ;Chamber # ; efficiency");
01202 g_effRHs_good->SetName("FINAL_Rechit_good_Efficiency_perChamber");
01203 g_effRHs_good->Write();
01204
01205 delete h_effRHs_good;
01206 delete h_allRHs_good;
01207
01208 TGraphAsymmErrors* g_effRHs_inSeg = new TGraphAsymmErrors(chRange);
01209 g_effRHs_inSeg->BayesDivide(h_effRHs_inSeg,h_allRHs_inSeg);
01210 g_effRHs_inSeg->SetTitle("Rechit (in Segment) efficiency per chamber ;Chamber # ; efficiency");
01211 g_effRHs_inSeg->SetName("FINAL_Rechit_inSegment_Efficiency_perChamber");
01212 g_effRHs_inSeg->Write();
01213
01214 delete h_effRHs_inSeg;
01215 delete h_allRHs_inSeg;
01216
01217
01218 sprintf(SpecName,"FINAL_dydz_Efficiency_ALCT");
01219 FINAL_dydz_Efficiency_ALCT=
01220 new TH1F(SpecName,"ALCT efficiency vs dy/dz of the segment in ref. system;dydz;efficiency", 30, -1.5, 1.5);
01221 FINAL_dydz_Efficiency_ALCT->Sumw2();
01222
01223 FINAL_dydz_Efficiency_ALCT->Divide(dydz_Eff_ALCT, dydz_All_ALCT, 1.,1.,"B");
01224
01225
01226 sprintf(SpecName,"FINAL_dydz_Efficiency_ALCT");
01227 FINAL_dydz_Efficiency_ALCT=
01228 new TH1F(SpecName,"ALCT efficiency vs dy/dz of the segment in ref. system;dydz;efficiency", 30, -1.5, 1.5);
01229 FINAL_dydz_Efficiency_ALCT->Sumw2();
01230
01231 FINAL_dydz_Efficiency_ALCT->Divide(dydz_Eff_ALCT, dydz_All_ALCT, 1.,1.,"B");
01232 FINAL_dydz_Efficiency_ALCT ->Write();
01233 dydz_Eff_ALCT->Write();
01234 dydz_All_ALCT->Write();
01235
01236
01237 sprintf(SpecName,"FINAL_Segment_Efficiency");
01238 FINAL_Segment_Efficiency =
01239 new TH1F(SpecName,"Segment Efficiency;chamber number;efficiency", NumCh, FirstCh-0.5, LastCh+0.5);
01240 FINAL_Segment_Efficiency->Sumw2();
01241 FINAL_Segment_Efficiency->
01242 Divide(EfficientSegments,
01243 AllSegments,
01244 1.,1.,"B");
01245 FINAL_Segment_Efficiency->Write();
01246
01247 sprintf(SpecName,"FINAL_Segment_Efficiency_theta");
01248 FINAL_Segment_Efficiency_theta =
01249 new TH1F(SpecName,"Segment Efficiency in theta;theta;efficiency", 79, 0., 3.16);
01250 FINAL_Segment_Efficiency_theta->Sumw2();
01251 FINAL_Segment_Efficiency_theta->
01252 Divide(EfficientSegments_theta,
01253 AllSegments_theta,
01254 1.,1.,"B");
01255 FINAL_Segment_Efficiency_theta->Write();
01256
01257 sprintf(SpecName,"FINAL_Rechit_inSegment_Efficiency");
01258 FINAL_Rechit_inSegment_Efficiency =
01259 new TH1F(SpecName,"Rechit in segment Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01260 readHisto = EfficientRechits_inSegment;
01261 writeHisto = FINAL_Rechit_inSegment_Efficiency;
01262 histoEfficiency(readHisto, writeHisto,10);
01263 FINAL_Rechit_inSegment_Efficiency->Write();
01264
01265
01266 sprintf(SpecName,"FINAL_Attachment_Efficiency");
01267 FINAL_Attachment_Efficiency =
01268 new TH1F(SpecName,"Attachment Efficiency (rechit to segment);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01269 FINAL_Attachment_Efficiency->Sumw2();
01270 sprintf(SpecName,"efficientSegments");
01271 TH1F * efficientSegments = new TH1F(SpecName,"Attachment Efficiency (rechit to segment);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01272 efficientSegments = (TH1F*)AllSingleHits->Clone();
01273 efficientSegments->Add(InefficientSingleHits,-1.);
01274 FINAL_Attachment_Efficiency->
01275 Divide(efficientSegments,
01276 AllSingleHits,
01277 1.,1.,"B");
01278 delete efficientSegments;
01279 FINAL_Attachment_Efficiency->Write();
01280
01281
01282 sprintf(SpecName,"FINAL_Rechit_Efficiency");
01283 FINAL_Rechit_Efficiency =
01284 new TH1F(SpecName,"Rechit Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01285 readHisto = EfficientRechits;
01286 writeHisto = FINAL_Rechit_Efficiency;
01287 histoEfficiency(readHisto, writeHisto,9);
01288 FINAL_Rechit_Efficiency->Write();
01289
01290
01291 sprintf(SpecName,"FINAL_Rechit_Efficiency_good");
01292 FINAL_Rechit_Efficiency_good =
01293 new TH1F(SpecName,"Rechit Efficiency - sensitive area only;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01294 readHisto = EfficientRechits_good;
01295 writeHisto = FINAL_Rechit_Efficiency_good;
01296 histoEfficiency(readHisto, writeHisto,9);
01297 FINAL_Rechit_Efficiency_good->Write();
01298
01299
01300 sprintf(SpecName,"FINAL_LCTs_Efficiency");
01301 FINAL_LCTs_Efficiency = new TH1F(SpecName,"LCTs Efficiency;1-a, 2-c, 3-corr (3 sets);efficiency",30,0.5,30.5);
01302 Nbins = EfficientLCTs->GetSize()-2;
01303 bins.clear();
01304 bins.resize(Nbins);
01305 Efficiency.clear();
01306 Efficiency.resize(Nbins);
01307 EffError.clear();
01308 EffError.resize(Nbins);
01309 bins[Nbins-1] = EfficientLCTs->GetBinContent(Nbins);
01310 bins[Nbins-2] = EfficientLCTs->GetBinContent(Nbins-1);
01311 bins[Nbins-3] = EfficientLCTs->GetBinContent(Nbins-2);
01312 for (int i=0;i<Nbins;i++){
01313 bins[i] = EfficientLCTs->GetBinContent(i+1);
01314 float Norm = bins[Nbins-1];
01315
01316 if(i>19){
01317 Norm = bins[Nbins-3];
01318 }
01319 getEfficiency(bins[i], Norm, eff);
01320 Efficiency[i] = eff[0];
01321 EffError[i] = eff[1];
01322 FINAL_LCTs_Efficiency->SetBinContent(i+1, Efficiency[i]);
01323 FINAL_LCTs_Efficiency->SetBinError(i+1, EffError[i]);
01324 }
01325 FINAL_LCTs_Efficiency->Write();
01326
01327
01328 sprintf(SpecName,"FINAL_Strip_Efficiency");
01329 FINAL_Strip_Efficiency =
01330 new TH1F(SpecName,"Strip Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01331 readHisto = EfficientStrips;
01332 writeHisto = FINAL_Strip_Efficiency;
01333 histoEfficiency(readHisto, writeHisto,9);
01334 FINAL_Strip_Efficiency->Write();
01335
01336
01337 sprintf(SpecName,"FINAL_WireGroup_Efficiency");
01338 FINAL_WireGroup_Efficiency =
01339 new TH1F(SpecName,"WireGroup Efficiency;layer (1-6);efficiency",nLayer_bins,Layer_min,Layer_max);
01340 readHisto = EfficientWireGroups;
01341 writeHisto = FINAL_WireGroup_Efficiency;
01342 histoEfficiency(readHisto, writeHisto,9);
01343 FINAL_WireGroup_Efficiency->Write();
01344
01345 for(int iLayer=0; iLayer<6;iLayer++){
01346 sprintf(SpecName,"FINAL_Y_RecHit_InSegment_Efficiency_L%d",iLayer);
01347 FINAL_Y_RecHit_InSegment_Efficiency.push_back
01348 (new TH1F(SpecName,"RecHit/layer in a segment efficiency (local system);Y, cm;entries",
01349 nYbins,Ymin, Ymax));
01350 FINAL_Y_RecHit_InSegment_Efficiency[iLayer]->Sumw2();
01351 sprintf(SpecName,"efficientRecHits_L%d",iLayer);
01352 TH1F *efficientRecHits_Y = new TH1F(SpecName,"RecHit/layer in a segment efficiency (local system, whole chamber);Y, cm;entries",
01353 nYbins,Ymin, Ymax);
01354 efficientRecHits_Y = (TH1F*)Y_AllRecHits_inSegment[iLayer]->Clone();
01355 efficientRecHits_Y->Add(Y_InefficientRecHits_inSegment[iLayer],-1.);
01356 FINAL_Y_RecHit_InSegment_Efficiency[iLayer]->
01357 Divide(efficientRecHits_Y,
01358 Y_AllRecHits_inSegment[iLayer],
01359 1.,1.,"B");
01360 delete efficientRecHits_Y;
01361 FINAL_Y_RecHit_InSegment_Efficiency[iLayer]->Write();
01362 }
01363
01364 sprintf(SpecName,"FINAL_SimRechit_Efficiency");
01365 FINAL_SimRechit_Efficiency =
01366 new TH1F(SpecName,"Rechit Efficiency (Nrechit layers / Nsimhit layers);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01367 FINAL_SimRechit_Efficiency->Sumw2();
01368 sprintf(SpecName,"SimRechits_tmp");
01369 TH1F * SimRechits_tmp = new TH1F(SpecName,"Rechit Efficiency (Nrechit layers / Nsimhit layers);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01370 SimRechits_tmp = (TH1F*)SimRechits->Clone();
01371 FINAL_SimRechit_Efficiency->
01372 Divide(SimRechits_tmp,
01373 SimSimhits,
01374 1.,1.,"B");
01375 delete SimRechits_tmp;
01376 FINAL_SimRechit_Efficiency->Write();
01377
01378 sprintf(SpecName,"FINAL_SimRechit_each_Efficiency");
01379 FINAL_SimRechit_each_Efficiency =
01380 new TH1F(SpecName,"Rechit Efficiency (Nrechits/ Nsimhits);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01381 FINAL_SimRechit_each_Efficiency->Sumw2();
01382 sprintf(SpecName,"SimRechits_each_tmp");
01383 TH1F * SimRechits_each_tmp = new TH1F(SpecName,"Rechit Efficiency (Nrechits/ Nsimhits);layer (1-6);efficiency",nLayer_bins+2,Layer_min,Layer_max+2.);
01384 SimRechits_each_tmp = (TH1F*)SimRechits_each->Clone();
01385 FINAL_SimRechit_each_Efficiency->
01386 Divide(SimRechits_each_tmp,
01387 SimSimhits_each,
01388 1.,1.,"B");
01389 delete SimRechits_each_tmp;
01390 FINAL_SimRechit_each_Efficiency->Write();
01391
01392
01393
01394 theFile->Close();
01395 }
01396
01397
01398 void CSCEfficiency::analyze(const Event & event, const EventSetup& eventSetup){
01399 DataFlow->Fill(0.);
01400
01401 nEventsAnalyzed++;
01402
01403 printalot = (nEventsAnalyzed < 100);
01404 int iRun = event.id().run();
01405 int iEvent = event.id().event();
01406 if(0==fmod(double (nEventsAnalyzed) ,double(100) )){
01407 printf("\n==enter==CSCEfficiency===== run %i\tevent %i\tn Analyzed %i\n",iRun,iEvent,nEventsAnalyzed);
01408 }
01409
01410
01411
01412 if (printalot) printf("\tget handles for digi collections\n");
01413 edm::Handle<CSCWireDigiCollection> wires;
01414 edm::Handle<CSCStripDigiCollection> strips;
01415
01416
01417
01418
01419 if (printalot) printf("\tpass handles\n");
01420 if(DATA){
01421
01422
01423 event.getByLabel( stripDigiTag_, strips);
01424 event.getByLabel( wireDigiTag_, wires);
01425 }
01426 else{
01427 theSimHitMap.fill(event);
01428
01429
01430 event.getByLabel( stripDigiTag_, strips);
01431 event.getByLabel( wireDigiTag_, wires);
01432 }
01433
01434
01435 if (printalot) printf("\tget the CSC geometry.\n");
01436 ESHandle<CSCGeometry> cscGeom;
01437 eventSetup.get<MuonGeometryRecord>().get(cscGeom);
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448 for(int iE=0;iE<2;iE++){
01449 for(int iS=0;iS<4;iS++){
01450 for(int iR=0;iR<4;iR++){
01451 for(int iC=0;iC<NumCh;iC++){
01452 for(int iL=0;iL<6;iL++){
01453 AllWG[iE][iS][iR][iC][iL].clear();
01454 AllStrips[iE][iS][iR][iC][iL].clear();
01455 }
01456 }
01457 }
01458 }
01459 }
01460
01461 for (CSCWireDigiCollection::DigiRangeIterator j=wires->begin(); j!=wires->end(); j++) {
01462 CSCDetId id = (CSCDetId)(*j).first;
01463 const CSCLayer *layer_p = cscGeom->layer (id);
01464 const CSCLayerGeometry *layerGeom = layer_p->geometry ();
01465 const std::vector<float> LayerBounds = layerGeom->parameters ();
01466 std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
01467 std::vector<CSCWireDigi>::const_iterator last = (*j).second.second;
01468
01469 for( ; digiItr != last; ++digiItr) {
01470 std::pair < int, float > WG_pos(digiItr->getWireGroup(), layerGeom->yOfWireGroup(digiItr->getWireGroup()));
01471 std::pair <std::pair < int, float >, int > LayerSignal(WG_pos, digiItr->getTimeBin());
01472
01473
01474 AllWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
01475 [id.layer()-1].push_back(LayerSignal);
01476 }
01477 }
01478
01479
01480 for (CSCStripDigiCollection::DigiRangeIterator j=strips->begin(); j!=strips->end(); j++) {
01481 CSCDetId id = (CSCDetId)(*j).first;
01482 const CSCLayer *layer_p = cscGeom->layer (id);
01483 const CSCLayerGeometry *layerGeom = layer_p->geometry ();
01484 const std::vector<float> LayerBounds = layerGeom->parameters ();
01485 int largestADCValue = -1;
01486 int largestStrip = -1;
01487 std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
01488 std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;
01489 for( ; digiItr != last; ++digiItr) {
01490 int maxADC=largestADCValue;
01491 int myStrip = digiItr->getStrip();
01492 std::vector<int> myADCVals = digiItr->getADCCounts();
01493 bool thisStripFired = false;
01494 float thisPedestal = 0.5*(float)(myADCVals[0]+myADCVals[1]);
01495 float threshold = 13.3 ;
01496 float diff = 0.;
01497 float peakADC = -1000.;
01498 int peakTime = -1;
01499 for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
01500 diff = (float)myADCVals[iCount]-thisPedestal;
01501 if (diff > threshold) {
01502 thisStripFired = true;
01503 if (myADCVals[iCount] > largestADCValue) {
01504 largestADCValue = myADCVals[iCount];
01505 largestStrip = myStrip;
01506 }
01507 }
01508 if (diff > threshold && diff > peakADC) {
01509 peakADC = diff;
01510 peakTime = iCount;
01511 }
01512 }
01513 if(largestADCValue>maxADC){
01514 maxADC = largestADCValue;
01515 std::pair <int, float> LayerSignal (myStrip, peakADC);
01516
01517
01518
01519 AllStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].clear();
01520 AllStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].push_back(LayerSignal);
01521 }
01522 }
01523 }
01524
01525
01526
01527
01528
01529
01530
01531
01532 if (printalot) printf("\tGet the recHits collection.\t");
01533 Handle<CSCRecHit2DCollection> recHits;
01534 event.getByLabel("csc2DRecHits",recHits);
01535 int nRecHits = recHits->size();
01536 if (printalot) printf(" The size is %i\n",nRecHits);
01537
01538 SetOfRecHits AllRecHits[2][4][4][ NumCh];
01539 std::vector<bool> InitVectBool6(6);
01540 std::vector<double> InitVectD6(6);
01541 map<int, std::vector <bool> > MyRecHits;
01542 map<int, std::vector <double> > MyRecHitsPosX;
01543 map<int, std::vector <double> > MyRecHitsPosY;
01544 map<int, std::vector <double> > MyRecHitsPosZ;
01545 for(int iSt=0;iSt<NumCh;iSt++){
01546 MyRecHits[iSt] = InitVectBool6;
01547 MyRecHitsPosX[iSt] = MyRecHitsPosY[iSt] = MyRecHitsPosZ[iSt] = InitVectD6;
01548 }
01549
01550
01551
01552 if (printalot) printf("\t...start loop over rechits...\n");
01553
01554 CSCRecHit2DCollection::const_iterator recIt;
01555
01556 for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
01557
01558 CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
01559 int kEndcap = idrec.endcap();
01560 int kRing = idrec.ring();
01561 int kStation = idrec.station();
01562 int kChamber = idrec.chamber();
01563 int kLayer = idrec.layer();
01564 if (printalot) printf("\t\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",kEndcap,kStation,kRing,kChamber,kLayer);
01565
01566 LocalPoint rhitlocal = (*recIt).localPosition();
01567 double xreco = rhitlocal.x();
01568 double yreco = rhitlocal.y();
01569 double zreco = rhitlocal.z();
01570 LocalError rerrlocal = (*recIt).localPositionError();
01571 double xxerr = rerrlocal.xx();
01572 double yyerr = rerrlocal.yy();
01573 double xyerr = rerrlocal.xy();
01574
01575
01576 const CSCLayer* csclayer = cscGeom->layer( idrec );
01577
01578
01579 GlobalPoint rhitglobal= csclayer->toGlobal(rhitlocal);
01580
01581 double grecx = rhitglobal.x();
01582 double grecy = rhitglobal.y();
01583 double grecz = rhitglobal.z();
01584
01585
01586 if(WorkInEndcap==kEndcap &&
01587 ExtrapolateToStation==kStation &&
01588 ExtrapolateToRing==kRing){
01589 if(kChamber>=FirstCh && kChamber<=LastCh){
01590 MyRecHits[kChamber-FirstCh][kLayer-1] = true;
01591 MyRecHitsPosX[kChamber-FirstCh][kLayer-1] = rhitglobal.x();
01592 MyRecHitsPosY[kChamber-FirstCh][kLayer-1] = rhitglobal.y();
01593 MyRecHitsPosZ[kChamber-FirstCh][kLayer-1] = rhitglobal.z();
01594 }
01595 }
01596
01597
01598 ChamberRecHits *ThisChamber = &AllRecHits[kEndcap-1][kStation-1][kRing-1][kChamber-FirstCh].sChamber;
01599 std::vector <double> *vec_p = &ThisChamber->RecHitsPosX[kLayer-1];
01600 vec_p->push_back(grecx);
01601 vec_p = &ThisChamber->RecHitsPosY[kLayer-1];
01602 vec_p->push_back(grecy);
01603 vec_p = &ThisChamber->RecHitsPosZ[kLayer-1];
01604 vec_p->push_back(grecz);
01605 vec_p = &ThisChamber->RecHitsPosXlocal[kLayer-1];
01606 vec_p->push_back(xreco);
01607 vec_p = &ThisChamber->RecHitsPosYlocal[kLayer-1];
01608 vec_p->push_back(yreco);
01609
01610
01611 if (printalot) printf("\t\t\tx,y,z: %f, %f, %f\texx,eey,exy: %f, %f, %f\tglobal x,y,z: %f, %f, %f \n",xreco,yreco,zreco,xxerr,yyerr,xyerr,grecx,grecy,grecz);
01612 }
01613
01614 for(int ii=0;ii<2;ii++){
01615 for(int jj=0;jj<4;jj++){
01616 for(int kk=0;kk<4;kk++){
01617 for(int ll=0;ll<LastCh-FirstCh+1;ll++){
01618 for(int mm = 0;mm<6;mm++){
01619 int new_size =
01620 AllRecHits[ii][jj][kk][ll].sChamber.RecHitsPosX[mm].size();
01621
01622 AllRecHits[ii][jj][kk][ll].sChamber.NRecHits[mm]=new_size;
01623
01624 if((WorkInEndcap-1) == ii &&
01625 (ExtrapolateToStation-1) == jj &&
01626 (ExtrapolateToRing-1)== kk){
01627
01628
01629 if( 1!=new_size){
01630 MyRecHits[ll][mm] = false;
01631 }
01632 }
01633 }
01634 }
01635 }
01636 }
01637 }
01638
01639
01640
01641 SetOfSimHits AllSimHits[2][4][4][ NumCh];
01642 if(!DATA){
01643
01644
01645 std::vector<int> detsWithHits = theSimHitMap.detsWithHits();
01646 int rawdetId = 0;
01647 if (printalot) printf("\t Loop over simhits... \n");
01648 for(unsigned int iDet=0;iDet<detsWithHits.size();iDet++){
01649 rawdetId = detsWithHits[iDet];
01650 edm::PSimHitContainer simHits = theSimHitMap.hits(rawdetId);
01651 CSCDetId simId = (CSCDetId)(*(simHits.begin())).detUnitId();
01652 int kEndcap = simId.endcap();
01653 int kRing = simId.ring();
01654 int kStation = simId.station();
01655 int kChamber = simId.chamber();
01656 int kLayer = simId.layer();
01657
01658
01659 AllSimHits[kEndcap-1][kStation-1][kRing-1][kChamber-FirstCh].sChamber.NSimHits[kLayer-1] = simHits.size();
01660 for(unsigned int simhit = 0; simhit <simHits.size(); ++simhit) {
01661 AllSimHits[kEndcap-1][kStation-1][kRing-1][kChamber-FirstCh].sChamber.SimHitsPosXlocal[kLayer-1].push_back(simHits[simhit].localPosition().x());
01662 AllSimHits[kEndcap-1][kStation-1][kRing-1][kChamber-FirstCh].sChamber.SimHitsPosYlocal[kLayer-1].push_back(simHits[simhit].localPosition().y());
01663 AllSimHits[kEndcap-1][kStation-1][kRing-1][kChamber-FirstCh].sChamber.SimHitsEnergy[kLayer-1].push_back(simHits[simhit].energyLoss());
01664 AllSimHits[kEndcap-1][kStation-1][kRing-1][kChamber-FirstCh].sChamber.SimHitsPID[kLayer-1].push_back(simHits[simhit].particleType());
01665 }
01666 }
01667 all_SimHits = &AllSimHits;
01668 }
01669
01670
01671
01672
01673
01674
01675
01676
01677 if (printalot) printf("\tGet CSC segment collection...\n");
01678 Handle<CSCSegmentCollection> cscSegments;
01679 event.getByLabel("cscSegments", cscSegments);
01680 int nSegments = cscSegments->size();
01681 if (printalot) printf(" The size is %i\n",nSegments);
01682
01683
01684 int iSegment = 0;
01685
01686
01687 int Couple = 2;
01688 std::vector<int> InitVect6(6);
01689 std::vector<int> InitVectCh(NumCh);
01690 std::vector<double> InitVect3(3);
01691
01692 map<int, std::vector<int> > ChambersWithSegments;
01693
01694 map<int, std::vector<double> > PosLocalCouple;
01695
01696 map<int, std::vector<double> > DirLocalCouple;
01697
01698 map<int, std::vector<double> > PosCouple;
01699
01700 map<int, std::vector<double> > DirCouple;
01701
01702 map<int, std::vector<double> > ChamberBoundsCouple;
01703
01704 for(int iCop=0;iCop<2;iCop++){
01705 PosLocalCouple[iCop] =
01706 DirLocalCouple[iCop] =
01707 PosCouple[iCop] =
01708 DirCouple[iCop] =
01709 ChamberBoundsCouple[iCop] = InitVect3;
01710 ChambersWithSegments[iCop] = InitVectCh;
01711 }
01712 std::vector<double> Chi2Couple(Couple);
01713 std::vector<double> NDFCouple(Couple);
01714 std::vector<double> NhitsCouple(Couple);
01715 std::vector<double> XchamberCouple(Couple);
01716 std::vector<double> YchamberCouple(Couple);
01717 std::vector<double> RotPhiCouple(Couple);
01718 std::vector<int> goodSegment(Couple);
01719 std::vector<int> NChamberCouple(Couple);
01720 std::vector<int> LayersInSecondCouple(6);
01721
01722 std::vector<int> SegmentInChamber;
01723 std::vector<int> SegmentInRing;
01724 std::vector<int> SegmentInStation;
01725 double rotationPhi = 0.;
01726 std::vector <double> DirGlobal_ThirdSegment;
01727
01728
01729 int thisEndcap = -99;
01730 int thisRing = -99;
01731 int thisStation = - 99;
01732 int thisChamber = -99;
01733 for(CSCSegmentCollection::const_iterator it=cscSegments->begin(); it != cscSegments->end(); it++) {
01734 CSCDetId id = (CSCDetId)(*it).cscDetId();
01735 if(thisEndcap==id.endcap() && thisRing==id.ring() &&
01736 thisStation == id.station() && thisChamber == id.chamber()){
01737 AllRecHits[thisEndcap-1][thisStation-1][thisRing-1][thisChamber-FirstCh].sChamber.nSegments++;
01738
01739 }
01740 else{
01741 AllRecHits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh].sChamber.nSegments = 1;
01742 thisEndcap=id.endcap();
01743 thisRing=id.ring();
01744 thisStation = id.station();
01745 thisChamber = id.chamber();
01746 }
01747 }
01748
01749 all_RecHits = &AllRecHits;
01750
01751
01752 for(CSCSegmentCollection::const_iterator it=cscSegments->begin(); it != cscSegments->end(); it++) {
01753 iSegment++;
01754 CSCDetId id = (CSCDetId)(*it).cscDetId();
01755 SegmentInChamber.push_back(id.chamber());
01756 SegmentInRing.push_back(id.ring());
01757 SegmentInStation.push_back(id.station());
01758
01759 std::vector <int> LayersInChamber(6);
01760 printf("\t iSegment = %i",iSegment);
01761 double chisq = (*it).chi2();
01762 int DOF = (*it).degreesOfFreedom();
01763 int nhits = (*it).nRecHits();
01764 LocalPoint localPos = (*it).localPosition();
01765 LocalVector localDir = (*it).localDirection();
01766 if (printalot){
01767 printf("\tendcap/station/ring/chamber: %i %i %i %i\n",
01768 id.endcap(),id.station(),id.ring(),id.chamber());
01769 cout<<"chi2/ndf = "<<chisq/DOF<<" nhits = "<<nhits <<std::endl;
01770 }
01771 Chi2->Fill(chisq/DOF);
01772 if(1==id.station()){
01773 if(1==id.ring()){
01774 Chi2_ME1_b->Fill(chisq/DOF);
01775 }
01776 else if(2==id.ring()){
01777 Chi2_ME1_2->Fill(chisq/DOF);
01778 }
01779 else if(3==id.ring()){
01780 Chi2_ME1_3->Fill(chisq/DOF);
01781 }
01782 else if(4==id.ring()){
01783 Chi2_ME1_a->Fill(chisq/DOF);
01784 }
01785 }
01786 else if(2==id.station()){
01787 if(1==id.ring()){
01788 Chi2_ME2_1->Fill(chisq/DOF);
01789 }
01790 else if(2==id.ring()){
01791 Chi2_ME2_2->Fill(chisq/DOF);
01792 }
01793 }
01794 else if(3==id.station()){
01795 if(1==id.ring()){
01796 Chi2_ME3_1->Fill(chisq/DOF);
01797 }
01798 else if(2==id.ring()){
01799 Chi2_ME3_2->Fill(chisq/DOF);
01800 }
01801 }
01802 else if(4==id.station()){
01803 if(1==id.ring()){
01804 Chi2_ME4_1->Fill(chisq/DOF);
01805 }
01806 }
01807
01808 if (printalot) printf("\tGet the recHits for this segment.\t");
01809 std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
01810 int nRH = (*it).nRecHits();
01811 if (printalot) printf(" nRH = %i\n",nRH);
01812 int jRH = 0;
01813 for ( vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
01814 jRH++;
01815 CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
01816 int kEndcap = idRH.endcap();
01817 int kRing = idRH.ring();
01818 int kStation = idRH.station();
01819 int kChamber = idRH.chamber();
01820 int kLayer = idRH.layer();
01821 LayersInChamber[kLayer-1] = 1;
01822
01823
01824 int iterations = 0;
01825 int RecHitCoincidence = 0;
01826 ChamberRecHits *sChamber_p=&AllRecHits[kEndcap-1][kStation-1][kRing-1][kChamber-FirstCh].sChamber;
01827 for(unsigned int hitsIn =0; hitsIn < (*sChamber_p).RecHitsPosXlocal[kLayer-1].size();hitsIn++){
01828
01829
01830 if( (*sChamber_p).RecHitsPosXlocal[kLayer-1][hitsIn] == (*iRH).localPosition().x() &&
01831 (*sChamber_p).RecHitsPosYlocal[kLayer-1][hitsIn] == (*iRH).localPosition().y() ){
01832 (*sChamber_p).TheRightRecHit[kLayer-1] = iterations;
01833 RecHitCoincidence++;
01834 }
01835 iterations++;
01836 }
01837 if(!RecHitCoincidence){
01838 (*sChamber_p).TheRightRecHit[kLayer-1] = -1;
01839 }
01840 if (printalot) printf("\t%i RH\tendcap/station/ring/chamber/layer: %i %i %i %i %i\n",jRH,kEndcap,kStation,kRing,kChamber,kLayer);
01841 if(printalot) std::cout<<"recHit number from the layer in the segment = "<< (*sChamber_p).TheRightRecHit[kLayer-1]<<std::endl;
01842 }
01843
01844
01845 double globX = 0.;
01846 double globY = 0.;
01847 double globZ = 0.;
01848 double globDirX = 0.;
01849 double globDirY = 0.;
01850 double globDirZ = 0.;
01851
01852 double globchamberPhi = 0.;
01853 double globchamberX =0.;
01854 double globchamberY =0.;
01855 double globchamberZ =0.;
01856
01857 const CSCChamber *cscchamber = cscGeom->chamber(id);
01858 if (cscchamber) {
01859 LocalPoint localCenter(0.,0.,0);
01860 GlobalPoint cscchamberCenter = cscchamber->toGlobal(localCenter);
01861 rotationPhi = globchamberPhi = cscchamberCenter.phi();
01862 globchamberX = cscchamberCenter.x();
01863 globchamberY = cscchamberCenter.y();
01864 globchamberZ = cscchamberCenter.z();
01865
01866 GlobalPoint globalPosition = cscchamber->toGlobal(localPos);
01867 globX = globalPosition.x();
01868 globY = globalPosition.y();
01869 globZ = globalPosition.z();
01870
01871 GlobalVector globalDirection = cscchamber->toGlobal(localDir);
01872 globDirX = globalDirection.x();
01873 globDirY = globalDirection.y();
01874 globDirZ = globalDirection.z();
01875 if(printalot) std::cout<<"SEGMENT: globDirX/globDirZ = "<<globDirX/globDirZ<<" globDirY/globDirZ = "<<globDirY/globDirZ<<std::endl;
01876 } else {
01877 if (printalot) printf("\tFailed to get a local->global segment tranformation.\n");
01878 }
01879
01880
01881 if (WorkInEndcap==id.endcap()) {
01882 if((ExtrapolateToStation==id.station()&& ExtrapolateToRing==id.ring() &&
01883 id.station()>=FirstCh &&id.station()<=LastCh) ||
01884 ExtrapolateFromStation==id.station()){
01885 int CoupleNum = id.station()-2;
01886
01887 if(id.station()==ExtrapolateToStation){
01888 CoupleNum = 1;
01889 }
01890 else if(id.station()==ExtrapolateFromStation){
01891 CoupleNum = 0;
01892 }
01893 else{
01894 cout<<"Wrong reference station!!! (shouldn't be)"<<endl;
01895 }
01896
01897 ChambersWithSegments[CoupleNum][LastCh-id.chamber()]++;
01898 PosLocalCouple[CoupleNum][0] = localPos.x();
01899 PosLocalCouple[CoupleNum][1] = localPos.y();
01900 PosLocalCouple[CoupleNum][2] = localPos.z();
01901
01902 DirLocalCouple[CoupleNum][0] = localDir.x();
01903 DirLocalCouple[CoupleNum][1] = localDir.y();
01904 DirLocalCouple[CoupleNum][2] = localDir.z();
01905
01906 PosCouple[CoupleNum][0] = globX;
01907 PosCouple[CoupleNum][1] = globY;
01908 PosCouple[CoupleNum][2] = globZ;
01909
01910 DirCouple[CoupleNum][0] = globDirX;
01911 DirCouple[CoupleNum][1] = globDirY;
01912 DirCouple[CoupleNum][2] = globDirZ;
01913
01914
01915
01916 RotPhiCouple[CoupleNum] = rotationPhi;
01917 Chi2Couple[CoupleNum] = chisq;
01918 NDFCouple[CoupleNum] = DOF;
01919 NhitsCouple[CoupleNum] = nhits;
01920
01921 XchamberCouple[CoupleNum] = globchamberX;
01922 YchamberCouple[CoupleNum] = globchamberY;
01923 NChamberCouple[CoupleNum] = id.chamber();
01924
01925 const CSCLayer *layer_p = cscchamber->layer(1);
01926 const CSCLayerGeometry *layerGeom = layer_p->geometry ();
01927 const std::vector<float> LayerBounds = layerGeom->parameters ();
01928 ChamberBoundsCouple[CoupleNum][0] = LayerBounds[0];
01929 ChamberBoundsCouple[CoupleNum][1] = LayerBounds[1];
01930 ChamberBoundsCouple[CoupleNum][2] = LayerBounds[3];
01931 if(ExtrapolateToStation==id.station()){
01932 LayersInSecondCouple = LayersInChamber;
01933 }
01934 }
01935 if(2 == ExtrapolateToStation &&
01936 (1==ExtrapolateFromStation || 3==ExtrapolateFromStation)){
01937 if(ExtrapolateFromStation != id.station() && ExtrapolateToStation != id.station() && 6==nhits && chisq/DOF<3){
01938 DirGlobal_ThirdSegment.push_back(globDirX);
01939 DirGlobal_ThirdSegment.push_back(globDirY);
01940 DirGlobal_ThirdSegment.push_back(globDirZ);
01941 }
01942 }
01943 }
01944 }
01945 printf("My nSegments: %i\n",nSegments);
01946
01947
01948 if(nSegments){
01949 DataFlow->Fill(2.);
01950 std::vector<int> SegmentsInStation(2);
01951 for(int iCh=0;iCh<NumCh;iCh++){
01952 SegmentsInStation[0]+=ChambersWithSegments[0][iCh];
01953 SegmentsInStation[1]+=ChambersWithSegments[1][iCh];
01954 }
01955
01956
01957 if(1==SegmentsInStation[0] ){
01958 DataFlow->Fill(4.);
01959
01960
01961 if(6==NhitsCouple[0] && (Chi2Couple[0]/NDFCouple[0])<3.){
01962
01963 DataFlow->Fill(6.);
01964 flag = false;
01965
01966 if(DirGlobal_ThirdSegment.size()){
01967 double XDirprime, YDirprime;
01968 Rotate(DirCouple[0][0], DirCouple[0][1], -(RotPhiCouple[1]+M_PI/2), XDirprime, YDirprime);
01969 double ZDirprime = DirCouple[0][2];
01970 double XDirsecond, YDirsecond;
01971 Rotate(DirGlobal_ThirdSegment[0], DirGlobal_ThirdSegment[1], -(RotPhiCouple[1]+M_PI/2), XDirsecond, YDirsecond);
01972 double ZDirsecond = DirGlobal_ThirdSegment[2];
01973 float diff_dxdz = XDirprime/ZDirprime - XDirsecond/ZDirsecond;
01974 if(fabs(diff_dxdz)<0.4){
01975 flag = true;
01976 seg_dydz = YDirprime/ZDirprime;
01977 }
01978 }
01979 int NSegFound = 0;
01980 double theta = acos(DirCouple[0][2]);
01981 for(int iSeg=0;iSeg<nSegments;iSeg++){
01982 if( NChamberCouple[1]==SegmentInChamber[iSeg]
01983 && ExtrapolateToRing==SegmentInRing[iSeg]
01984 && ExtrapolateToStation==SegmentInStation[iSeg]){
01985 NSegFound++;
01986 }
01987 }
01988
01989
01990 CalculateEfficiencies( event, eventSetup, PosCouple[0], DirCouple[0],NSegFound, theta);
01991
01992
01993 if(1==NSegFound){
01994 DataFlow->Fill(25.);
01995
01996
01997 for (int iLayer=0; iLayer<6; iLayer++) {
01998
01999 if(MyRecHits[NChamberCouple[1]-FirstCh][iLayer]){
02000
02001 if(!LayersInSecondCouple[iLayer]){
02002 ChHist[NChamberCouple[1]-FirstCh].InefficientSingleHits->Fill(iLayer+1);
02003 }
02004 ChHist[NChamberCouple[1]-FirstCh].AllSingleHits->Fill(iLayer+1);
02005 }
02006 }
02007
02008
02009
02010 double Xsecond, Ysecond;
02011 Rotate(PosCouple[1][0], PosCouple[1][1], -(RotPhiCouple[1]+M_PI/2), Xsecond, Ysecond);
02012 double XDirprime, YDirprime;
02013 Rotate(DirCouple[0][0], DirCouple[0][1], -(RotPhiCouple[1]+M_PI/2), XDirprime, YDirprime);
02014 double ZDirprime = DirCouple[0][2];
02015 double XDirsecond, YDirsecond;
02016 Rotate(DirCouple[1][0], DirCouple[1][1], -(RotPhiCouple[1]+M_PI/2), XDirsecond, YDirsecond);
02017 double ZDirsecond = DirCouple[1][2];
02018 double dxdz_diff = XDirprime/ZDirprime - XDirsecond/ZDirsecond;
02019
02020 if(abs(dxdz_diff)<0.15){
02021 DataFlow->Fill(26.);
02022 if(6!=NhitsCouple[1]){
02023 ChHist[NChamberCouple[1]-FirstCh].XvsY_InefficientSegments->Fill(PosLocalCouple[1][0],PosLocalCouple[1][1]);
02024 }
02025 double Xchambersecond, Ychambersecond;
02026 Rotate(XchamberCouple[1], YchamberCouple[1], -(RotPhiCouple[1]+M_PI/2.), Xchambersecond, Ychambersecond);
02027
02028 double Yup, Ydown, LineSlope , LineConst, Yright;
02029 Yup = Ychambersecond + ChamberBoundsCouple[1][2];
02030 Ydown = Ychambersecond - ChamberBoundsCouple[1][2];
02031 LineSlope = (Yup - Ydown)/(ChamberBoundsCouple[1][0]-ChamberBoundsCouple[1][1]);
02032 LineConst = Yup - LineSlope*ChamberBoundsCouple[1][0];
02033 Yright = LineSlope*abs(Xsecond) + LineConst;
02034 double XBound1Shifted = ChamberBoundsCouple[1][0]-20.;
02035 double XBound2Shifted = ChamberBoundsCouple[1][1]-20.;
02036 LineSlope = (Yup - Ydown)/(XBound1Shifted-XBound2Shifted);
02037 LineConst = Yup - LineSlope*XBound1Shifted;
02038 Yright = LineSlope*abs(Xsecond) + LineConst;
02039
02040
02041 if(GoodRegion(Ysecond, Yright, ExtrapolateToStation, ExtrapolateToRing, 0)){
02042 DataFlow->Fill(27.);
02043 if(6!=NhitsCouple[1]){
02044 ChHist[NChamberCouple[1]-FirstCh].XvsY_InefficientSegments_good->Fill(PosLocalCouple[1][0],PosLocalCouple[1][1]);
02045 }
02046 }
02047
02048 ChamberRecHits *sChamber_p =
02049 &(*all_RecHits)[WorkInEndcap-1][ExtrapolateToStation-1][ExtrapolateToRing-1][NChamberCouple[1]-FirstCh].sChamber;
02050
02051 double Zlayer = 0.;
02052 int ChosenLayer = -1;
02053 if(sChamber_p->RecHitsPosZ[3-1].size()){
02054 ChosenLayer = 2;
02055 Zlayer = sChamber_p->RecHitsPosZ[3-1][0];
02056 }
02057 else if(sChamber_p->RecHitsPosZ[4-1].size()){
02058 ChosenLayer = 3;
02059 Zlayer = sChamber_p->RecHitsPosZ[4-1][0];
02060 }
02061 else if(sChamber_p->RecHitsPosZ[5-1].size()){
02062 ChosenLayer = 4;
02063 Zlayer = sChamber_p->RecHitsPosZ[5-1][0];
02064 }
02065 else if(sChamber_p->RecHitsPosZ[2-1].size()){
02066 ChosenLayer = 1;
02067 Zlayer = sChamber_p->RecHitsPosZ[2-1][0];
02068 }
02069
02070 if(-1!=ChosenLayer){
02071
02072
02073 double dist = 2.54;
02074 for (int iLayer=0; iLayer<6; iLayer++) {
02075
02076
02077 double z1Position = PosCouple[1][2];
02078 double z2Position = Zlayer;
02079 double z1Direction = DirLocalCouple[1][2];
02080 double initPosition = PosLocalCouple[1][1];
02081 double initDirection = DirLocalCouple[1][1];
02082 double ParamLine = LineParam(z1Position, z2Position, z1Direction);
02083
02084
02085 double y = Extrapolate1D(initPosition, initDirection, ParamLine);
02086
02087
02088 initPosition = PosLocalCouple[1][0];
02089 initDirection = DirLocalCouple[1][0];
02090 double x = Extrapolate1D(initPosition, initDirection, ParamLine);
02091
02092 int sign;
02093
02094 if( z1Position>z2Position){
02095 sign = -1;
02096 }
02097 else{
02098 sign = 1;
02099 }
02100 z1Position = z2Position;
02101 int diffLayer = abs(ChosenLayer - iLayer);
02102 z2Position = z1Position + float(sign)*float(diffLayer)*dist;
02103
02104 ParamLine = LineParam(z1Position, z2Position, z1Direction);
02105 initPosition = y;
02106 initDirection = DirLocalCouple[1][1];
02107 y = Extrapolate1D(initPosition, initDirection, ParamLine);
02108
02109 initPosition = x;
02110 initDirection = DirLocalCouple[1][0];
02111 x = Extrapolate1D(initPosition, initDirection, ParamLine);
02112
02113 if(GoodRegion(Ysecond, Yright, ExtrapolateToStation, ExtrapolateToRing, 0)){
02114 if(sChamber_p->NRecHits[iLayer]>0){
02115 ChHist[NChamberCouple[1]-FirstCh].EfficientRechits_inSegment->Fill(iLayer+1);
02116 }
02117 else{
02118 ChHist[NChamberCouple[1]-FirstCh].XvsY_InefficientRecHits_inSegment[iLayer]->Fill(x,y);
02119 }
02120 }
02121 if(sChamber_p->NRecHits[iLayer]<1){
02122 ChHist[NChamberCouple[1]-FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Fill(y);
02123 }
02124 ChHist[NChamberCouple[1]-FirstCh].Y_AllRecHits_inSegment[iLayer]->Fill(y);
02125 }
02126
02127 if(GoodRegion(Ysecond, Yright, ExtrapolateToStation, ExtrapolateToRing, 0)){
02128 ChHist[NChamberCouple[1]-FirstCh].EfficientRechits_inSegment->Fill(9);
02129 }
02130 }
02131 }
02132 }
02133 }
02134 }
02135 }
02136
02137 if (printalot) printf("==exit===CSCEfficiency===== run %i\tevent %i\n\n",iRun,iEvent);
02138 }
02139
02140 void Rotate(double Xinit, double Yinit, double angle, double & Xrot, double & Yrot){
02141
02142 Xrot = Xinit*cos(angle) - Yinit*sin(angle);
02143 Yrot = Xinit*sin(angle) + Yinit*cos(angle);
02144 }
02145 bool CSCEfficiency::GoodRegion(double Y, double Yborder, int Station, int Ring, int Chamber){
02146
02147
02148 bool pass = false;
02149 double Ycenter = 99999.;
02150 float y_center = 99999.;
02151 if(Station>1 && Station<5){
02152 if(2==Ring){
02153 y_center = -0.95;
02154 Ycenter = 338.0/2+360.99-3.49+y_center;
02155 }
02156 else if(1==Ring){
02157 if(2==Station){
02158 y_center = -0.955;
02159 Ycenter = 204.60/2 + 143.89 - 3.49 + y_center;
02160 }
02161 else if(3==Station){
02162 y_center = -0.97;
02163 Ycenter = 184.61/2 + 163.89 - 3.49 + y_center;
02164 }
02165 else if(4==Station){
02166 y_center = -0.94;
02167 Ycenter = 164.67/2 + 183.79 - 3.49 + y_center;
02168 }
02169 }
02170 }
02171 else{
02172 if(3==Ring){
02173 y_center = -1.075;
02174 Ycenter = 179.30/2+508.99-3.49+y_center;
02175 }
02176 else if(2==Ring){
02177 y_center = -0.96;
02178 Ycenter = 189.41/2+278.49-3.49+y_center;
02179 }
02180 else if(1==Ring){
02181 Ycenter = 0.;
02182 }
02183 }
02184 Ycenter = -Ycenter;
02185 double Y_local = -(Y - Ycenter);
02186 double Yborder_local = -(Yborder - Ycenter);
02187 bool withinChamberOnly = false;
02188 pass = CheckLocal(Y_local, Yborder_local, Station, Ring, withinChamberOnly);
02189 return pass;
02190 }
02191 bool CSCEfficiency::GoodLocalRegion(double X, double Y, int Station, int Ring, int Chamber){
02192
02193 bool pass = false;
02194 std::vector <double> ChamberBoundsCouple(3);
02195 float y_center = 99999.;
02196
02197 if(Station>1 && Station<5){
02198 if(2==Ring){
02199 ChamberBoundsCouple[0] = 66.46/2;
02200 ChamberBoundsCouple[1] = 127.15/2;
02201 ChamberBoundsCouple[2] = 323.06/2;
02202 y_center = -0.95;
02203 }
02204 else{
02205 if(2==Station){
02206 ChamberBoundsCouple[0] = 54.00/2;
02207 ChamberBoundsCouple[1] = 125.71/2;
02208 ChamberBoundsCouple[2] = 189.66/2;
02209 y_center = -0.955;
02210 }
02211 else if(3==Station){
02212 ChamberBoundsCouple[0] = 61.40/2;
02213 ChamberBoundsCouple[1] = 125.71/2;
02214 ChamberBoundsCouple[2] = 169.70/2;
02215 y_center = -0.97;
02216 }
02217 else if(4==Station){
02218 ChamberBoundsCouple[0] = 69.01/2;
02219 ChamberBoundsCouple[1] = 125.65/2;
02220 ChamberBoundsCouple[2] = 149.42/2;
02221 y_center = -0.94;
02222 }
02223 }
02224 }
02225 else if(1==Station){
02226 if(3==Ring){
02227 ChamberBoundsCouple[0] = 63.40/2;
02228 ChamberBoundsCouple[1] = 92.10/2;
02229 ChamberBoundsCouple[2] = 164.16/2;
02230 y_center = -1.075;
02231 }
02232 else if(2==Ring){
02233 ChamberBoundsCouple[0] = 51.00/2;
02234 ChamberBoundsCouple[1] = 83.74/2;
02235 ChamberBoundsCouple[2] = 174.49/2;
02236 y_center = -0.96;
02237 }
02238 else{
02239 ChamberBoundsCouple[0] = 40./2;
02240 ChamberBoundsCouple[1] = 100./2;
02241 ChamberBoundsCouple[2] = 142./2;
02242 y_center = 0.;
02243 }
02244 }
02245 double Yup = ChamberBoundsCouple[2] + y_center;
02246 double Ydown = - ChamberBoundsCouple[2] + y_center;
02247 double XBound1Shifted = ChamberBoundsCouple[0]-20.;
02248 double XBound2Shifted = ChamberBoundsCouple[1]-20.;
02249 double LineSlope = (Yup - Ydown)/(XBound2Shifted-XBound1Shifted);
02250 double LineConst = Yup - LineSlope*XBound2Shifted;
02251 double Yborder = LineSlope*abs(X) + LineConst;
02252 bool withinChamberOnly = false;
02253 pass = CheckLocal(Y, Yborder, Station, Ring, withinChamberOnly);
02254 return pass;
02255 }
02256
02257 bool CSCEfficiency::CheckLocal(double Y, double Yborder, int Station, int Ring, bool withinChamberOnly){
02258
02259 bool pass = false;
02260
02261 std::vector <float> DeadZoneCenter(6);
02262 float CutZone = 10.;
02263
02264 if(!withinChamberOnly){
02265 if(Station>1 && Station<5){
02266 if(2==Ring){
02267 DeadZoneCenter[0]= -162.48 ;
02268 DeadZoneCenter[1] = -81.8744;
02269 DeadZoneCenter[2] = -21.18165;
02270 DeadZoneCenter[3] = 39.51105;
02271 DeadZoneCenter[4] = 100.2939;
02272 DeadZoneCenter[5] = 160.58;
02273
02274 if(Y >Yborder &&
02275 ((Y> DeadZoneCenter[0] + CutZone && Y< DeadZoneCenter[1] - CutZone) ||
02276 (Y> DeadZoneCenter[1] + CutZone && Y< DeadZoneCenter[2] - CutZone) ||
02277 (Y> DeadZoneCenter[2] + CutZone && Y< DeadZoneCenter[3] - CutZone) ||
02278 (Y> DeadZoneCenter[3] + CutZone && Y< DeadZoneCenter[4] - CutZone) ||
02279 (Y> DeadZoneCenter[4] + CutZone && Y< DeadZoneCenter[5] - CutZone))){
02280 pass = true;
02281 }
02282 }
02283 else if(1==Ring){
02284 if(2==Station){
02285 DeadZoneCenter[0]= -95.80 ;
02286 DeadZoneCenter[1] = -27.47;
02287 DeadZoneCenter[2] = 33.67;
02288 DeadZoneCenter[3] = 90.85;
02289 }
02290 else if(3==Station){
02291 DeadZoneCenter[0]= -89.305 ;
02292 DeadZoneCenter[1] = -39.705;
02293 DeadZoneCenter[2] = 20.195;
02294 DeadZoneCenter[3] = 77.395;
02295 }
02296 else if(4==Station){
02297 DeadZoneCenter[0]= -75.645;
02298 DeadZoneCenter[1] = -26.055;
02299 DeadZoneCenter[2] = 23.855;
02300 DeadZoneCenter[3] = 70.575;
02301 }
02302 if(Y >Yborder &&
02303 ((Y> DeadZoneCenter[0] + CutZone && Y< DeadZoneCenter[1] - CutZone) ||
02304 (Y> DeadZoneCenter[1] + CutZone && Y< DeadZoneCenter[2] - CutZone) ||
02305 (Y> DeadZoneCenter[2] + CutZone && Y< DeadZoneCenter[3] - CutZone))){
02306 pass = true;
02307 }
02308
02309 }
02310 }
02311 else if(1==Station){
02312 if(3==Ring){
02313 DeadZoneCenter[0]= -83.155 ;
02314 DeadZoneCenter[1] = -22.7401;
02315 DeadZoneCenter[2] = 27.86665;
02316 DeadZoneCenter[3] = 81.005;
02317 if(Y > Yborder &&
02318 ((Y> DeadZoneCenter[0] + CutZone && Y< DeadZoneCenter[1] - CutZone) ||
02319 (Y> DeadZoneCenter[1] + CutZone && Y< DeadZoneCenter[2] - CutZone) ||
02320 (Y> DeadZoneCenter[2] + CutZone && Y< DeadZoneCenter[3] - CutZone))){
02321 pass = true;
02322 }
02323 }
02324 else if(2==Ring){
02325 DeadZoneCenter[0]= -86.285 ;
02326 DeadZoneCenter[1] = -32.88305;
02327 DeadZoneCenter[2] = 32.867423;
02328 DeadZoneCenter[3] = 88.205;
02329 if(Y > (Yborder) &&
02330 ((Y> DeadZoneCenter[0] + CutZone && Y< DeadZoneCenter[1] - CutZone) ||
02331 (Y> DeadZoneCenter[1] + CutZone && Y< DeadZoneCenter[2] - CutZone) ||
02332 (Y> DeadZoneCenter[2] + CutZone && Y< DeadZoneCenter[3] - CutZone))){
02333 pass = true;
02334 }
02335 }
02336 else{
02337 DeadZoneCenter[0]= -81.0;
02338 DeadZoneCenter[1] = 81.0;
02339 if(Y > (Yborder) &&
02340 ((Y> DeadZoneCenter[0] + CutZone && Y< DeadZoneCenter[1] - CutZone) )){
02341 pass = true;
02342 }
02343 }
02344 }
02345 }
02346 else{
02347 if(Station>1 && Station<5){
02348 if(2==Ring){
02349 if(Y >Yborder && fabs(Y+0.95)<151.53){
02350 pass = true;
02351 }
02352 }
02353 else if(1==Ring){
02354 float theCenter = 0.95;
02355 float theCutEdge = 151.53;
02356 if(2==Station){
02357 theCenter = - 0.955;
02358 theCutEdge = 84.83;
02359 }
02360 else if(3==Station){
02361 theCenter = - 0.97;
02362 theCutEdge = 74.85;
02363 }
02364 else if(4==Station){
02365 theCenter = -0.94;
02366 theCutEdge = 64.71;
02367 }
02368 if(Y >Yborder && fabs(Y-theCenter)<theCutEdge){
02369 pass = true;
02370 }
02371 }
02372 }
02373 else if(1==Station){
02374 if(3==Ring){
02375 if(Y > Yborder && fabs(Y+1.075)<72.08){
02376 pass = true;
02377 }
02378 }
02379 else if(2==Ring){
02380 if(Y > (Yborder) && fabs(Y+0.96)<77.245){
02381 pass = true;
02382 }
02383 }
02384 else{
02385 if(Y > (Yborder) && fabs(Y)<71.){
02386 pass = true;
02387 }
02388 }
02389 }
02390 }
02391 return pass;
02392 }
02393
02394 void CSCEfficiency::CalculateEfficiencies(const Event & event, const EventSetup& eventSetup,
02395 std::vector<double> &Pos , std::vector<double> &Dir, int NSegFound, double theta){
02396 DataFlow->Fill(7.);
02397 edm::Handle<CSCALCTDigiCollection> alcts;
02398 edm::Handle<CSCCLCTDigiCollection> clcts;
02399
02400 edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts;
02401 if(DATA){
02402 event.getByLabel("muonCSCDigis","MuonCSCALCTDigi",alcts);
02403 event.getByLabel("muonCSCDigis","MuonCSCCLCTDigi",clcts);
02404
02405 event.getByLabel("muonCSCDigis","MuonCSCCorrelatedLCTDigi",correlatedlcts);
02406 }
02407 ESHandle<CSCGeometry> cscGeom;
02408 eventSetup.get<MuonGeometryRecord>().get(cscGeom);
02409
02410 const std::vector<CSCChamber*> ChamberContainer = cscGeom->chambers();
02411
02412 if(!DATA){
02413 RecSimHitEfficiency();
02414 }
02415
02416
02417 for(unsigned int nCh=0;nCh<ChamberContainer.size();nCh++){
02418 const CSCChamber *cscchamber = ChamberContainer[nCh];
02419 CSCDetId id = cscchamber->id();
02420 if(id.chamber() > (FirstCh-1) && id.chamber() < (LastCh+1) &&
02421 id.station() == ExtrapolateToStation &&
02422 id.ring() == ExtrapolateToRing && id.endcap() == WorkInEndcap){
02423 LocalPoint localCenter(0.,0.,0);
02424 GlobalPoint cscchamberCenter = cscchamber->toGlobal(localCenter);
02425 float ZStation = cscchamberCenter.z();
02426 double ParLine = LineParam(Pos[2], ZStation, Dir[2]);
02427 double Xextrapolated = Extrapolate1D(Pos[0],Dir[0], ParLine );
02428 double Yextrapolated = Extrapolate1D(Pos[1],Dir[1], ParLine );
02429
02430 GlobalPoint ExtrapolatedSegment(Xextrapolated, Yextrapolated, ZStation);
02431 LocalPoint ExtrapolatedSegmentLocal = cscchamber->toLocal(ExtrapolatedSegment);
02432 const CSCLayer *layer_p = cscchamber->layer(1);
02433 const CSCLayerGeometry *layerGeom = layer_p->geometry ();
02434 const std::vector<float> LayerBounds = layerGeom->parameters ();
02435 float y_center = 0.;
02436 float ShiftFromEdge = 20.;
02437 double Yup = LayerBounds[3] + y_center;
02438 double Ydown = - LayerBounds[3] + y_center;
02439 double XBound1Shifted = LayerBounds[0] - ShiftFromEdge;
02440 double XBound2Shifted = LayerBounds[1] - ShiftFromEdge;
02441 double LineSlope = (Yup - Ydown)/(XBound2Shifted-XBound1Shifted);
02442 double LineConst = Yup - LineSlope*XBound2Shifted;
02443 double Yborder = LineSlope*abs(ExtrapolatedSegmentLocal.x()) + LineConst;
02444 CSCDetId id = cscchamber->id();
02445 bool withinChamberOnly = false;
02446 if( CheckLocal(ExtrapolatedSegmentLocal.y(), Yborder, id.station(), id.ring(), withinChamberOnly)){
02447 DataFlow->Fill(9.);
02448 int cond = 0;
02449 if(flag){
02450 cond = 1;
02451 }
02452
02453
02454
02455
02456 bool LCTflag = false;
02457 if(DATA){
02458 LCTflag = LCT_Efficiencies(alcts, clcts, correlatedlcts, id.chamber(), cond);
02459 }
02460 if(!LCTflag){
02461 XY_ALCTmissing->Fill(ExtrapolatedSegmentLocal.x(),ExtrapolatedSegmentLocal.y());
02462 std::cout<<"NO ALCT when ME1 and ME3"<<std::endl;
02463 }
02464 StripWire_Efficiencies(id.chamber());
02465 Segment_Efficiency(id.chamber(), NSegFound, theta);
02466 }
02467
02468
02469 withinChamberOnly = true;
02470 ShiftFromEdge = 10.;
02471 XBound1Shifted = LayerBounds[0] - ShiftFromEdge;
02472 XBound2Shifted = LayerBounds[1] - ShiftFromEdge;
02473 LineSlope = (Yup - Ydown)/(XBound2Shifted-XBound1Shifted);
02474 LineConst = Yup - LineSlope*XBound2Shifted;
02475 Yborder = LineSlope*abs(ExtrapolatedSegmentLocal.x()) + LineConst;
02476 if( CheckLocal(ExtrapolatedSegmentLocal.y(), Yborder, id.station(), id.ring(), withinChamberOnly)){
02477 RecHitEfficiency(ExtrapolatedSegmentLocal.x() , ExtrapolatedSegmentLocal.y(), id.chamber());
02478 }
02479 }
02480 }
02481 }
02482
02483 double CSCEfficiency::Extrapolate1D(double initPosition, double initDirection, double ParameterOfTheLine){
02484 double ExtrapolatedPosition = initPosition + initDirection*ParameterOfTheLine;
02485 return ExtrapolatedPosition;
02486 }
02487
02488 double CSCEfficiency::LineParam(double z1Position, double z2Position, double z1Direction){
02489 double ParamLine = (z2Position-z1Position)/z1Direction;
02490 return ParamLine;
02491 }
02492
02493 void CSCEfficiency::RecHitEfficiency(double X, double Y, int iCh){
02494 ChamberRecHits *sChamber_p =
02495 &(*all_RecHits)[WorkInEndcap-1][ExtrapolateToStation-1][ExtrapolateToRing-1][iCh-FirstCh].sChamber;
02496 ChamberRecHits *sChamberLeft_p = sChamber_p;
02497 ChamberRecHits *sChamberRight_p = sChamber_p;
02498
02499
02500 if(iCh-FirstCh-1>=0){
02501 sChamberLeft_p =
02502 &(*all_RecHits)[WorkInEndcap-1][ExtrapolateToStation-1][ExtrapolateToRing-1][iCh-FirstCh-1].sChamber;
02503 }
02504
02505 if(iCh-FirstCh+1<NumCh){
02506 sChamberRight_p =
02507 &(*all_RecHits)[WorkInEndcap-1][ExtrapolateToStation-1][ExtrapolateToRing-1][iCh-FirstCh+1].sChamber;
02508 }
02509 int missingLayersLeft = 0;
02510 int missingLayersRight = 0;
02511 int missingLayers = 0;
02512 for(int iLayer=0;iLayer<6;iLayer++){
02513 if(!sChamber_p->RecHitsPosX[iLayer].size()){
02514 missingLayers++;
02515 }
02516 if(!sChamberLeft_p->RecHitsPosX[iLayer].size()){
02517 missingLayersLeft++;
02518 }
02519 if(!sChamberRight_p->RecHitsPosX[iLayer].size()){
02520 missingLayersRight++;
02521 }
02522 }
02523 if(missingLayers>missingLayersLeft || missingLayers>missingLayersRight){
02524
02525 }
02526 else{
02527
02528 if(GoodLocalRegion(X, Y, ExtrapolateToStation, ExtrapolateToRing, iCh)){
02529 if(6==missingLayers){
02530 if(printalot) std::cout<<"missing all layers"<<std::endl;
02531 }
02532 for(int iLayer=0;iLayer<6;iLayer++){
02533 if(missingLayers){
02534 for(unsigned int iRH=0; iRH<sChamber_p->RecHitsPosX[iLayer].size();iRH++){
02535 ChHist[iCh-FirstCh].XvsY_InefficientRecHits->Fill(sChamber_p->RecHitsPosXlocal[iLayer][iRH],
02536 sChamber_p->RecHitsPosYlocal[iLayer][iRH]);
02537 }
02538 }
02539
02540 if(sChamber_p->NRecHits[iLayer]>0){
02541 ChHist[iCh-FirstCh].EfficientRechits->Fill(iLayer+1);
02542 }
02543 }
02544 if(6!=missingLayers){
02545 ChHist[iCh-FirstCh].EfficientRechits->Fill(8);
02546 }
02547 ChHist[iCh-FirstCh].EfficientRechits->Fill(9);
02548 }
02549
02550 int badhits = 0;
02551 int realhit = 0;
02552 for(int iLayer=0;iLayer<6;iLayer++){
02553 for(unsigned int iRH=0; iRH<sChamber_p->RecHitsPosX[iLayer].size();iRH++){
02554 realhit++;
02555
02556 if(!GoodLocalRegion(sChamber_p->RecHitsPosXlocal[iLayer][iRH],
02557 sChamber_p->RecHitsPosYlocal[iLayer][iRH],
02558 ExtrapolateToStation, ExtrapolateToRing, iCh)){
02559 badhits++;
02560 }
02561 }
02562 }
02563
02564
02565 if(0!=realhit && 0==badhits ){
02566 if(printalot) std::cout<<"good rechits"<<std::endl;
02567 for(int iLayer=0;iLayer<6;iLayer++){
02568 if(missingLayers){
02569 for(unsigned int iRH=0; iRH<sChamber_p->RecHitsPosX[iLayer].size();iRH++){
02570 ChHist[iCh-FirstCh].XvsY_InefficientRecHits_good->Fill(sChamber_p->RecHitsPosXlocal[iLayer][iRH],
02571 sChamber_p->RecHitsPosYlocal[iLayer][iRH]);
02572 }
02573 }
02574 if(sChamber_p->NRecHits[iLayer]>0){
02575 ChHist[iCh-FirstCh].EfficientRechits_good->Fill(iLayer+1);
02576 }
02577 }
02578 if(6!=missingLayers){
02579 ChHist[iCh-FirstCh].EfficientRechits_good->Fill(8);
02580 }
02581 ChHist[iCh-FirstCh].EfficientRechits_good->Fill(9);
02582 }
02583 }
02584 }
02585
02586 void CSCEfficiency::RecSimHitEfficiency(void){
02587 for(int iCh=0;iCh<LastCh;iCh++){
02588
02589 ChamberRecHits *sChamber_p =
02590 &(*all_RecHits)[WorkInEndcap-1][ExtrapolateToStation-1][ExtrapolateToRing-1][iCh].sChamber;
02591
02592 ChamberSimHits *sSimChamber_p =
02593 &(*all_SimHits)[WorkInEndcap-1][ExtrapolateToStation-1][ExtrapolateToRing-1][iCh].sChamber;
02594 for(int iLayer=0; iLayer<6;iLayer++){
02595 if(sSimChamber_p->NSimHits[iLayer]){
02596 ChHist[iCh].SimSimhits->Fill(iLayer+1);
02597 if(sChamber_p->NRecHits[iLayer]){
02598 ChHist[iCh].SimRechits->Fill(iLayer+1);
02599 }
02600
02601 for(unsigned int iSimHits=0;iSimHits<sSimChamber_p->SimHitsPosXlocal[iLayer].size();iSimHits++){
02602 ChHist[iCh].SimSimhits_each->Fill(iLayer+1);
02603 }
02604 for(unsigned int iRecHits=0;iRecHits<sChamber_p->RecHitsPosXlocal[iLayer].size();iRecHits++){
02605 ChHist[iCh].SimRechits_each->Fill(iLayer+1);
02606 }
02607
02608 }
02609 }
02610 }
02611 }
02612
02613 bool CSCEfficiency::LCT_Efficiencies(edm::Handle<CSCALCTDigiCollection> alcts,
02614 edm::Handle<CSCCLCTDigiCollection> clcts,
02615 edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts, int iCh, int cond){
02616 bool result = true;
02617 int Nalcts = 0;
02618 int Nclcts = 0;
02619 int Ncorr = 0;
02620 int Nalcts_1ch = 0;
02621 int Nclcts_1ch = 0;
02622 int Ncorr_1ch = 0;
02623
02624
02625 for (CSCALCTDigiCollection::DigiRangeIterator j=alcts->begin(); j!=alcts->end(); j++) {
02626 const CSCDetId& id = (*j).first;
02627 const CSCALCTDigiCollection::Range& range =(*j).second;
02628 for (CSCALCTDigiCollection::const_iterator digiIt =
02629 range.first; digiIt!=range.second;
02630 ++digiIt){
02631
02632
02633 if(0!=(*digiIt).isValid() &&
02634 id.station()==ExtrapolateToStation && id.ring()==ExtrapolateToRing){
02635 if(id.chamber()==iCh){
02636
02637
02638 Nalcts++;
02639 Nalcts_1ch++;
02640 }
02641 else if(1==abs(id.chamber()-iCh)){
02642 Nalcts++;
02643 }
02644 }
02645 }
02646 }
02647
02648
02649 for (CSCCLCTDigiCollection::DigiRangeIterator j=clcts->begin(); j!=clcts->end(); j++) {
02650 const CSCDetId& id = (*j).first;
02651 std::vector<CSCCLCTDigi>::const_iterator digiIt = (*j).second.first;
02652 std::vector<CSCCLCTDigi>::const_iterator last = (*j).second.second;
02653 for( ; digiIt != last; ++digiIt) {
02654
02655
02656 if(0!=(*digiIt).isValid() &&
02657 id.station()==ExtrapolateToStation && id.ring()==ExtrapolateToRing){
02658 if(id.chamber()==iCh){
02659
02660
02661 Nclcts++;
02662 Nclcts_1ch++;
02663 }
02664 else if(1==abs(id.chamber()-iCh)){
02665 Nclcts++;
02666 }
02667 }
02668 }
02669 }
02670
02671
02672 for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator j=correlatedlcts->begin(); j!=correlatedlcts->end(); j++) {
02673 const CSCDetId& id = (*j).first;
02674 std::vector<CSCCorrelatedLCTDigi>::const_iterator digiIt = (*j).second.first;
02675 std::vector<CSCCorrelatedLCTDigi>::const_iterator last = (*j).second.second;
02676 for( ; digiIt != last; ++digiIt) {
02677
02678
02679 if(0!=(*digiIt).isValid() &&
02680 id.station()==ExtrapolateToStation && id.ring()==ExtrapolateToRing){
02681 if(id.chamber()==iCh){
02682
02683
02684 Ncorr++;
02685 Ncorr_1ch++;
02686 }
02687 else if(1==abs(id.chamber()-iCh)){
02688 Ncorr++;
02689 }
02690 }
02691 }
02692 }
02693
02694 if(Nalcts){
02695 if(!Nclcts){
02696 if(printalot) std::cout<<"No alct-clct coincidence!"<<std::endl;
02697 }
02698
02699 ChHist[iCh-FirstCh].EfficientLCTs->Fill(1);
02700 if(Nalcts_1ch){
02701 ChHist[iCh-FirstCh].EfficientLCTs->Fill(11);
02702 if(cond){
02703 ChHist[iCh-FirstCh].EfficientLCTs->Fill(21);
02704 dydz_Eff_ALCT->Fill(seg_dydz);
02705 }
02706 }
02707 }
02708 else{
02709 if(cond){
02710 result = false;
02711 }
02712
02713 }
02714
02715 if(Nclcts){
02716 ChHist[iCh-FirstCh].EfficientLCTs->Fill(3);
02717 if(Nclcts_1ch){
02718 ChHist[iCh-FirstCh].EfficientLCTs->Fill(13);
02719 if(cond){
02720 ChHist[iCh-FirstCh].EfficientLCTs->Fill(23);
02721 }
02722 }
02723 }
02724
02725 if(Ncorr){
02726 ChHist[iCh-FirstCh].EfficientLCTs->Fill(5);
02727 if(Ncorr_1ch){
02728 ChHist[iCh-FirstCh].EfficientLCTs->Fill(15);
02729 if(cond){
02730 ChHist[iCh-FirstCh].EfficientLCTs->Fill(25);
02731 }
02732 }
02733 }
02734 ChHist[iCh-FirstCh].EfficientLCTs->Fill(30);
02735 if(cond){
02736 if(!Nalcts) if(printalot) std::cout<<"NO ALCT!"<<std::endl;
02737 ChHist[iCh-FirstCh].EfficientLCTs->Fill(28);
02738 dydz_All_ALCT->Fill(seg_dydz);
02739 }
02740 return result;
02741 }
02742
02743 void CSCEfficiency::StripWire_Efficiencies( int iCh){
02744 int EndCap = WorkInEndcap -1;
02745 int Ring = ExtrapolateToRing - 1;
02746 int Station = ExtrapolateToStation - 1;
02747
02748 int missingLayers_s = 0;
02749 int missingLayers_wg = 0;
02750 int badhits_s = 0;
02751 int badhits_wg = 0;
02752
02753
02754 for(int iLayer=0;iLayer<6;iLayer++){
02755 if(!AllStrips[EndCap][Station][Ring][iCh-FirstCh][iLayer].size()){
02756 missingLayers_s++;
02757 }
02758 for(unsigned int iStrip=0;
02759 iStrip<AllStrips[EndCap][Station][Ring][iCh-FirstCh][iLayer].size();
02760 iStrip++){
02761
02762 int Nstrips =80;
02763 int Step = 10;
02764 if(1==ExtrapolateToStation && 3==ExtrapolateToRing){
02765 Nstrips =64;
02766 }
02767
02768 if(AllStrips[EndCap][Station][Ring][iCh-FirstCh][iLayer][iStrip].first<Step ||
02769 AllStrips[EndCap][Station][Ring][iCh-FirstCh][iLayer][iStrip].first>(Nstrips-Step) ){
02770 badhits_s++;
02771 }
02772 }
02773 }
02774
02775
02776 for(int iLayer=0;iLayer<6;iLayer++){
02777 if(!AllWG[EndCap][Station][Ring][iCh-FirstCh][iLayer].size()){
02778 missingLayers_wg++;
02779 }
02780 for(unsigned int iWG=0;
02781 iWG<AllWG[EndCap][Station][Ring][iCh-FirstCh][iLayer].size();
02782 iWG++){
02783 if(!GoodLocalRegion(0.,
02784 AllWG[EndCap][Station][Ring][iCh-FirstCh][iLayer][iWG].first.second,
02785 ExtrapolateToStation, ExtrapolateToRing, iCh)){
02786 badhits_wg++;
02787 }
02788 }
02789 }
02790
02791
02792
02793
02794
02795 if(0==badhits_s && 0==badhits_wg){
02796 if(printalot){
02797 if(6!=missingLayers_wg){
02798 if(printalot) std::cout<<"good strips"<<std::endl;
02799 }
02800 }
02801 for(int iLayer=0;iLayer<6;iLayer++){
02802 if(AllStrips[EndCap][Station][Ring][iCh-FirstCh][iLayer].size()>0){
02803 ChHist[iCh-FirstCh].EfficientStrips->Fill(iLayer+1);
02804 }
02805 else{
02806 if(6!=missingLayers_s){
02807 if(printalot) std::cout<<"missing strips iLayer+1 = "<<iLayer+1<<std::endl;
02808 }
02809 }
02810 }
02811 if(6!=missingLayers_s){
02812 ChHist[iCh-FirstCh].EfficientStrips->Fill(8);
02813 }
02814 ChHist[iCh-FirstCh].EfficientStrips->Fill(9);
02815 }
02816
02817
02818
02819 if(0==badhits_wg && 0==badhits_s){
02820 if(printalot){
02821 if(6!=missingLayers_wg){
02822 if(printalot) std::cout<<"good WG"<<std::endl;
02823 }
02824 }
02825 for(int iLayer=0;iLayer<6;iLayer++){
02826 if(AllWG[EndCap][Station][Ring][iCh-FirstCh][iLayer].size()>0){
02827 ChHist[iCh-FirstCh].EfficientWireGroups->Fill(iLayer+1);
02828 }
02829 else{
02830 if(6!=missingLayers_wg){
02831 if(printalot) std::cout<<"missing WG iLayer+1 = "<<iLayer+1<<std::endl;
02832 }
02833 }
02834 }
02835 if(6!=missingLayers_wg){
02836 ChHist[iCh-FirstCh].EfficientWireGroups->Fill(8);
02837 }
02838 ChHist[iCh-FirstCh].EfficientWireGroups->Fill(9);
02839 }
02840 }
02841
02842 void CSCEfficiency::Segment_Efficiency(int iCh, int NSegmentsFound, double theta){
02843 ChamberRecHits *sChamber_p =
02844 &(*all_RecHits)[WorkInEndcap-1][ExtrapolateToStation-1][ExtrapolateToRing-1][iCh-FirstCh].sChamber;
02845 int missingLayers = 0;
02846 for(int iLayer=0;iLayer<6;iLayer++){
02847 if(!sChamber_p->RecHitsPosX[iLayer].size()){
02848 missingLayers++;
02849 }
02850 }
02851 if(missingLayers<5){
02852 if(NSegmentsFound){
02853 EfficientSegments->Fill(iCh);
02854 EfficientSegments_theta->Fill(theta);
02855 }
02856 AllSegments->Fill(iCh);
02857 AllSegments_theta->Fill(theta);
02858 }
02859 }
02860
02861 void CSCEfficiency::getEfficiency(float bin, float Norm, std::vector<float> &eff){
02862
02863 float Efficiency = 0.;
02864 float EffError = 0.;
02865 if(fabs(Norm)>0.000000001){
02866 Efficiency = bin/Norm;
02867 if(bin<Norm){
02868 EffError = sqrt( (1.-Efficiency)*Efficiency/Norm );
02869 }
02870 }
02871 eff[0] = Efficiency;
02872 eff[1] = EffError;
02873 }
02874
02875 void CSCEfficiency::histoEfficiency(TH1F *readHisto, TH1F *writeHisto, int flag){
02876 std::vector<float> eff(2);
02877 int Nbins = readHisto->GetSize()-2;
02878 std::vector<float> bins(Nbins);
02879 std::vector<float> Efficiency(Nbins);
02880 std::vector<float> EffError(Nbins);
02881 float Norm = 1;
02882 if(flag<=Nbins){
02883 Norm = readHisto->GetBinContent(flag);;
02884 }
02885 for (int i=0;i<Nbins;i++){
02886 bins[i] = readHisto->GetBinContent(i+1);
02887 getEfficiency(bins[i], Norm, eff);
02888 Efficiency[i] = eff[0];
02889 EffError[i] = eff[1];
02890 writeHisto->SetBinContent(i+1, Efficiency[i]);
02891 writeHisto->SetBinError(i+1, EffError[i]);
02892 }
02893 }
02894
02895 const char* CSCEfficiency::ChangeTitle(const char * name){
02896 std::string str = to_string(name);
02897 std::string searchString( "Ch1" );
02898 std::string replaceString( "AllCh" );
02899
02900 assert( searchString != replaceString );
02901
02902 std::string::size_type pos = 0;
02903 while ( (pos = str.find(searchString, pos)) != string::npos ) {
02904 str.replace( pos, searchString.size(), replaceString );
02905 pos++;
02906 }
02907 const char* NewName = str.c_str();
02908 return NewName;
02909 }
02910 DEFINE_FWK_MODULE(CSCEfficiency);
02911
02912
1.5.4