CMSSW Reference Manual

00001 // Analysis of HF LED/Laser run: 
00002 // SPE calibration for low light intensity or raw SPE calibration for high light intensity
00003 // and HF performance based on this analysis
00004 //
00005 // Igor Vodopiyanov. Oct-2007
00006 //
00007 #include <memory>
00008 #include <string>
00009 #include <iostream>
00010 
00011 #include "TH1F.h"
00012 #include "TH2F.h"
00013 #include "TFile.h"
00014 #include "math.h"
00015 #include "TMath.h"
00016 #include "TF1.h"
00017 
00018 #include "CalibCalorimetry/HcalStandardModules/interface/HFLightCal.h"
00019 
00020 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00021 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00022 
00023 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
00024 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
00025 #include "CondFormats/HcalObjects/interface/HcalQIEShape.h"
00026 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
00027 
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00030 #include "FWCore/Framework/interface/EventSetup.h"
00031 #include "FWCore/Framework/interface/ESHandle.h"
00032 
00033 using namespace std;
00034 Int_t Nev, runN=0, eventN=0;
00035 Int_t itsmax[26][36][2];
00036 Int_t itspinmax[8][3];
00037 
00038 namespace {
00039   //bool verbose = true;
00040   bool verbose = false;
00041 }
00042 
00043 HFLightCal::HFLightCal (const edm::ParameterSet& fConfiguration) {
00044   //std::string histfile = fConfiguration.getUntrackedParameter<string>("rootFile");
00045   histfile = fConfiguration.getUntrackedParameter<string>("rootFile");
00046   textfile = fConfiguration.getUntrackedParameter<string>("textFile");
00047   prefile = fConfiguration.getUntrackedParameter<string>("preFile");
00048 }
00049 
00050 HFLightCal::~HFLightCal () {
00051   //delete mFile;
00052 }
00053 
00054 void HFLightCal::beginJob(const edm::EventSetup& fSetup) {
00055 
00056   char htit[64];
00057   Int_t neta,nphi,ndepth,nmax,nquad,npin;
00058 
00059   std::cout<<std::endl<<"HFLightCal beginJob: --> "<<std::endl;
00060 
00061   // Read info about signal timing in TS from PreAnalysis
00062   mFile = new TFile (histfile.c_str(),"RECREATE");
00063   if ((tFile = fopen(textfile.c_str(),"w"))==NULL) {
00064     printf("\nNo output textfile open\n\n");
00065     std::cout<<"Problem with output textFILE => exit"<<std::endl;
00066     exit(1);
00067   }
00068   //if ((preFile = fopen("hf_preanal.txt","r"))==NULL){
00069   if ((preFile = fopen(prefile.c_str(),"r"))==NULL){
00070     printf("\nNo input pre-file open\n\n");
00071     std::cout<<"Problem with input textFILE => exit"<<std::endl;
00072     exit(1);
00073   }
00074   rewind(preFile);
00075   for (int i=0; i<1728; i++) {
00076     fscanf(preFile,"%d%d%d%d\r",&neta,&nphi,&ndepth,&nmax);
00077     //std::cout<<neta<<"  "<<nphi<<"  "<<ndepth<<"  "<<nmax<<std::endl;
00078     if (neta>=29 && neta<=41 && nphi<72 && nphi>0 && ndepth>0 && ndepth<=2) 
00079       itsmax[neta-29][(nphi-1)/2][ndepth-1] = nmax;
00080     else if (neta<=-29 && neta>=-41 && nphi<72 && nphi>0 && ndepth>0 && ndepth<=2) 
00081       itsmax[13-neta-29][(nphi-1)/2][ndepth-1] = nmax;
00082     else {
00083       std::cout<<"Input pre-file: wrong channel record:"<<std::endl;
00084       std::cout<<"eta="<<neta<<"  phi="<<nphi<<"  depth="<<ndepth<<"  max="<<nmax<<std::endl;
00085     }
00086   }
00087   for (int i=0; i<24; i++) {
00088     fscanf(preFile,"%d%d%d\r",&nquad,&npin,&nmax);
00089     //std::cout<<nquad<<"  "<<npin<<"  "<<nmax<<std::endl;
00090     if (nquad>0 && nquad<=4 && npin<=3 && npin>0) 
00091       itspinmax[nquad-1][npin-1] = nmax;
00092     else if (nquad<0 && nquad>=-4 && npin<=3 && npin>0) 
00093       itspinmax[4-nquad-1][npin-1] = nmax;
00094     else {
00095       std::cout<<"Input pre-file: wrong PIN record:"<<std::endl;
00096       std::cout<<"quad="<<nquad<<"  pin="<<npin<<"  max="<<nmax<<std::endl;
00097     }
00098   }
00099 
00100   // General Histos
00101   htmax = new TH1F("htmax","Max TS",10,-0.5,9.5);
00102   htmean = new TH1F("htmean","Mean signal TS",100,0,10);
00103   hsignalmean = new TH1F("hsignalmean","Mean ADC 4maxTS",1201,-25,30000);
00104   hsignalrms = new TH1F("hsignalrms","RMS ADC 4maxTS",500,0,500);
00105   hpedmean = new TH1F("hpedmean","Mean ADC 4lowTS",200,-10,90);
00106   hpedrms = new TH1F("hpedrms","RMS ADC 4lowTS",200,0,100);
00107   hspes = new TH1F("hspes","SPE if measured",200,0,40);
00108   hnpevar = new TH1F("hnpevar","~N PE input",500,0,500);
00109   hsignalmapP = new TH2F("hsignalmapP","Mean(Response) - Mean(Pedestal) HFP;#eta;#phi",26,28.5,41.5,36,0,72);
00110   hsignalRMSmapP = new TH2F("hsignalRMSmapP","RMS Response HFP;#eta;#phi",26,28.5,41.5,36,0,72);
00111   hnpemapP = new TH2F("hnpemapP","~N PE input HFP;#eta;#phi",26,28.5,41.5,36,0,72);
00112   hnpemapP->SetOption("COLZ");hsignalmapP->SetOption("COLZ");hsignalRMSmapP->SetOption("COLZ");
00113   hsignalmapM = new TH2F("hsignalmapM","Mean(Response) - Mean(Pedestal) HFM;#eta;#phi",26,-41.5,-28.5,36,0,72);
00114   hsignalRMSmapM = new TH2F("hsignalRMSmapM","RMS Response HFM;#eta;#phi",26,-41.5,-28.5,36,0,72);
00115   hnpemapM = new TH2F("hnpemapM","~N PE input HFM;#eta;#phi",26,-41.5,-28.5,36,0,72);
00116   hnpemapM->SetOption("COLZ");hsignalmapM->SetOption("COLZ");hsignalRMSmapM->SetOption("COLZ");
00117   // Channel-by-channel histos
00118   for (int i=0;i<13;i++) for (int j=0;j<36;j++) for (int k=0;k<2;k++) {
00119     if (i>10 && j%2==0) continue;
00120     sprintf(htit,"ts_+%d_%d_%d",i+29,j*2+1,k+1);
00121     hts[i][j][k] = new TH1F(htit,htit,10,-0.5,9.5); // TimeSlices (pulse shape)
00122     sprintf(htit,"tsmean_+%d_%d_%d",i+29,j*2+1,k+1);
00123     htsm[i][j][k] = new TH1F(htit,htit,100,0,10);   // Mean signal time estimated from TS 
00124     sprintf(htit,"sp_+%d_%d_%d",i+29,j*2+1,k+1);
00125     hsp[i][j][k] = new TH1F(htit,htit,1201,-25,30000); // Big-scale spectrum (linear ADC)
00126     sprintf(htit,"spe_+%d_%d_%d",i+29,j*2+1,k+1);
00127     hspe[i][j][k] = new TH1F(htit,htit,200,-9.5,190.5); // Small-scale spectrum (linear ADC)
00128     sprintf(htit,"ped_+%d_%d_%d",i+29,j*2+1,k+1);
00129     hped[i][j][k] = new TH1F(htit,htit,200,-9.5,190.5); // Pedestal spectrum
00130     sprintf(htit,"ts_-%d_%d_%d",i+29,j*2+1,k+1);
00131     hts[i+13][j][k] = new TH1F(htit,htit,10,-0.5,9.5);
00132     sprintf(htit,"tsmean_-%d_%d_%d",i+29,j*2+1,k+1);
00133     htsm[i+13][j][k] = new TH1F(htit,htit,100,0,10);  
00134     sprintf(htit,"sp_-%d_%d_%d",i+29,j*2+1,k+1);
00135     hsp[i+13][j][k] = new TH1F(htit,htit,1201,-25,30000);
00136     sprintf(htit,"spe_-%d_%d_%d",i+29,j*2+1,k+1);
00137     hspe[i+13][j][k] = new TH1F(htit,htit,200,-9.5,190.5); 
00138     sprintf(htit,"ped_-%d_%d_%d",i+29,j*2+1,k+1);
00139     hped[i+13][j][k] = new TH1F(htit,htit,200,-9.5,190.5); 
00140   } 
00141   // PIN-diodes histos
00142   for (int i=0;i<4;i++) for (int j=0;j<3;j++) {
00143     sprintf(htit,"ts_PIN%d_+Q%d",j+1,i+1);
00144     htspin[i][j] = new TH1F(htit,htit,10,-0.5,9.5);
00145     sprintf(htit,"sp_PIN%d_+Q%d",j+1,i+1);
00146     hsppin[i][j] = new TH1F(htit,htit,1601,-25,40000);
00147     sprintf(htit,"spe_PIN%d_+Q%d",j+1,i+1);
00148     hspepin[i][j] = new TH1F(htit,htit,200,-9.5,190.5);
00149     sprintf(htit,"ped_PIN%d_+Q%d",j+1,i+1);
00150     hpedpin[i][j] = new TH1F(htit,htit,200,-9.5,190.5);
00151     sprintf(htit,"tsmean_PIN%d_+Q%d",j+1,i+1);
00152     htsmpin[i][j] = new TH1F(htit,htit,100,0,10);  
00153     sprintf(htit,"ts_PIN%d_-Q%d",j+1,i+1);
00154     htspin[i+4][j] = new TH1F(htit,htit,10,-0.5,9.5);
00155     sprintf(htit,"sp_PIN%d_-Q%d",j+1,i+1);
00156     hsppin[i+4][j] = new TH1F(htit,htit,1601,-25,40000);
00157     sprintf(htit,"spe_PIN%d_-Q%d",j+1,i+1);
00158     hspepin[i+4][j] = new TH1F(htit,htit,200,-9.5,190.5);
00159     sprintf(htit,"ped_PIN%d_-Q%d",j+1,i+1);
00160     hpedpin[i+4][j] = new TH1F(htit,htit,200,-9.5,190.5);
00161     sprintf(htit,"tsmean_PIN%d_-Q%d",j+1,i+1);
00162     htsmpin[i+4][j] = new TH1F(htit,htit,100,0,10);  
00163   }
00164   std::cout<<std::endl<<"histfile="<<histfile.c_str()<<"  textfile="<<textfile.c_str()<<std::endl;
00165   return;
00166 }
00167 
00168 void HistSpecs(TH1F* hist, Double_t &mean, Double_t &rms, Double_t range=4) {
00169   Double_t xmin,xmax;
00170   mean=hist->GetMean();
00171   rms=hist->GetRMS();
00172   xmin=hist->GetXaxis()->GetXmin();
00173   xmax=hist->GetXaxis()->GetXmax();
00174   hist->SetAxisRange(mean-range*rms-100,mean+range*rms+100);
00175   mean=hist->GetMean();
00176   rms=hist->GetRMS();
00177   hist->SetAxisRange(mean-range*rms-100,mean+range*rms+100);
00178   mean=hist->GetMean();
00179   rms=hist->GetRMS();
00180   hist->SetAxisRange(xmin,xmax);
00181   return;
00182 }
00183 
00184 Double_t FitFun(Double_t *x, Double_t *par) { 
00185 // Spectra fit function: Pedestal Gaussian + asymmetric 1PE and 2PE peaks
00186 
00187   Double_t sum,xx,A0,C0,r0,sigma0,mean1,sigma1,A1,C1,r1,mean2,sigma2,A2,C2,r2;
00188   const Double_t k0=2.0,k1=1.0, k2=1.2;
00189 
00190   xx=x[0];
00191   sigma0 = par[2];
00192   A0 = 2*Nev/(2+2*par[0]+par[0]*par[0]);
00193   r0 = ((xx-par[1])/sigma0);
00194   C0 = 1/(sigma0* TMath::Exp(-k0*k0/2)/k0 +
00195           sigma0*sqrt(2*3.14159)*0.5*(1+TMath::Erf(k0/1.41421)));
00196   //sum = 1/(sqrt(2*3.14159)*par[2])*A0*TMath::Exp(-0.5*r0*r0);
00197   if(r0 < k0) sum = C0*A0*TMath::Exp(-0.5*r0*r0);
00198   else sum = C0*A0*TMath::Exp(0.5*k0*k0-k0*r0);
00199 
00200   mean1 = par[1]+par[3];
00201   //sigma1 = par[4];
00202   sigma1 = 1.547+0.125*par[3]+0.004042*par[3]*par[3];
00203   sigma1 = (sigma1+(9.1347e-3+3.845e-2*par[3])*par[4]*1.0)*par[2];
00204   A1 = A0*par[0];
00205   C1 = 1/(sigma1* TMath::Exp(-k1*k1/2)/k1 +
00206           sigma1*sqrt(2*3.14159)*0.5*(1+TMath::Erf(k1/1.41421)));
00207   r1 = ((xx-mean1)/sigma1);
00208   if(r1 < k1) sum += C1*A1*TMath::Exp(-0.5*r1*r1);
00209   else sum += C1*A1*TMath::Exp(0.5*k1*k1-k1*r1);
00210 
00211   mean2 = 2*par[3]+par[1];
00212   sigma2 = sqrt(2*sigma1*sigma1 - par[2]*par[2]);
00213   A2 = A0*par[0]*par[0]/2;
00214   C2 = 1/(sigma2* TMath::Exp(-k2*k2/2)/k2 +
00215           sigma2*sqrt(2*3.14159)*0.5*(1+TMath::Erf(k2/1.41421)));
00216   r2 = ((xx-mean2)/sigma2);
00217   if(r2 < k2) sum += C2*A2*TMath::Exp(-0.5*r2*r2);
00218   else sum += C2*A2*TMath::Exp(0.5*k2*k2-k2*r2);
00219 
00220   return sum;
00221 }
00222 
00223 void HFLightCal::endJob(void)
00224 {
00225   Double_t mean,rms,meanped,rmsped,maxc,npevar,npevarm;
00226   Double_t par[5],dspe=0,dnpe;
00227   Int_t tsmax,intspe;
00228   std::cout<<std::endl<<"HFLightCal endJob --> ";
00229   fprintf(tFile,"#RunN %d   Events processed %d",runN,eventN);
00230 
00231   for (int i=0;i<26;i++) for (int j=0;j<36;j++) for (int k=0;k<2;k++) {
00232     if (i>10 && i<13 && j%2==0) continue;
00233     if (i>23 && j%2==0) continue;
00234     meanped=rmsped=mean=rms=0;
00235     if (hsp[i][j][k]->Integral()>0) {
00236       HistSpecs(hped[i][j][k],meanped,rmsped);
00237       HistSpecs(hsp[i][j][k],mean,rms);
00238       if (hspe[i][j][k]->Integral()>hsp[i][j][k]->Integral()*0.9 || mean<100) {
00239         HistSpecs(hspe[i][j][k],mean,rms);
00240       }
00241       hsignalmean->Fill(mean);
00242       hsignalrms->Fill(rms);
00243       hpedmean->Fill(meanped);
00244       hpedrms->Fill(rmsped);
00245     }
00246   }
00247 
00248   meanped=hpedmean->GetMean();
00249   rmsped=hpedrms->GetMean();
00250   mean=hsignalmean->GetMean();
00251   rms=hsignalrms->GetMean();
00252   fprintf(tFile,"   MeanInput=<%.2f> [linADCcount]   RMS=<%.2f>\n",mean,rms);
00253   fprintf(tFile,"#eta/phi/depth  sum4maxTS     RMS      ~N_PE  sum4lowTS     RMS  maxTS  SPE +/- Err   Comment\n");
00254   TF1* fPed = new TF1("fPed","gaus",0,120);
00255   fPed->SetNpx(200);
00256   TF1 *fTot = new TF1("fTot",FitFun ,0,200,5);
00257   fTot->SetNpx(800);
00258   for (int i=0;i<26;i++) for (int j=0;j<36;j++) for (int k=0;k<2;k++) {
00259     if (i>10 && i<13 && j%2==0) continue;
00260     if (i>23 && j%2==0) continue;
00261     HistSpecs(hped[i][j][k],meanped,rmsped);
00262     HistSpecs(hsp[i][j][k],mean,rms);
00263     par[3]=0;
00264     if (hspe[i][j][k]->Integral()>hsp[i][j][k]->Integral()*0.9 || mean<100) {
00265       HistSpecs(hspe[i][j][k],mean,rms);
00266       if (hspe[i][j][k]->Integral(1,(int) (meanped+3*rmsped+12))/Nev>0.1) {
00267         //if (hspe[i][j][k]->Integral()>100 && mean-meanped<100) {
00268         if (mean+rms*3-meanped-rmsped*3>2 && rmsped>0) { // SPE fit if low intensity>0
00269           par[1] = meanped;
00270           par[2] = rmsped;
00271           par[0] = hped[i][j][k]->GetMaximum();
00272           fPed->SetParameters(par);
00273           hped[i][j][k]->Fit(fPed,"BQ0");
00274           fPed->GetParameters(&par[0]);
00275           hped[i][j][k]->Fit(fPed,"B0Q","",par[1]-par[2]*3,par[1]+par[2]*3);
00276           fPed->GetParameters(par);
00277           hped[i][j][k]->Fit(fPed,"BLIQ","",par[1]-par[2]*3,par[1]+par[2]*3);
00278           fPed->GetParameters(&par[0]);
00279           Nev = (int) hspe[i][j][k]->Integral();
00280           par[0]=0.1;
00281           par[3]=10;
00282           par[4]=6;
00283           fTot->SetParameters(par);
00284           fTot->SetParLimits(0,0,2);
00285           //fTot->FixParameter(1,par[1]);
00286           fTot->SetParLimits(1,par[1]-1,par[1]+1);
00287           fTot->FixParameter(2,par[2]);
00288           fTot->SetParLimits(3,1.2,100);
00289           fTot->SetParLimits(4,-1.64,1.64);
00290           hspe[i][j][k]->Fit(fTot,"BLEQ","");
00291           fTot->GetParameters(par);
00292           dspe=fTot->GetParError(3);
00293           dnpe=fTot->GetParError(0);
00294           if (par[3]<1.21 || dnpe>par[0]) par[3]=-1;
00295           else if (par[0]>1.96 || par[3]>49) par[3]=0;
00296           else {
00297             hspes->Fill(par[3]);
00298           }
00299         } 
00300       }
00301     }
00302 
00303     // NPE
00304     npevar=0;
00305     if (par[3]>0) npevar=par[0];                          // NPE from SPE fit
00306     else {                                                // NPE from high intensity signal
00307       if (hspe[i][j][k]->Integral()>hsp[i][j][k]->Integral()*0.98) {
00308         HistSpecs(hspe[i][j][k],mean,rms,3);
00309       }
00310       else {
00311         HistSpecs(hsp[i][j][k],mean,rms,3);
00312       }
00313       if (rmsped>0) {
00314         if (rms*rms-rmsped*rmsped>1 && mean>meanped) {
00315           npevar=(mean-meanped)*(mean-meanped)/(rms*rms-rmsped*rmsped);
00316         }
00317         else if (mean<100) {
00318           intspe=hspe[i][j][k]->Integral();
00319           hspe[i][j][k]->SetAxisRange(meanped+rmsped*4,300);
00320           npevar=hspe[i][j][k]->Integral()/intspe;
00321           if (npevar>0.01) npevar=-1;
00322           else npevar=0;
00323           hspe[i][j][k]->SetAxisRange(-20,300);
00324         }
00325       }
00326     }
00327     if (npevar>5.0e-5) hnpevar->Fill(npevar);
00328 
00329     if (i<13) {
00330       hsignalmapP->Fill(i+28.6+k/2.0,j*2+1,mean-meanped); 
00331       hsignalRMSmapP->Fill(i+28.6+k/2.0,j*2+1,rms);
00332       if (npevar>0) hnpemapP->Fill(i+28.6+k/2.0,j*2+1,npevar);
00333       fprintf(tFile,"%3d%4d%5d  %11.2f%8.2f",i+29,j*2+1,k+1,mean,rms);
00334     }
00335     else {
00336       fprintf(tFile,"%3d%4d%5d  %11.2f%8.2f",13-i-29,j*2+1,k+1,mean,rms);
00337       hsignalmapM->Fill(13-i-28.6-k/2.0,j*2+1,mean-meanped);
00338       hsignalRMSmapM->Fill(13-i-28.6-k/2.0,j*2+1,rms);
00339       if (npevar>0) hnpemapM->Fill(13-i-28.6-k/2.0,j*2+1,npevar);
00340     }
00341     if (npevar>0) fprintf(tFile,"  %9.4f",npevar);
00342     else  fprintf(tFile,"      0    ");
00343     fprintf(tFile,"   %8.2f%8.2f",meanped,rmsped);
00344     tsmax=hts[i][j][k]->GetMaximumBin()-1;
00345     fprintf(tFile," %4d",tsmax);
00346     if (par[3]>0 && par[3]<99)  fprintf(tFile,"%8.2f%7.2f",par[3],dspe);
00347     else if (npevar>0) fprintf(tFile,"%8.2f    0  ",(mean-meanped)/npevar);
00348     else fprintf(tFile,"     0      0  ");
00349 
00350     // Diagnostics 
00351     fprintf(tFile,"    ");
00352     if (hsp[i][j][k]->GetEntries()<=0) fprintf(tFile,"NoSignal\n");
00353     else if (hsp[i][j][k]->GetEntries()<=10) fprintf(tFile,"Nev<10\n");
00354     else {
00355       if (hsp[i][j][k]->Integral()<=10 || mean>12000)  fprintf(tFile,"SignalOffRange\n");
00356       else {
00357         if (hsp[i][j][k]->Integral()<100)  fprintf(tFile,"Nev<100/");
00358         if (npevar>0 && par[3]>0 && (npevar*Nev<10 || npevar<0.001)) 
00359           fprintf(tFile,"LowSignal/");
00360         else if (npevar==0 && fabs(mean-meanped)<3) fprintf(tFile,"LowSignal/");
00361         if (par[3]<0)  fprintf(tFile,"BadFit/");
00362         else if (par[3]==0)  fprintf(tFile,"NoSPEFit/");
00363         else if (par[3]>0 && npevar>1)   fprintf(tFile,"NPE>1/");
00364         if (npevar<0)   fprintf(tFile,"Problem/");
00365         if (mean<2) fprintf(tFile,"LowMean/");
00366         if (rms<0.5) fprintf(tFile,"LowRMS/"); 
00367         if (meanped<-1) fprintf(tFile,"LowPed/");
00368         else if (meanped>25) fprintf(tFile,"HighPed/"); 
00369         if (rmsped<0.5 && rmsped>0) fprintf(tFile,"NarrowPed/"); 
00370         else if (rmsped>10) fprintf(tFile,"WidePed/");
00371         if (hped[i][j][k]->GetBinContent(201)>10) fprintf(tFile,"PedOffRange"); 
00372         fprintf(tFile,"-\n");
00373       }
00374     }
00375   }
00376 
00377   for (int i=0;i<8;i++) for (int j=0;j<3;j++) {
00378     HistSpecs(hpedpin[i][j],meanped,rmsped);
00379     HistSpecs(hsppin[i][j],mean,rms);
00380     if (hspepin[i][j]->Integral()>hsppin[i][j]->Integral()*0.9 || mean<100) {
00381       HistSpecs(hspepin[i][j],mean,rms);
00382     }
00383     if (i<4) fprintf(tFile," PIN%d  +Q%d  %12.2f  %6.2f",j+1,i+1,mean,rms);
00384     else fprintf(tFile," PIN%d  -Q%d  %12.2f  %6.2f",j+1,i-3,mean,rms);
00385     fprintf(tFile,"  %15.2f  %6.2f",meanped,rmsped);
00386     tsmax=htspin[i][j]->GetMaximumBin()-1;
00387     fprintf(tFile,"  %4d\n",tsmax);
00388   } 
00389 
00390   mFile->Write();
00391   mFile->Close();
00392   fclose(tFile);
00393   std::cout<<std::endl<<" --endJob-- done"<<std::endl;
00394   return;
00395 }
00396 
00397 void HFLightCal::analyze(const edm::Event& fEvent, const edm::EventSetup& fSetup) {
00398 
00399   // event ID
00400   edm::EventID eventId = fEvent.id();
00401   int runNer = eventId.run ();
00402   int eventNumber = eventId.event ();
00403   if (runN==0) runN=runNer;
00404   eventN++;
00405   if (verbose) std::cout << "========================================="<<std::endl
00406                          << "run/event: "<<runNer<<'/'<<eventNumber<<std::endl;
00407 
00408   Double_t buf[20];
00409   Double_t maxADC,signal,ped,meant;
00410   Int_t ii,maxisample=0,i1=3,i2=6;
00411 
00412   // HF PIN-diodes
00413   edm::Handle<HcalCalibDigiCollection> calib;  
00414   fEvent.getByType(calib);
00415   if (verbose) std::cout<<"Analysis-> total CAL digis= "<<calib->size()<<std::endl;
00416   /* COMMENTED OUT by J. Mans (7-28-2008) as major changes needed with new Calib DetId 
00417   for (unsigned j = 0; j < calib->size (); ++j) {
00418     const HcalCalibDataFrame digi = (*calib)[j];
00419     HcalElectronicsId elecId = digi.elecId();
00420     HcalCalibDetId calibId = digi.id();
00421     if (verbose) std::cout<<calibId.sectorString().c_str()<<" "<<calibId.rbx()<<" "<<elecId.fiberChanId()<<std::endl;
00422     int isector = calibId.rbx()-1;
00423     int ipin = elecId.fiberChanId();
00424     int iside = -1;
00425     if (calibId.sectorString() == "HFP") iside = 0; 
00426     else if (calibId.sectorString() == "HFM") iside = 4;
00427     maxisample = itspinmax[isector+iside][ipin]-1;
00428 
00429     if (iside != -1) {
00430       for (int isample = 0; isample < digi.size(); ++isample) {
00431         int adc = digi[isample].adc();
00432         int capid = digi[isample].capid ();
00433         double linear_ADC = digi[isample].nominal_fC();
00434         if (verbose) std::cout<<"PIN linear_ADC = "<<linear_ADC<<"  atMAXTS="<<maxisample<<std::endl;
00435         htspin[isector+iside][ipin]->Fill(isample,linear_ADC);
00436         buf[isample]=linear_ADC;
00437       }
00438       i1=maxisample-1;
00439       i2=maxisample+2;
00440       if (i1<0) {i1=0;i2=3;}
00441       else if (i2>9) {i1=6;i2=9;} 
00442       if      (i1==0) ped=buf[8]+buf[9]+buf[6]+buf[7];
00443       else if (i1==1) ped=buf[8]+buf[9]+buf[6]+buf[7];
00444       else if (i1==2) ped=buf[0]+buf[1]+buf[6]+buf[7];
00445       else if (i1==3) ped=buf[0]+buf[1]+buf[2]+buf[7];
00446       else if (i1>=4) ped=buf[0]+buf[1]+buf[2]+buf[3];
00447       signal=0;
00448       for (ii=0;ii<4;ii++) signal+=TMath::Max(buf[ii+i1],ped/4.0);
00449       hsppin[isector+iside][ipin]->Fill(signal);
00450       hspepin[isector+iside][ipin]->Fill(signal);
00451       hpedpin[isector+iside][ipin]->Fill(ped);
00452 
00453       // Mean signal time estimation
00454       ped=ped/4;
00455       meant=0;
00456       for (ii=0;ii<4;ii++) meant+=(TMath::Max(buf[ii+i1],ped)-ped)*(ii+i1);
00457       if (signal-ped*4>0) meant/=(signal-ped*4); 
00458       else meant=i1+1;
00459       htsmpin[isector+iside][ipin]->Fill(meant);
00460     }
00461   }
00462   */  
00463   // HF
00464   edm::Handle<HFDigiCollection> hf_digi;
00465   fEvent.getByType(hf_digi);
00466   if (verbose) std::cout<<"Analysis-> total HF digis= "<<hf_digi->size()<<std::endl;
00467 
00468   for (unsigned ihit = 0; ihit < hf_digi->size (); ++ihit) {
00469     const HFDataFrame& frame = (*hf_digi)[ihit];
00470     HcalDetId detId = frame.id();
00471     int ieta = detId.ieta();
00472     int iphi = detId.iphi();
00473     int depth = detId.depth();
00474     if (verbose) std::cout <<"HF digi # " <<ihit<<": eta/phi/depth: "
00475                            <<ieta<<'/'<<iphi<<'/'<< depth << std::endl;
00476 
00477     if (ieta>0) ieta = ieta-29;
00478     else ieta = 13-ieta-29;
00479 
00480     maxADC=-99;
00481     for (int isample = 0; isample < frame.size(); ++isample) {
00482       int adc = frame[isample].adc();
00483       int capid = frame[isample].capid ();
00484       double linear_ADC = frame[isample].nominal_fC();
00485       double nominal_fC = detId.subdet () == HcalForward ? 2.6 *  linear_ADC : linear_ADC;
00486 
00487       if (verbose) std::cout << "Analysis->     HF sample # " << isample 
00488                              << ", capid=" << capid 
00489                              << ": ADC=" << adc 
00490                              << ", linearized ADC=" << linear_ADC
00491                              << ", nominal fC=" << nominal_fC <<std::endl;
00492 
00493       hts[ieta][(iphi-1)/2][depth-1]->Fill(isample,linear_ADC);
00494       buf[isample]=linear_ADC;
00495       /*
00496       if (maxADC<linear_ADC) {
00497         maxADC=linear_ADC;
00498         maxisample=isample;
00499       }
00500       */
00501     }
00502 
00503     maxADC=-99;
00504     for (int ii=0; ii<10; ii++) {
00505       signal=buf[ii];
00506       if      (ii<2) signal -= (buf[ii+4]+buf[ii+8])/2.0;
00507       else if (ii<4) signal -= buf[ii+4];
00508       else if (ii<6) signal -= (buf[ii+4]+buf[ii-4])/2.0;
00509       else if (ii<8) signal -= buf[ii-8];
00510       else           signal -= (buf[ii-4]+buf[ii-8])/2.0;
00511       if (signal>maxADC) {
00512         maxADC=signal;
00513         maxisample=ii;
00514       }
00515     }
00516 
00517     if (abs(maxisample-itsmax[ieta][(iphi-1)/2][depth-1]-1)>1)  maxisample=itsmax[ieta][(iphi-1)/2][depth-1]-1;
00518     if (verbose) std::cout<<eventNumber<<"/"<<ihit<<": maxTS="<<maxisample<<endl;
00519 
00520     // Signal = four capIDs found by PreAnal, Pedestal = four capIDs off the signal
00521     htmax->Fill(maxisample);
00522     i1=maxisample-1;
00523     i2=maxisample+2;
00524     if (i1<0) {i1=0;i2=3;}
00525     else if (i2>9) {i1=6;i2=9;} 
00526     signal=buf[i1]+buf[i1+1]+buf[i1+2]+buf[i1+3];
00527     hsp[ieta][(iphi-1)/2][depth-1]->Fill(signal);
00528     hspe[ieta][(iphi-1)/2][depth-1]->Fill(signal);
00529     /*
00530     if      (i1==0) ped=(buf[4]+buf[8])/2.0+(buf[5]+buf[9])/2.0+buf[6]+buf[7];
00531     else if (i1==1) ped=(buf[0]+buf[8])/2.0+(buf[5]+buf[9])/2.0+buf[6]+buf[7];
00532     else if (i1==2) ped=(buf[0]+buf[8])/2.0+(buf[1]+buf[9])/2.0+buf[6]+buf[7];
00533     else if (i1==3) ped=(buf[0]+buf[8])/2.0+(buf[1]+buf[9])/2.0+buf[2]+buf[7];
00534     else if (i1==4) ped=(buf[0]+buf[8])/2.0+(buf[1]+buf[9])/2.0+buf[2]+buf[3];
00535     else if (i1==5) ped=(buf[0]+buf[4])/2.0+(buf[1]+buf[9])/2.0+buf[2]+buf[3];
00536     else if (i1==6) ped=(buf[0]+buf[4])/2.0+(buf[1]+buf[5])/2.0+buf[2]+buf[3];
00537     */
00538     
00539     if      (i1<2) ped=buf[8]+buf[9]+buf[6]+buf[7];
00540     else if (i1==2) ped=buf[6]+buf[9]+buf[7]+buf[0];
00541     else if (i1==3) ped=buf[0]+buf[1]+buf[7]+buf[2];
00542     else if (i1>=4) ped=buf[0]+buf[1]+buf[2]+buf[3];
00543     
00544     hped[ieta][(iphi-1)/2][depth-1]->Fill(ped);
00545 
00546     // Mean signal time estimation
00547     ped=ped/4;
00548     meant=(buf[i1]-ped)*i1+(buf[i1+1]-ped)*(i1+1)+(buf[i1+2]-ped)*(i1+2)+(buf[i1+3]-ped)*(i1+3);
00549     meant /= (buf[i1]-ped)+(buf[i1+1]-ped)+(buf[i1+2]-ped)+(buf[i1+3]-ped);
00550     htmean->Fill(meant);
00551     htsm[ieta][(iphi-1)/2][depth-1]->Fill(meant);
00552   }
00553 }
00554
HFLightCal.cc
