CMSSW Reference Manual

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