Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:26:03

0001 {
0002 #include <sstream>
0003 #include <iostream>
0004 #include <string.h>
0005 #include <vector.h>
0006 #include <math.h>
0007 
0008 gROOT->Reset();
0009 gROOT->SetStyle("Plain"); // to get rid of gray color of pad and have it white
0010 gStyle->SetPalette(1,0); // 
0011 std::ostringstream ss,ss1;
0012 
0013 /// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
0014 
0015 // select run number
0016 std::string run="62232";
0017 Int_t nminentries=50;
0018 
0019 //  input file with histograms
0020 ss.str("");
0021 ss<<"validationHists_"<<run<<".root";
0022 TFile f_in(ss.str().c_str());
0023 
0024 // output file with histograms
0025 ss.str("");
0026 ss<<"gas_gain_fit_"<<run<<".root";
0027 TFile f_out(ss.str().c_str(),"RECREATE");
0028 
0029 // folder in input file
0030 std::string in_folder="GasGain/";
0031 
0032 /// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
0033 
0034 /// common names and titles
0035 
0036 // for CSC in given ME+- station/ring
0037 std::vector<std::string> xTitle_ME_CSC;
0038 xTitle_ME_CSC.push_back("ME-1/1 and ME+1/1 CSC #");
0039 xTitle_ME_CSC.push_back("ME-1/2 and ME+1/2 CSC #");
0040 xTitle_ME_CSC.push_back("ME-1/3 and ME+1/3 CSC #"); 
0041 xTitle_ME_CSC.push_back("ME-2/1 and ME+2/1 CSC #");
0042 xTitle_ME_CSC.push_back("ME-2/2 and ME+2/2 CSC #");
0043 xTitle_ME_CSC.push_back("ME-3/1 and ME+3/1 CSC #");
0044 xTitle_ME_CSC.push_back("ME-3/2 and ME+3/2 CSC #");
0045 xTitle_ME_CSC.push_back("ME-4/1 and ME+4/1 CSC #");
0046 xTitle_ME_CSC.push_back("ME-4/2 and ME+4/2 CSC #");
0047 
0048 std::vector<std::string> xLabel_A_ME_CSC, xLabel_B_ME_CSC;
0049 Int_t flag=2;
0050 for(Int_t i=1;i<=36;i++) {
0051   
0052    ss.str("");
0053    if(flag==2) {
0054    Int_t cscnmb;
0055    if(i<19) cscnmb=i-19;
0056    if(i>18) cscnmb=i-18;
0057    if(cscnmb !=1) ss<<cscnmb;
0058    flag=0;
0059    if(cscnmb ==1) flag=1;
0060    }
0061    xLabel_A_ME_CSC.push_back(ss.str().c_str());
0062    flag++;
0063 }
0064 
0065 flag=4;
0066 for(Int_t i=1;i<=72;i++) {
0067   
0068    ss.str("");
0069    if(flag==5) {
0070    Int_t cscnmb;
0071    if(i<37) cscnmb=i-37;
0072    if(i>36) cscnmb=i-36;
0073    if(cscnmb !=1) ss<<cscnmb;
0074    flag=0;
0075    if(cscnmb ==1) flag=1;
0076    }
0077    xLabel_B_ME_CSC.push_back(ss.str().c_str());
0078    flag++;
0079 }
0080 
0081 // for CSC in all ME
0082 std::vector<std::string> yTitle_ME;
0083 yTitle_ME.push_back("ME- 4/2"); yTitle_ME.push_back("ME- 4/1");
0084 yTitle_ME.push_back("ME- 3/2"); yTitle_ME.push_back("ME- 3/1");
0085 yTitle_ME.push_back("ME- 2/2"); yTitle_ME.push_back("ME- 2/1");
0086 yTitle_ME.push_back("ME- 1/3"); yTitle_ME.push_back("ME- 1/2");
0087 yTitle_ME.push_back("ME- 1/1");
0088 yTitle_ME.push_back("ME+ 1/1");
0089 yTitle_ME.push_back("ME+ 1/2"); yTitle_ME.push_back("ME+ 1/3");
0090 yTitle_ME.push_back("ME+ 2/1"); yTitle_ME.push_back("ME+ 2/2");
0091 yTitle_ME.push_back("ME+ 3/1"); yTitle_ME.push_back("ME+ 3/2");
0092 yTitle_ME.push_back("ME+ 4/1"); yTitle_ME.push_back("ME+ 4/2");
0093 
0094 // input hist names
0095   std::string input_histName = "gas_gain_rechit_adc_3_3_sum_location_ME_";  
0096   std::string input_title_X="Location=(layer-1)*nsegm+segm";
0097   std::string input_title_Y="3X3 ADC Sum";
0098 
0099   std::string slice_title_X="3X3 ADC Sum Location";
0100 
0101   Int_t ny=30;
0102   Float_t ylow=1.0, yhigh=31.0;
0103   std::string result_histName = "landau_peak_position_vs_location_csc_ME_";
0104   std::string result_histTitle="Landau peak position";
0105   std::string result_title_Y="Location=(layer-1)*nsegm+segm";
0106 
0107   std::string result_histNameError = "landau_peak_position_error_vs_location_csc_ME_";
0108   std::string result_histTitleError="Landau peak position error";
0109 
0110   std::string result_histNameEntries = "entries_gas_gain_vs_location_csc_ME_";
0111   std::string result_histTitleEntries="Entries 3X3 ADC Sum";
0112 
0113   std::string result_histNameLandauPeakPositionME = "landau_peak_position_ME";
0114   std::string result_histTitleLandauPeakPositionME="Landau peak position for ME";
0115   std::string result_title_X_LandauPeakPositionME="Landau peak position";
0116   std::string result_title_Y_LandauPeakPositionME="Entries";
0117 
0118   std::string result_histNameADCSum = "adcsum_ME_";
0119   std::string result_histTitleADCSum="3X3 ADC Sum for ME";
0120   std::string result_title_X_ADCSum="3X3 ADC Sum";
0121   std::string result_title_Y_ADCSum="Entries";
0122 
0123   std::string result_histNameDeltaHV = "delta_hv";
0124   std::string result_histTitleDeltaHV="Gas Gain equalizing Delta(HV)";
0125   std::string result_title_X_DeltaHV="Delta(HV), volts";
0126   std::string result_title_Y_DeltaHV="Entries";
0127 
0128 std::string result_graphNameLandauCSCME = "graph_landau_peak_position_csc_ME_";
0129 std::string result_graphTitleLandauCSCME = "Landau peak position vs CSC ME";
0130 
0131 f_out.cd();
0132 f_out.mkdir("Input_hists");
0133 f_out.mkdir("Y_projections");
0134 f_out.mkdir("Slices");
0135 f_out.mkdir("Results");
0136 
0137 
0138 TH2F *h2;
0139 TH2F *h2mod;
0140 TH2F *h;
0141 TH2F *hlpp[9],*hlpperror[9],*hentr[9];
0142 TH1F *hlppme[18],*hadcsum[9];
0143 TH1F *hdeltahv, *widthlppratioME11,*widthlppratioME11, *statusflag;
0144 TH2F *h_min_csc_me,*h_max_csc_me, *h2mod;
0145 
0146 TGraphErrors *gr_lpp_csc_me[9];
0147 TCanvas *cgraph[9];
0148 
0149 vector<Float_t> gr_y[9],gr_y_er[9],gr_x[9];
0150 
0151 Int_t esr[18]={111,112,113,121,122,131,132,141,142,
0152                211,212,213,221,222,231,232,241,242};
0153 Int_t entries[18]={0,0,0,0,0,0,0,0,0,
0154                    0,0,0,0,0,0,0,0,0};
0155 
0156 // station, ring sequence, # of CSC and HV segments
0157 Int_t sr[9]={11,12,13,21,22,31,32,41,42};
0158 Int_t ncsc[9]={36,36,36,18,36,18,36,18,36};
0159 Int_t nhvsegm[9]={6,18,18,18,30,18,30,18,30};   // number of HV segments per CSC              
0160 Float_t dhv_coeff[9]={158.0,190.0,190.0,190.0,190.0,190.0,190.0,190.0,190.0};
0161 Int_t nentr[9]={0};
0162 Int_t ncscmis=0;
0163 Int_t indsegm=0;
0164 vector<Int_t> cscmis; cscmis.push_back(0);
0165 
0166 vector<Int_t> selflag,locid, nentrloc;
0167 selflag.clear(); locid.clear(); nentrloc.clear();
0168 
0169 vector<Float_t> peak, peakpos, peakwidth, peakposer, dhvv, meanadclpp;
0170 peak.clear(); peakpos.clear(); peakwidth.clear(); peakposer.clear(); 
0171 dhvv.clear(); meanadclpp.clear();
0172 
0173 Int_t locid_cur, nentrloc_cur, selflag_cur;
0174 Float_t peak_cur, peakpos_cur, peakwidth_cur, peakposer_cur, dhv_cur;
0175 
0176 
0177 /// book output histograms with all CSCs for given ME+- station,ring
0178 
0179 for(Int_t isr=0;isr<8;isr++) {
0180   std::vector<std::string> xLabel_ME_CSC; xLabel_ME_CSC.clear();
0181      Int_t ny=nhvsegm[isr];
0182      Float_t ylow=1.0; Float_t yhigh=ny; yhigh=yhigh+1.0;
0183      if(ncsc[isr]==36) {
0184        Int_t nx=72; Float_t xlow=-35.0; Float_t xhigh=37.0;
0185        xLabel_ME_CSC=xLabel_B_ME_CSC;
0186      }
0187      if(ncsc[isr]==18) {
0188        Int_t nx=36; Float_t xlow=-17.0; Float_t xhigh=19.0;
0189        xLabel_ME_CSC=xLabel_A_ME_CSC;
0190      }
0191 
0192   // 2D hists for Landau peak in 3X3 ADC Sum and its error vs hv segment and 
0193   // csc in ME+-
0194      
0195      ss.str("");
0196      ss<<result_histName.c_str()<<sr[isr];
0197      ss1.str("");
0198      ss1<<result_histTitle<<" in run "<<run.c_str();
0199      hlpp[isr]=new TH2F(ss.str().c_str(),ss1.str().c_str(),nx,xlow,xhigh,ny,ylow,yhigh);
0200      for(Int_t i=1;i<=nx;i++) hlpp[isr]->GetXaxis()->SetBinLabel(i,xLabel_ME_CSC[i-1].c_str());
0201      hlpp[isr]->SetStats(kFALSE);
0202      hlpp[isr]->GetXaxis()->SetTitle(xTitle_ME_CSC[isr].c_str());
0203      hlpp[isr]->GetYaxis()->SetTitle(result_title_Y.c_str());
0204      hlpp[isr]->GetZaxis()->SetLabelSize(0.03);
0205      hlpp[isr]->SetOption("COLZ");
0206      hlpp[isr]->SetMinimum(0.0);
0207      hlpp[isr]->SetMaximum(2000.0);
0208 
0209      ss.str("");
0210      ss<<result_histNameError.c_str()<<sr[isr];
0211      ss1.str("");
0212      ss1<<result_histTitleError<<" in run "<<run.c_str();
0213      hlpperror[isr]=new TH2F(ss.str().c_str(),ss1.str().c_str(),nx,xlow,xhigh,ny,ylow,yhigh);
0214      for(Int_t i=1;i<=nx;i++) hlpperror[isr]->GetXaxis()->SetBinLabel(i,xLabel_ME_CSC[i-1].c_str());
0215      hlpperror[isr]->SetStats(kFALSE);
0216      hlpperror[isr]->GetXaxis()->SetTitle(xTitle_ME_CSC[isr].c_str());
0217      hlpperror[isr]->GetYaxis()->SetTitle(result_title_Y.c_str());
0218      hlpperror[isr]->GetZaxis()->SetLabelSize(0.03);
0219      hlpperror[isr]->SetOption("COLZ");
0220      hlpperror[isr]->SetMinimum(0.0);
0221      hlpperror[isr]->SetMaximum(100.0);
0222 
0223 
0224   // 2D hists for entries 3X3 ADC Sum vs hv segment and csc in given ME+-
0225      ss.str("");
0226      ss<<result_histNameEntries.c_str()<<sr[isr];
0227      ss1.str("");
0228      ss1<<result_histTitleEntries<<" in run "<<run.c_str();
0229      hentr[isr]=new TH2F(ss.str().c_str(),ss1.str().c_str(),nx,xlow,xhigh,ny,ylow,yhigh);
0230      for(Int_t i=1;i<=nx;i++) hentr[isr]->GetXaxis()->SetBinLabel(i,xLabel_ME_CSC[i-1].c_str());
0231      hentr[isr]->SetStats(kFALSE);
0232      hentr[isr]->GetXaxis()->SetTitle(xTitle_ME_CSC[isr].c_str());
0233      hentr[isr]->GetYaxis()->SetTitle(result_title_Y.c_str());
0234      hentr[isr]->GetZaxis()->SetLabelSize(0.03);
0235      hentr[isr]->SetOption("COLZ");
0236 
0237      // 1D hists for 3X3 ADC Sum distribution in given ME 
0238      ss.str("");
0239      ss<< result_histNameADCSum <<sr[isr];
0240      ss1.str("");
0241      ss1<< result_histTitleADCSum<<sr[isr] <<" in run "<<run.c_str();
0242      hadcsum[isr]=new TH1F(ss.str().c_str(),ss1.str().c_str(),50,0.0,2000.0);
0243      
0244      hadcsum[isr]->GetXaxis()->SetTitle(result_title_X_ADCSum.c_str());
0245      hadcsum[isr]->GetYaxis()->SetTitle(result_title_Y_ADCSum.c_str());
0246      hadcsum[isr]->SetFillColor(4);
0247 
0248 } // end of for(Int_t isr=0
0249 
0250   
0251   // 1D hists for gas gain as Landau peak position per csc in given ME 
0252 for(Int_t jesr=0;jesr<18;jesr++) {
0253   if(esr[jesr] != 142 && esr[jesr] != 242) {
0254     if(esr[jesr] < 200) {Int_t isr=jesr;   std::string endcapsign="+";}
0255     if(esr[jesr] > 200) {Int_t isr=jesr-9; std::string endcapsign="-";}
0256      ss.str("");
0257      ss<< result_histNameLandauPeakPositionME <<endcapsign.c_str()<<sr[isr];
0258      ss1.str("");
0259      ss1<< result_histTitleLandauPeakPositionME<<endcapsign.c_str()<<sr[isr] <<" in run "<<run.c_str();
0260      hlppme[jesr]=new TH1F(ss.str().c_str(),ss1.str().c_str(),80,0.0,1600.0);
0261      
0262      hlppme[jesr]->GetXaxis()->SetTitle(result_title_X_LandauPeakPositionME.c_str());
0263      hlppme[jesr]->GetYaxis()->SetTitle(result_title_Y_LandauPeakPositionME.c_str());
0264      hlppme[jesr]->SetFillColor(4);
0265   }
0266 }
0267 // book two output 2D hists for min. and max. gas gain vs CSC and ME
0268 
0269   ss.str("");
0270   ss<<"min_landau_peak_position_csc_ME";
0271   ss1.str("");
0272   ss1<<"Min. Landau peak position vs CSC and ME"<<" in run "<<run.c_str();
0273   h_min_csc_me=new TH2F(ss.str().c_str(),ss1.str().c_str(),36,1.0,37.0,18,1.0,19.0);
0274   h_min_csc_me->SetStats(kFALSE);
0275   h_min_csc_me->GetXaxis()->SetTitle("CSC #");
0276   for(Int_t i=1;i<=18;i++) h_min_csc_me->GetYaxis()->SetBinLabel(i,yTitle_ME[i-1].c_str());
0277   h_min_csc_me->GetZaxis()->SetLabelSize(0.03);
0278   h_min_csc_me->SetOption("COLZ");
0279   h_min_csc_me->SetMinimum(0.0);
0280   h_min_csc_me->SetMaximum(2000.0);
0281  
0282   ss.str("");
0283   ss<<"max_landau_peak_position_csc_ME";
0284   ss1.str("");
0285   ss1<<"Max. Landau peak position vs CSC and ME"<<" in run "<<run.c_str();
0286   h_max_csc_me=new TH2F(ss.str().c_str(),ss1.str().c_str(),36,1.0,37.0,18,1.0,19.0);
0287   h_max_csc_me->SetStats(kFALSE);
0288   h_max_csc_me->GetXaxis()->SetTitle("CSC #");
0289   for(Int_t i=1;i<=18;i++) h_max_csc_me->GetYaxis()->SetBinLabel(i,yTitle_ME[i-1].c_str());
0290   h_max_csc_me->GetZaxis()->SetLabelSize(0.03);
0291   h_max_csc_me->SetOption("COLZ");
0292   h_max_csc_me->SetMinimum(0.0);
0293   h_max_csc_me->SetMaximum(2000.0);
0294 
0295   
0296   // 1D hist for delta hv  
0297   ss.str("");
0298   ss<< result_histNameDeltaHV;
0299   ss1.str("");
0300   ss1<< result_histTitleDeltaHV<<" in run "<<run.c_str();
0301   hdeltahv=new TH1F(ss.str().c_str(),ss1.str().c_str(),80,-200.0,200.0);
0302   hdeltahv->GetXaxis()->SetTitle(result_title_X_DeltaHV.c_str());
0303   hdeltahv->GetYaxis()->SetTitle(result_title_Y_DeltaHV.c_str());
0304   hdeltahv->SetFillColor(4);
0305 
0306   // 1D hist for ratio in Landau fit for ME11
0307   ss.str("");
0308   ss<< "widthlppratioME11";
0309   ss1.str("");
0310   ss1<< "Width to Landau Peak Ratio ME11 "<<" in run "<<run.c_str();
0311   widthlppratioME11=new TH1F(ss.str().c_str(),ss1.str().c_str(),100,0.0,1.0);
0312   widthlppratioME11->GetXaxis()->SetTitle("Width to Landau Peak Ratio ME11");
0313   widthlppratioME11->GetYaxis()->SetTitle("Entries");
0314   widthlppratioME11->SetFillColor(4);  
0315 
0316   // 1D hist for ratio in Landau fit for the rest of ME
0317   ss.str("");
0318   ss<< "widthlppratioME";
0319   ss1.str("");
0320   ss1<< "Width to Landau Peak Ratio ME "<<" in run "<<run.c_str();
0321   widthlppratioME=new TH1F(ss.str().c_str(),ss1.str().c_str(),100,0.0,1.0);
0322   widthlppratioME->GetXaxis()->SetTitle("Width to Landau Peak Ratio ME");
0323   widthlppratioME->GetYaxis()->SetTitle("Entries");
0324   widthlppratioME->SetFillColor(4);  
0325 
0326   // 1D hist for status flag
0327   ss.str("");
0328   ss<< "StatusFlag";
0329   ss1.str("");
0330   ss1<< "Status Flag "<<" in run "<<run.c_str();
0331   statusflag=new TH1F(ss.str().c_str(),ss1.str().c_str(),10,0.0,10.0);
0332   statusflag->GetXaxis()->SetTitle("Status Flag");
0333   statusflag->GetYaxis()->SetTitle("Entries");
0334   statusflag->SetFillColor(4); 
0335 
0336 Int_t k=0;
0337 TCanvas *c1=new TCanvas("c1","canvas");
0338 c1->cd();
0339 
0340 for(Int_t jesr=0;jesr<18;jesr++) { 
0341 
0342   if(esr[jesr] != 142 && esr[jesr] != 242) {
0343   Int_t me;
0344   if(jesr<9) me=10+jesr;
0345   if(jesr>8) me=18-jesr;
0346 
0347   Int_t indisr;
0348   if(esr[jesr] < 200) indisr=esr[jesr]-100;
0349   if(esr[jesr] > 200) indisr=esr[jesr]-200;
0350 
0351   // define index for vector sr[9] (station,ring)
0352   for(Int_t i=0;i<8;i++) if(sr[i]==indisr) indisr=i;
0353  
0354   // cycle over CSC 
0355   for(Int_t csc=1;csc<(ncsc[indisr]+1);csc++) {
0356 
0357      Int_t cscbin;
0358      // define bin cscbin in hists for ME+-
0359      if(esr[jesr]<200) cscbin=ncsc[indisr]+csc;
0360      if(esr[jesr]>200) cscbin=ncsc[indisr]-csc+1;
0361 
0362      Int_t idchamber=esr[jesr]*100+csc;
0363      ss.str("");
0364      ss<<in_folder.c_str()<<input_histName.c_str()<<idchamber;
0365      f_in.cd();
0366      TH2F *h2 = (TH2F*)f_in.Get(ss.str().c_str());
0367      if(h2 ==NULL) {
0368        ncscmis++;
0369        cscmis.push_back(idchamber);
0370      }
0371    if(h2 != NULL) {
0372 
0373      ss.str(""); ss1.str("");
0374      ss<<input_histName.c_str()<<idchamber;
0375      ss1<<h2->GetTitle()<<" in run "<<run.c_str();
0376      ny=h2->GetNbinsY(); Float_t ylow=h2->GetYaxis()->GetBinLowEdge(1);
0377                          Float_t yhigh=h2->GetYaxis()->GetBinLowEdge(ny)+
0378                                        h2->GetYaxis()->GetBinWidth(ny);
0379      Int_t nx=nhvsegm[indisr]; Float_t xlow=1.0; Float_t xhigh=xlow+nx;
0380      h2mod=new TH2F(ss.str().c_str(),ss1.str().c_str(),nx,xlow,xhigh,ny,ylow,yhigh);
0381      h2mod->SetStats(kFALSE);
0382 
0383      // saving original, adding X,Y titles, color and "BOX" option
0384      Float_t w;
0385      for(Int_t i=1;i<=h2->GetNbinsX();i++) { 
0386        if(i <= nhvsegm[indisr]) {
0387          for(Int_t j=1;j<=h2->GetNbinsY();j++) {
0388         w=h2->GetBinContent(i,j);
0389         h2mod->SetBinContent(i,j,w);
0390             Float_t x=h2->GetYaxis()->GetBinLowEdge(j);
0391             nentr[indisr]=nentr[indisr]+1;
0392             hadcsum[indisr]->Fill(x,w);
0393      }
0394        }     
0395      }
0396      h2mod->GetXaxis()->SetTitle(input_title_X.c_str());
0397      h2mod->GetYaxis()->SetTitle(input_title_Y.c_str());
0398      h2mod->GetYaxis()->SetTitleOffset(1.2);
0399      h2mod->SetOption("COLZ");
0400      f_out.cd("Input_hists");     
0401      h2mod->Write();
0402 
0403      // saving Y projection of the whole 2D hist for given chamber
0404      ss.str("");
0405      ss<<input_histName.c_str()<<idchamber<<"_Y_all";
0406      TH1D *h1d = h2->ProjectionY(ss.str().c_str(),1,nhvsegm[indisr],"");
0407      h1d->GetXaxis()->SetTitle(input_title_Y.c_str());
0408      h1d->GetYaxis()->SetTitle("Entries");
0409      h1d->GetYaxis()->SetTitleOffset(1.2);
0410      h1d->SetFillColor(4);
0411      gStyle->SetOptStat(1001111);
0412      f_out.cd("Y_projections");
0413 
0414         // Fit by Landau peak the 3X3 ADC Sum distribution in the whole csc
0415      TF1 *fcsc=new TF1("fcsc","landau",0.0,2000.0);
0416      Float_t inpeak=5.0*h1d->GetMaximum(99999.0);
0417      Float_t inpeakpos=0.6*h1d->GetMean();
0418      Float_t inpeakwidth=0.3*inpeakpos;
0419      fcsc->SetParameters(inpeak,inpeakpos,inpeakwidth);
0420      //h1d->SetNdivisions(506);
0421      h1d->Fit("fcsc","BQ","",0.0,2000.0);
0422      Float_t lppcsc=fcsc->GetParameter(1);
0423      Float_t lppcscerr=fcsc->GetParError(1);
0424 
0425      h1d->Write(); 
0426 
0427      gr_y[indisr].push_back(lppcsc);
0428      gr_y_er[indisr].push_back(lppcscerr);
0429      Float_t xcsc=csc;
0430      if(esr[jesr] > 200) xcsc=-csc;
0431      gr_x[indisr].push_back(xcsc);
0432 
0433      delete h1d;   
0434 
0435      // saving slices, fitting Landau peak in each slice, fill 2D hist
0436      f_out.cd("Slices");
0437      Float_t mingasgain=10000.0;
0438      Float_t maxgasgain=0.0;
0439      for(Int_t j=1;j<=h2->GetNbinsX();j++) {
0440        if(j <= nhvsegm[indisr]) {
0441         Int_t n=j;
0442     locid_cur=idchamber*100+n;
0443         selflag_cur=0;
0444     peak_cur=0.0; peakpos_cur=0.0; peakwidth_cur=0.0; peakposer_cur=0.0;
0445         dhv_cur=0.0;
0446 
0447         ss.str("");
0448         ss<<input_histName.c_str()<<idchamber<<"_Y_"<<n;
0449         TH1D *h1d = h2->ProjectionY(ss.str().c_str(),j,j,"");
0450         nentrloc_cur=h1d->GetEntries();
0451         if( nentrloc_cur == 0) selflag_cur=1;
0452         if(nentrloc_cur > 0 && nentrloc_cur <= nminentries) selflag_cur=2;
0453 
0454     //  if(h1d->GetEntries() > 0 ) {
0455         if(nentrloc_cur > nminentries) {   // since Sep. 29, 2008
0456           /* Fit by Landau peak the 3X3 ADC Sum distribution in the slice*/
0457           TF1 *f1=new TF1("f1","landau",0.0,2000.0);
0458           Float_t inpeak=5.0*h1d->GetMaximum(99999.0);
0459           Float_t inpeakpos=0.6*h1d->GetMean();
0460           Float_t inpeakwidth=0.3*inpeakpos;
0461           f1->SetParameters(inpeak,inpeakpos,inpeakwidth);
0462           //h1d->SetNdivisions(506);
0463           h1d->Fit("f1","BQ","",0.0,2000.0);
0464           peak_cur=f1->GetParameter(0);
0465           peakpos_cur= f1->GetParameter(1);
0466           peakwidth_cur=f1->GetParameter(2);
0467           peakposer_cur=f1->GetParError(1);
0468 
0469           Float_t lpp=f1->GetParameter(1);
0470       Float_t lpperr=f1->GetParError(1);
0471 
0472           Float_t ratio=0.0;
0473           if(peakpos_cur != 0.0) ratio=peakwidth_cur/peakpos_cur;
0474           if(sr[indisr] == 11) widthlppratioME11->Fill(ratio,1.0);
0475           if(sr[indisr] != 11) widthlppratioME->Fill(ratio,1.0);
0476 
0477           if(peakpos_cur < 0.0 || peakpos_cur > 2000.0 || 
0478         ratio <= 0.0 || ratio>=1.0) selflag_cur=selflag_cur+4; // Fit failed
0479 
0480           if(lpperr>100.0) lpperr=100.0;
0481           Float_t entr=h1d->GetEntries();
0482           entries[jesr]=entries[jesr]+1;
0483           hlpp[indisr]->SetBinContent(cscbin,j,lpp);
0484           hlpperror[indisr]->SetBinContent(cscbin,j,lpperr);
0485           hentr[indisr]->SetBinContent(cscbin,j,entr);
0486           if(selflag_cur == 0 ) {
0487             hlppme[jesr]->Fill(lpp,1.0);
0488             if(lpp < mingasgain) mingasgain=lpp;
0489             if(lpp > maxgasgain) maxgasgain=lpp;
0490       }
0491           ss.str("");
0492           ss<<slice_title_X<<" "<<n;
0493           h1d->GetXaxis()->SetTitle(ss.str().c_str());
0494           h1d->GetYaxis()->SetTitle("Entries");
0495           h1d->GetYaxis()->SetTitleOffset(1.2);
0496           h1d->SetFillColor(4);
0497           gStyle->SetOptStat(1001111);
0498           h1d->Write();
0499 
0500     } // end of if(h1d->GetEntries() > 0)
0501         delete h1d;
0502         Float_t dum=selflag_cur;
0503         statusflag->Fill(dum,1.0);
0504 
0505         locid.push_back(locid_cur);
0506         nentrloc.push_back(nentrloc_cur);
0507         selflag.push_back(selflag_cur);
0508         peak.push_back(peak_cur);
0509         peakpos.push_back(peakpos_cur);
0510         peakwidth.push_back(peakwidth_cur);
0511         peakposer.push_back(peakposer_cur);
0512         dhvv.push_back(dhv_cur);  // dhv will be modified later
0513 
0514        }  // end of if(j <= nhvsegm[indisr]) {
0515      }    // end of for(Int_t j=1;j<=h2->GetNbinsX();j++) {
0516 
0517    if(mingasgain < 10000.0) h_min_csc_me->SetBinContent(csc,me,mingasgain);
0518    if(maxgasgain > 0.0) h_max_csc_me->SetBinContent(csc,me,maxgasgain);
0519    delete h2;
0520    delete h2mod; 
0521    } // end of if(h2 != NULL)
0522   }  //  for(Int_t csc=1;csc<(ncsc[indisr]+1);csc++) {
0523   }  // end of   if(esr[jesr] != 142 && esr[jesr] != 242) {
0524 }    // end of for(Int_t jesr=0;jesr<18;jesr++) { 
0525 
0526 
0527 f_out.cd("Results");
0528 for(Int_t isr=0;isr<8;isr++) {
0529   if(hlpp[isr] != NULL) {
0530        hlpp[isr]->SetStats(kFALSE);
0531        hlpp[isr]->Write();
0532   }
0533   if(hlpperror[isr] != NULL) {
0534        hlpperror[isr]->SetStats(kFALSE);
0535        hlpperror[isr]->Write();
0536   }
0537   if(hentr[isr] != NULL) {
0538        hentr[isr]->SetStats(kFALSE);
0539        hentr[isr]->Write();
0540   }
0541   if(hadcsum[isr] != NULL) {
0542        hadcsum[isr]->Write();
0543   }
0544 }
0545 
0546 for(Int_t jesr=0;jesr<18;jesr++) if(hlppme[jesr] != NULL) 
0547     hlppme[jesr]->Write();
0548 
0549 h_min_csc_me->Write();
0550 h_max_csc_me->Write();
0551 widthlppratioME11->Write();
0552 widthlppratioME->Write();
0553 statusflag->Write();
0554 
0555 /// Report missing CSCs
0556 std::cout<<"Missing csc "<<ncscmis<<std::endl;
0557 if(ncscmis > 0) {
0558   for(Int_t jesr=0;jesr<18;jesr++) {
0559     Int_t n=0,npr=0;
0560     for(Int_t i=0;i<cscmis.size();i++) {
0561       Int_t idesr=cscmis[i]/100;
0562       if(esr[jesr]==idesr && n<10) {n++;npr++;std::cout<<cscmis[i]<<"   ";}
0563       if(n==10) {std::cout<<std::endl; n=0;} 
0564     }
0565     if(npr>0) std::cout<<std::endl;
0566   }
0567 }
0568 
0569 // Print table
0570 
0571 std::cout<<std::endl;
0572 std::cout<<"Tables for CSC HV adjustment"<<std::endl;
0573 std::cout<<std::endl;
0574 
0575 Int_t station,ring,minentries=1000000,maxentries=0;
0576 for(Int_t jesr=0;jesr<18;jesr++) {
0577   if(hlppme[jesr] != NULL) { 
0578   if(esr[jesr] < 200) {Int_t isr=jesr;   std::string endcapsign="+";}
0579   if(esr[jesr] > 200) {Int_t isr=jesr-9; std::string endcapsign="-";}
0580   station=sr[isr]/10; ring=sr[isr]-station*10;
0581   ss.str("");
0582   ss<<"ME"<<endcapsign.c_str()<<station<<"/"<<ring<< 
0583   "      MEAN(Landau peak position) "<<hlppme[jesr]->GetMean();
0584   std::cout<<ss.str().c_str()<<std::endl;
0585   }
0586 }
0587 std::cout<<std::endl;
0588 
0589 std::cout<<"Total hv segments analyzed "<<locid.size()<<std::endl;
0590 std::cout<<std::endl;
0591 
0592 for(Int_t ind=0; ind<locid.size(); ind++) {
0593   Int_t jesr=locid[ind]/10000;
0594   for(Int_t i=0;i<18;i++) if(esr[i] == jesr) jesr=i;
0595   Float_t mean_adc_lpp=hlppme[jesr]->GetMean();
0596   meanadclpp.push_back(mean_adc_lpp);
0597 
0598   Int_t isr=locid[ind]/1000000;
0599   isr=locid[ind]-isr*1000000;
0600   isr=isr/10000;
0601   for(Int_t i=0;i<8;i++) if(sr[i] == isr) isr=i;
0602 
0603   if(selflag[ind] == 0) {   // statistics is acceptable and fit was OK 
0604     Float_t dhvc=dhv_coeff[isr];
0605     Float_t adc_lpp=peakpos[ind];
0606     Float_t dhv=-dhvc*log(adc_lpp/mean_adc_lpp);
0607     dhvv[ind]=dhv;
0608     hdeltahv->Fill(dhv,1.0);
0609   }
0610   printf("%7i%7.0f%3i%5i%6i%8.0f%8.0f%8.0f%8.0f%8.0f\n", locid[ind],meanadclpp[ind],selflag[ind],nminentries,nentrloc[ind],peak[ind],peakpos[ind],peakwidth[ind],peakposer[ind],dhvv[ind]);
0611 }
0612 
0613 f_out.cd("Results");
0614 hdeltahv->Write();
0615 
0616 /// Plot and save graphs of the mean Landau peak position vs CSC in ME
0617 Float_t xgr[72],ygr[72],erygr[72],erxgr[72];
0618 for(Int_t isr=0;isr<8;isr++) {
0619   if(gr_x[isr].size()>0) {
0620     Int_t np=gr_x[isr].size();
0621     for(Int_t i=0;i<np; i++) {
0622       xgr[i]=gr_x[isr][i];
0623       ygr[i]=gr_y[isr][i];
0624       erxgr[i]=0.0;
0625       erygr[i]=gr_y_er[isr][i];
0626     }
0627     ss.str("");
0628     ss<<result_graphNameLandauCSCME.c_str()<<sr[isr];
0629 
0630     cgraph[isr]= new TCanvas(ss.str().c_str());
0631     cgraph[isr]->cd();
0632     cgraph[isr]->SetGrid();
0633     gr_lpp_csc_me[isr]=new TGraphErrors(np,xgr,ygr,erxgr,erygr);
0634     ss1.str("");
0635     ss1<<result_graphTitleLandauCSCME<<sr[isr]<<" in run "<<run.c_str();
0636     gr_lpp_csc_me[isr]->SetTitle(ss1.str().c_str());
0637     gr_lpp_csc_me[isr]->GetXaxis()->SetTitle(xTitle_ME_CSC[isr].c_str());
0638     gr_lpp_csc_me[isr]->GetYaxis()->SetTitle("Landau peak position");
0639     gr_lpp_csc_me[isr]->SetMinimum(0.0);
0640     gr_lpp_csc_me[isr]->SetMaximum(500.0);
0641     if(isr==0) {
0642       gr_lpp_csc_me[isr]->SetMinimum(250.0);
0643       gr_lpp_csc_me[isr]->SetMaximum(1250.0);
0644     }
0645     gr_lpp_csc_me[isr]->SetMarkerStyle(20);
0646     gr_lpp_csc_me[isr]->SetMarkerColor(4);
0647     gr_lpp_csc_me[isr]->SetMarkerSize(1.2);
0648     gr_lpp_csc_me[isr]->SetLineColor(4);    // Blue error bar
0649     gr_lpp_csc_me[isr]->Draw("APZ"); 
0650     cgraph[isr]->Update();
0651     f_out.cd("Results");
0652     cgraph[isr]->Write();
0653     delete gr_lpp_csc_me[isr];
0654     delete cgraph[isr];
0655   }
0656 }
0657 //f_out.close();
0658 
0659 FILE *File_out_recom;
0660 ss.str("");
0661 ss<<"dhv_recom_fit_"<<run<<".txt";
0662 File_out_recom=fopen(ss.str().c_str(),"w");
0663 
0664 for(Int_t ind=0; ind<locid.size(); ind++) {
0665   if(selflag[ind] == 0) {
0666     Int_t d=dhvv[ind];
0667     if(d > 100 ) d=100;
0668     if((d >= -15) && (d <= 15)) d=0; 
0669     if(d < -15 || d > 15) {
0670       Float_t df;
0671       if(d<0) {df=d-5; d=df/10; d=d*10;}
0672       if(d>0) {df=d+5; d=df/10; d=d*10;}
0673     }    
0674   fprintf(File_out_recom,"%7i%7.0f%3i%5i%6i%8.0f%8.0f%8.0f%8.0f%8.0f%8i\n", locid[ind],meanadclpp[ind],selflag[ind],nminentries,nentrloc[ind],peak[ind],peakpos[ind],peakwidth[ind],peakposer[ind],dhvv[ind],d);
0675   }
0676 }
0677 fclose(File_out_recom);
0678 
0679 FILE *File_out_largedhvpos;
0680 ss.str("");
0681 ss<<"hv_segm_large_dhv_pos_"<<run<<".txt";
0682 File_out_largedhvpos=fopen(ss.str().c_str(),"w");
0683 
0684 for(Int_t ind=0; ind<locid.size(); ind++) {
0685   if(selflag[ind] == 0) {
0686     Int_t d=dhvv[ind];
0687     if(d > 100) {
0688       d=100;
0689       fprintf(File_out_largedhvpos,"%7i%7.0f%3i%5i%6i%8.0f%8.0f%8.0f%8.0f%8.0f%8i\n", locid[ind],meanadclpp[ind],selflag[ind],nminentries,nentrloc[ind],peak[ind],peakpos[ind],peakwidth[ind],peakposer[ind],dhvv[ind],d);
0690     }
0691   }
0692 }
0693 fclose(File_out_largedhvpos);
0694 
0695 FILE *File_out_largedhvneg;
0696 ss.str("");
0697 ss<<"hv_segm_large_dhv_neg_"<<run<<".txt";
0698 File_out_largedhvneg=fopen(ss.str().c_str(),"w");
0699 
0700 for(Int_t ind=0; ind<locid.size(); ind++) {
0701   if(selflag[ind] == 0) {
0702     Int_t d=dhvv[ind];
0703     if(d <-100) {
0704       Float_t df;
0705       df=d-5; d=df/10; d=d*10;  
0706       fprintf(File_out_largedhvneg,"%7i%7.0f%3i%5i%6i%8.0f%8.0f%8.0f%8.0f%8.0f%8i\n", locid[ind],meanadclpp[ind],selflag[ind],nminentries,nentrloc[ind],peak[ind],peakpos[ind],peakwidth[ind],peakposer[ind],dhvv[ind],d);
0707     }
0708   }
0709 }
0710 fclose(File_out_largedhvneg);
0711 
0712 
0713 FILE *File_out_nodata;
0714 ss.str("");
0715 ss<<"hv_segm_no_data_"<<run<<".txt";
0716 File_out_nodata=fopen(ss.str().c_str(),"w");
0717 
0718 for(Int_t ind=0; ind<locid.size(); ind++) {
0719   if(selflag[ind] == 1) 
0720   fprintf(File_out_nodata,"%7i%7.0f%3i%5i%6i%8.0f%8.0f%8.0f%8.0f%8i\n", locid[ind],meanadclpp[ind],selflag[ind],nminentries,nentrloc[ind],peak[ind],peakpos[ind],peakwidth[ind],peakposer[ind],dhvv[ind]);
0721   
0722 }
0723 fclose(File_out_nodata);
0724 
0725 FILE *File_out_lowstat;
0726 ss.str("");
0727 ss<<"hv_segm_low_stat_"<<run<<".txt";
0728 File_out_lowstat=fopen(ss.str().c_str(),"w");
0729 
0730 for(Int_t ind=0; ind<locid.size(); ind++) {
0731   if((selflag[ind] & 2) > 0) 
0732   fprintf(File_out_lowstat,"%7i%7.0f%3i%5i%6i%8.0f%8.0f%8.0f%8.0f%8i\n", locid[ind],meanadclpp[ind],selflag[ind],nminentries,nentrloc[ind],peak[ind],peakpos[ind],peakwidth[ind],peakposer[ind],dhvv[ind]);
0733   
0734 }
0735 fclose(File_out_lowstat);
0736 
0737 FILE *File_out_failedfit;
0738 ss.str("");
0739 ss<<"hv_segm_failed_fit_"<<run<<".txt";
0740 File_out_failedfit=fopen(ss.str().c_str(),"w");
0741 
0742 for(Int_t ind=0; ind<locid.size(); ind++) {
0743   if((selflag[ind] & 4) > 0) 
0744   fprintf(File_out_failedfit,"%7i%7.0f%3i%5i%6i%8.0f%8.0f%8.0f%8.0f%8i\n", locid[ind],meanadclpp[ind],selflag[ind],nminentries,nentrloc[ind],peak[ind],peakpos[ind],peakwidth[ind],peakposer[ind],dhvv[ind]); 
0745 }
0746 fclose(File_out_failedfit);
0747 }