Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:14:40 

0001 
0002 
0003 // **************************************************************************
0004 
0005 tileHelper::tileHelper(){
0006   if(ACTIVECHANNELS) {
0007     cout<<"<tileHelper>:  Wrong instantiation  sequence  - hLabHelper already exists.  EXITING"<<endl;
0008     //    return;
0009   }
0010   //  ONLY ACTIVE channels will be copyed from .prdf to adc
0011   CHTOTAL        = TILECHANNELS+TILETRIGGERCH;
0012   ACTIVECHANNELS = CHTOTAL;
0013   channels       = ACTIVECHANNELS;
0014   samples        = NSAMPLES;
0015   fibers         = TILEFIBERS;
0016   detchannels    = TILECHANNELS;
0017     //  assumption is that out of all stored ADC channels the first detchannels 
0018     //  correspon to detector we study, the rest is everything else
0019 
0020   tileCalib    = new Double_t [detchannels];
0021   tileCalPeak  = new Double_t [detchannels];
0022   fiberLY      = new Double_t [fibers];
0023   fiberY       = new Double_t [fibers];
0024   yFit         = NULL;
0025   fLY          = NULL;
0026   for(int ifb = 0;   ifb<fibers;   ifb++) {fiberY[ifb] = 15.-2.5*(ifb+1); if(ifb) fiberY[ifb] -= 2.5;}
0027 
0028 }
0029 
0030 // **************************************************************************
0031 
0032 int tileHelper::evLoop(int run, int evToProcess, int fToProcess){
0033   hLabHelper * hlHelper = hLabHelper::getInstance();
0034   hlHelper             -> setPRDFRun(run);
0035   tileTree::getInstance();
0036 
0037   int OK;
0038 
0039 
0040   Eventiterator *it =  new fileEventiterator(hlHelper->prdfName, OK);
0041   if (OK)
0042     {
0043       cout <<"<evLoop> Couldn't open input file " << hlHelper->prdfName << endl;
0044       delete(it);
0045       return 0;
0046     }
0047 
0048   Event *evt;
0049   int eventseq(0);
0050   // Loop over events
0051   while ( (evt = it->getNextEvent()) != 0 )
0052     {
0053       hlHelper->eventseq = evt->getEvtSequence();
0054 
0055       if ( evt->getEvtType() != 1 ) {
0056     cout << "eventseq " << hlHelper->eventseq << " event type = " << evt->getEvtType() << endl;
0057     continue;
0058       }
0059       // ev_display->Divide(nx_c,ny_c);
0060       // ev_display->SetBit(TPad::kClearAfterCR);
0061 
0062       // Get sPHENIX Packet
0063       Packet_hbd_fpgashort *spacket = dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101) );
0064       if ( spacket )
0065     {
0066       spacket->setNumSamples( NSAMPLES );
0067       //      int nmod_per_fem = spacket->iValue(0,"NRMODULES");
0068       //cout << "nmod_per_fem " << nmod_per_fem << endl;
0069       if(hlHelper->eventsread%1000==0)  cout <<"RUN  "<<hlHelper->runnumber<<"  EVENT "<<hlHelper->eventseq <<" - "<<hlHelper->eventsread<<endl;
0070       hlHelper->eventsread++;
0071       for (int ich=0; ich<ACTIVECHANNELS; ich++)    {
0072         for (int is=0; is<NSAMPLES; is++)       {
0073           hlHelper->adc[ich][is] =  spacket->iValue(hlHelper->active[ich],is);
0074           //          cout<<"<evLoop  Event "<< hclt->eventseq<<"  ich/is/adc  "<<ich<<" - "<<is<<" - "<< hlHelper->adc[ich][is]<<" active "<<hlHelper->active[ich]<<endl;
0075         }
0076       }
0077       hlHelper->collect();
0078       collectTileSummary();
0079       tileTiming();
0080       delete spacket;
0081     }
0082 
0083       delete evt;
0084       hlHelper-> thcl -> Fill();
0085       if ((evToProcess>0&&hlHelper->eventsread>=evToProcess)||hlHelper->eventsread>270000) break;
0086     }
0087   // Normalize trigger histogramm to get per event hit count (number of channels above threshold) 
0088   int nbs = runsum.trhits->GetNbinsX()-1;
0089   for(int bx = 1; bx<=nbs; bx++){
0090     float normx = runsum.trhits->GetBinContent(bx);
0091     if(normx<=0.) continue;
0092     for (int threshold = 1; threshold<=CHANNELTHRESHOLDS; threshold++) {
0093       float cont = runsum.hits_tpc->GetBinContent(bx,threshold);
0094       runsum.hits_tpc->SetBinContent(bx, threshold, cont/normx);    
0095     }
0096   }
0097   
0098   if (hlHelper->eventsread>1) runsum.hits_tpc->Scale(1./(Double_t)hlHelper->eventsread);
0099   // status();
0100   cout<<"ALLDONE for RUN  "<<hlHelper->runnumber<<"  Events "<<hlHelper->eventsread<<endl;
0101   delete it;
0102   hlHelper->thcl -> Write();
0103   hlHelper->fhcl -> Write();
0104   hlHelper->fhcl -> Close();
0105 
0106   //  closeLoop();
0107   return eventseq;
0108 
0109 }  
0110 
0111 // **************************************************************************
0112 void tileHelper::status(){
0113   hLabHelper * hlHelper = hLabHelper::getInstance();
0114   int nx_c = ACTIVECHANNELS<=12? 2  :  4;             
0115   int ny_c = ACTIVECHANNELS/nx_c+ACTIVECHANNELS%nx_c;
0116   if(!(fiberDisplay=(TCanvas*)(gROOT->FindObject("fiberDisplay")))) {
0117 
0118     TString fD = "Fiber Amplitudes (black - raw, red - fit) Run "; fD += runId.Data();
0119     //  TString fD = "Fiber Amplitudes (black - raw, red - fit, blue - integral) Run "; fD += hlHelper->prdfName;
0120     fiberDisplay = new TCanvas("fiberDisplay",fD,400*nx_c,200*ny_c);
0121     fiberDisplay->Divide(nx_c, ny_c);
0122   }
0123   Double_t ymx[ACTIVECHANNELS], xmx[ACTIVECHANNELS];
0124   for (int ich=0; ich<ACTIVECHANNELS; ich++){
0125     Double_t rvmax, fvmax, ivmax(0), rvrms, fvrms, ivrms(0), rvmean, fvmean, ivmean(0);
0126     rvmax = hlHelper->rpeak[ich]->GetMaximum(); 
0127     fvmax = hlHelper->fm[ich]->GetMaximum();
0128     //    ivmax = hlHelper->fint[ich]->GetMaximum();  
0129     rvrms = hlHelper->rpeak[ich]->GetRMS(); 
0130     fvrms = hlHelper->fm[ich]->GetRMS();
0131     //    ivrms = hlHelper->fint[ich]->GetRMS();
0132     rvmean = hlHelper->rpeak[ich]->GetMean(); 
0133     fvmean = hlHelper->fm[ich]->GetMean();
0134     //    ivmean = hlHelper->fint[ich]->GetMean();
0135 
0136     Double_t ym = std::max(max(rvmax, fvmax), ivmax);                            
0137     ymx[ich] = (int)(log10(ym)); ymx[ich] = pow(10., ymx[ich]);
0138     while(ymx[ich]<ym) ymx[ich] *=2;
0139     Double_t xm = std::max(max(rvmean+4*rvrms, fvmean+4*fvrms), ivmean+4*ivrms); xmx[ich] = (int)(log10(xm)); xmx[ich] = pow(10., xmx[ich]);
0140     while(xmx[ich]<xm) xmx[ich] *=2;
0141     cout<<"chan  "<<ich<<"  rvmax  "<<rvmax<<" fvmax "<<fvmax<<" ivmax "<<ivmax<<" xm "<<xm<<" xmx "<<xmx[ich]<<endl;  
0142   }
0143   Double_t xm(0.);
0144   for (int ich=0; ich<ACTIVECHANNELS; ich++){
0145     //  just to make x-scales on left and right for the same fiber identical
0146     hlHelper->rpeak[ich]->SetLineColor(1);    hlHelper->rpeak[ich]->SetMaximum(ymx[ich]);
0147     hlHelper->fm[ich]   ->SetLineColor(2);    hlHelper->fm[ich]    ->SetMaximum(ymx[ich]);
0148     hlHelper->fint[ich] ->SetLineColor(4);    hlHelper->fint[ich]  ->SetMaximum(ymx[ich]);
0149     if(!(ich%2)) xm = max(xmx[ich], xmx[ich+1]);
0150     cout<<"chan  "<<ich<<" xrange "<<xm<<endl;
0151     fiberDisplay->cd(ich+1); hlHelper->rpeak[ich]->SetAxisRange(0., xm); hlHelper->rpeak[ich]->Draw();
0152     hlHelper->fm[ich]->Draw("same");
0153     //    hlHelper->fint[ich]->Draw("same");
0154   }
0155   fiberDisplay->Update();
0156 }
0157 
0158 // **************************************************************************
0159 void tileHelper::updateMap(){
0160   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0161   int runnumber = hlHelper->runnumber;
0162   ACTIVECHANNELS  = 0;
0163   for (int ich = 0; ich<TILECHANNELS;    ich++) hlHelper->active[ich] = runnumber<803?   feech1[ich]   :   feech2[ich];
0164   ACTIVECHANNELS  += TILECHANNELS;
0165     for (int ich = 0; ich<TILETRIGGERCH; ich++) {
0166       if(runnumber<1125)                      {
0167     hlHelper->active[TILECHANNELS+ich] = -1;
0168       }      else if(runnumber>=1125&&runnumber<1152) {
0169     hlHelper->active[TILECHANNELS+ich]    =    trch1125[ich];
0170     if(hlHelper->active[TILECHANNELS+ich]>=0)  ACTIVECHANNELS ++;
0171       }      else                                     {
0172     hlHelper->active[TILECHANNELS+ich] = trch1152[ich];
0173     if(hlHelper->active[TILECHANNELS+ich]>=0)  ACTIVECHANNELS ++;
0174       }
0175     }
0176 }
0177 // **************************************************************************
0178 void tileHelper::updateCalibration(){
0179   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0180   int runnumber = hlHelper->runnumber;
0181   if(runnumber<864) {
0182     for (int ich = 0; ich<TILECHANNELS; ich++) tileCalib[ich] = sc_779[ich];
0183   } else if (runnumber<1061)  {
0184     for (int ich = 0; ich<TILECHANNELS; ich++) tileCalib[ich] = sc_900[ich];
0185   } else if (runnumber<1123)  {
0186     for (int ich = 0; ich<TILECHANNELS; ich++) tileCalib[ich] = sc_1061[ich];
0187   } else {
0188     for (int ich = 0; ich<TILECHANNELS; ich++) tileCalib[ich] = sc_1123[ich];
0189   }
0190 }
0191 
0192 // **************************************************************************
0193 
0194 void tileHelper::evreset(){
0195   evtsum.reset();
0196   uSum   = 0.;     cSum   = 0.;
0197   luSum  = 0.;     lcSum  = 0.;   
0198   ruSum  = 0.;     rcSum  = 0.;
0199   YF     = 0.;
0200   YU     = 0.;     YC     = 0.;
0201   XU     = 0.;     XC     = 0.;
0202   sumFU  = 0.;     sumFC = 0.;
0203   muxU   = -100.;  muxC  = -100.;  muxFU  = -100.;  muxFC  = -100.;
0204   muxUCh = 0;      muxCCh = 0;     muxUFiber = 0;   muxCFiber = 0;
0205   for(int ith = 0;   ith<CHANNELTHRESHOLDS; ith++)     hitcnt[ith]  = 0;
0206   for(int ifb = 0;   ifb<fibers;            ifb++)     fiberLY[ifb] = 0.;
0207   if(yFit) {delete yFit; yFit = NULL;}
0208   if(fLY)  {delete fLY;  fLY  = NULL;}
0209   
0210 }
0211 
0212 
0213 // **************************************************************************
0214 
0215 void tileHelper::tileTrigger(){
0216   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0217   for(int ich = 0; ich<detchannels; ich++) { 
0218     for(int ith = 0; ith<CHANNELTHRESHOLDS; ith++) {
0219       if(int(hlHelper->rawPeak[ich]/TRIGGERRES) >= ith) hitcnt[ith]++; else break;
0220     }
0221   }
0222   //  trigger hit multiplicity plot
0223   runsum.trhits -> Fill(uSum/TRIGGERRES);
0224   for(int ith=0; ith<CHANNELTHRESHOLDS; ith++) {if(hitcnt[ith]) runsum.hits_tpc->Fill(uSum/TRIGGERRES, ith,  hitcnt[ith]) ; else break;}
0225 
0226 
0227 }
0228 
0229 // **************************************************************************
0230 
0231 void tileHelper::tileTiming(){
0232   //  build eventsum (adc sum without pedestal subtruction)
0233   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0234   for(int is = 0; is < NSAMPLES; is++){
0235     evtsum.evtadcsum[is] = 0.;
0236     for (int ich = 0; ich < detchannels; ich++) { evtsum.evtadcsum[is] += hlHelper->adc[ich][is];}
0237   }
0238   for (int is = 0; is < NSAMPLES; is++)  evtsum.evtGraph->SetPoint(is,(Double_t)is,evtsum.evtadcsum[is]);
0239   fitTileSignal();
0240   //  do the fit and compare parameter distribution in individual channels and total sum
0241   //  minimize the number of free parameters and compare timing RMS for "all free" anf "Some constrained"
0242   ;
0243 }
0244 // **************************************************************************
0245 
0246 void tileHelper::tileDisplay(){
0247   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0248   //  build eventsum (adc sum without pedestal subtruction)
0249   for(int is = 0; is < NSAMPLES; is++){ 
0250     evtsum.evtadcsum[is] = 0.;
0251     for (int ich = 0; ich < detchannels; ich++) {
0252       evtsum.evtadcsum[is] += hlHelper->adc[ich][is];
0253     }
0254   }
0255   for (int is = 0; is < NSAMPLES; is++)  evtsum.evtGraph->SetPoint(is,(Double_t)is,evtsum.evtadcsum[is]);
0256   //  cout<<"<tileDisplay>   call to fitTileSignal"<<endl;
0257   fitTileSignal();
0258   TCanvas * tiledisplay = (TCanvas *)(gROOT->FindObject("tiledisplay"));
0259   if(!tiledisplay) tiledisplay  = new TCanvas("tiledisplay","sPhenix tile display", 300, 300, 800, 500);
0260   else tiledisplay->Clear();
0261   //  tiledisplay->Divide(1,2);
0262   //  tiledisplay->cd(1);
0263     //  tiledisplay->SetEditable(kTRUE);
0264   evtsum.evtGraph->SetMarkerStyle(20);
0265   evtsum.evtGraph->SetMarkerSize(0.4);
0266 
0267   evtsum.evtGraph->Draw("acp");
0268   tiledisplay->Update();
0269    //  tiledisplay->cd(2);
0270   evtsum.evtShape->Draw("same");
0271   tiledisplay->Update();
0272   //  tiledisplay->Clear();
0273 }
0274 
0275 // **************************************************************************
0276 
0277 void tileHelper::tileTriggerDisplay(){
0278   //  hLabHelper  * hlHelper = hLabHelper ::getInstance();
0279   int nx_c = ACTIVECHANNELS<=12?   2   :  4;             
0280   int ny_c = ACTIVECHANNELS/nx_c+ACTIVECHANNELS%nx_c;
0281   if(!(triggerDisplay=(TCanvas*)(gROOT->FindObject("triggerDisplay")))) {
0282     TString fD = "Fiber Amplitudes (black - raw, red - fit) Run "; fD += runId.Data();
0283     //  TString fD = "Fiber Amplitudes (black - raw, red - fit, blue - integral) Run "; fD += hlHelper->prdfName;
0284     triggerDisplay = new TCanvas("triggerDisplay",fD,400*nx_c,200*ny_c);
0285     triggerDisplay->Divide(nx_c, ny_c);
0286   }
0287 }
0288 // **************************************************************************
0289 
0290 //  let's now see if all this game with patterning was worth trying
0291 void tileHelper::tilePattern(Int_t nx, Int_t ny, Int_t run, Int_t mod){
0292   TH2  * pcPat(NULL);           //  pattern of pixel conts measured in different tile areas
0293   TH1 ** pcPatH(NULL);          //  Collections of  hits from different tile areas 
0294 //TF1 ** pcPatF(NULL);          //  Fits to the histograms of hits in different tile areas
0295   xdivisions            = nx;
0296   ydivisions            = ny;
0297   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0298   TFile *f;
0299   if(!(f=hlHelper->attachrootrun(run))) {
0300     cout<<"<ERROR> pattern(): root file for runNumber "<<hlHelper->runnumber<<" not in the list"<<endl;
0301     return;
0302   }
0303   f->cd();
0304   gStyle         ->SetOptFit();
0305   // hn         = "Total Pixel Count vs impact position  RunId = "; hn += Id;      hn += ";Calibated X [cm];Calibrated Y [cm];Pixels";
0306   // pc_imp     = new TH3F("pc_imp",  hn,   25, 0., 25., 15, 0., 15., 1000, 0, 200.);  
0307   if((lyPattern  = (TCanvas*)(gROOT->FindObject("lyPattern")))) delete lyPattern;
0308   if((lyFits     = (TCanvas*)(gROOT->FindObject("lyFits"))))    delete lyFits;
0309 
0310   TString lyp = "Pattern of light yield "; lyp += runId.Data();
0311   TString lyd = "LYD "; lyd += runId.Data();
0312   lyPattern = new TCanvas(lyd,   lyp, 400, 400);
0313   TString lyv = "LIV "; lyv += runId.Data();
0314   TString lyf = "Fits to the measurenments of light yield "; lyf += (mod==0? "(pc_imp) " : "(pc_fimp)"); lyf += runId.Data();
0315   lyFits    = new TCanvas("lyv", lyf, 200*nx, 200*ny);
0316   lyFits    ->Divide(xdivisions, ydivisions);
0317 
0318   TH3 * pcimp = (TH3F *) (gROOT->FindObject(mod==0? "pc_imp" : "pc_fimp"));
0319     //  TH3 * pc_imp = (TH3F *) (gROOT->FindObject("pc_imp"));
0320   if(!pcimp) {
0321     cout<<"<ERROR> pattern(): "<<(mod==0? "(pc_imp) " : "(pc_fimp)")<<" Scatterplot is not found"<<endl;
0322     return;
0323   }
0324   Int_t nximp  = pcimp->GetNbinsX();                Int_t    nyimp = pcimp->GetNbinsY(); 
0325   Int_t nzimp  = pcimp->GetNbinsZ();
0326   //  Double_t dx  = pcimp->GetXaxis()->GetBinWidth(1); Double_t dy    = pcimp->GetYaxis()->GetBinWidth(1);
0327   Double_t dxr = tileSizeX / nx;                     Double_t dyr   = tileSizeY / ny;
0328   //  total number of cells to combine 
0329   Int_t dnx    = nximp/nx;                              Int_t    dny   = nyimp/ny;
0330   if(dnx*nx!=nximp || dny*ny!=nyimp) {
0331     cout<<"Patterning request does not match your binning in the pixel count distribution pcimp"<<endl;
0332     return;
0333   }
0334   pcPat = (TH2F *) (gROOT->FindObject("pcPat"));
0335   if(pcPat) delete pcPat;
0336   pcPat       = new TH2F("pcPat", " fitted positions of maxima in pcimp pixel distributions", nx, 0., 25., ny, 0., 15.);
0337   pcPatH      = new TH1 * [nx*ny];
0338   for (int iy = 0; iy<ny; iy++){
0339     for (int ix = 0; ix<nx; ix++){
0340       Int_t  cpad = nx*(ny-iy-1) + ix +1;
0341       TString hn = "h_"; hn += ix; hn += "_"; hn += iy;
0342       TString ht = "Pixel count. Range X = "; ht += (dxr*ix);  ht += " | ";  ht += (dxr*(ix+1)); 
0343       ht += " Y = ";  ht += (dyr*iy);  ht += " | "; ht += (dyr*(iy+1));
0344       if((pcPatH[iy*nx+ix]=(TH1 *) (gROOT->FindObject(hn)))) delete pcPatH[iy*nx+ix];
0345       pcPatH[iy*nx+ix] = new TH1D(hn, ht, nzimp, 0., pcimp->GetZaxis()->GetXmax()); 
0346       //  fill that histogramm
0347       pcimp->ProjectionZ(hn, 1+dnx*ix, dnx*(ix+1), 1+dny*iy, dny*(iy+1));
0348       // draw that histogramm
0349       lyFits          -> cd(cpad);
0350       pcPatH[iy*nx+ix]-> Draw();
0351       // check for total number of entries in the histogram, if less then minimum - go to next one
0352       Int_t   entries = pcPatH[iy*nx+ix]->GetEntries();
0353       if(entries<minProjEntries) continue;
0354       //  rebin the histogramm to get at least 50 hits at the peak
0355       Double_t norm;
0356       while (((norm=pcPatH[iy*nx+ix]->GetMaximum())<50.||
0357           pcPatH[iy*nx+ix]->GetBinCenter(pcPatH[iy*nx+ix]->GetMaximumBin())<=10)&&
0358          pcPatH[iy*nx+ix]->GetNbinsX()>50 ) {
0359     Int_t combine = 2; 
0360     while((pcPatH[iy*nx+ix]->GetNbinsX()/combine)*combine!=pcPatH[iy*nx+ix]->GetNbinsX()) combine++;
0361     pcPatH[iy*nx+ix] = pcPatH[iy*nx+ix]->Rebin(combine);
0362       }
0363       //  redrow the histogramm
0364       pcPatH[iy*nx+ix]-> Draw();
0365       //  fit that histogramm
0366       Double_t mean = pcPatH[iy*nx+ix]->GetMean();
0367       Double_t rms  = pcPatH[iy*nx+ix]->GetRMS();
0368       Double_t p[6];
0369       p[0]    = norm;  p[1]  = mean;  p[2]  = rms;  p[3] = 0.;  p[4] = 0.;  p[5] = 0.;
0370       TF1 * fg  = new TF1("fg","[0]*TMath::Gaus(x,[1],[2])", 15., pcimp->GetZaxis()->GetXmax());
0371       fg       -> SetParameters(p);
0372       fg       -> SetParLimits(0, 0.1*norm,10.*norm);
0373       fg       -> SetParLimits(1, mean-rms, mean+rms);
0374       fg       -> SetParLimits(2, rms/2., 2.*rms);
0375       pcPatH[iy*nx+ix]->Fit(fg, "r");
0376       fg       ->GetParameters(p);
0377       TF1 * fgb = new TF1("fgb","[0]*TMath::Gaus(x,[1],[2])+[3]+[4]*x+[5]*x*x", 15., pcimp->GetZaxis()->GetXmax());
0378       fgb      ->SetParameters(p);
0379       fgb      ->SetParLimits(1, p[1]-rms,p[1]+rms);
0380       fgb      ->SetParLimits(2, rms/2,2.*rms);
0381       fgb      ->SetParLimits(3, 0., p[0]);
0382       pcPatH[iy*nx+ix]->Fit(fgb, "r");
0383       lyFits          ->Update();
0384       fgb       ->GetParameters(p);
0385       //  fill the patten 2D plot
0386       pcPat->Fill(dxr*(ix+0.5), dyr*(iy+0.5), p[1]);      
0387     }
0388   }
0389   lyPattern->cd();
0390   pcPat->Draw("lego2");
0391 }
0392 
0393 // **************************************************************************
0394 
0395 void tileHelper::collectTileSummary(){
0396   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0397   //  tile specific data
0398   for (int ich=0;  ich<detchannels;  ich++) {
0399     tileCalPeak[ich]      = hlHelper->fitPeak[ich]/tileCalib[ich];
0400     evtsum.times[ich]     = hlHelper->fitTime[ich];
0401     evtsum.chi2[ich]      = hlHelper->fitChi2[ich];
0402     uSum  += hlHelper->rawPeak[ich];
0403     cSum  += tileCalPeak[ich];
0404     runsum.cdata ->Fill(ich, tileCalPeak[ich]);
0405     if(ich%2) {
0406 
0407       // right ends
0408       ruSum += hlHelper->rawPeak[ich];
0409       rcSum += tileCalPeak[ich];
0410     } else {
0411       // lefy ends
0412       luSum += hlHelper->rawPeak[ich];
0413       lcSum += tileCalPeak[ich];
0414     }
0415   }
0416   //  Maximum amplitude in detector
0417   for (int ich=0; ich<detchannels; ich++){
0418     if(hlHelper->rawPeak[ich] >=muxU)   {muxUCh = ich; muxU = hlHelper->rawPeak[ich];}
0419     if(tileCalPeak[ich]       >=muxC)   {muxCCh = ich; muxC = tileCalPeak[ich];}
0420     if(ich%2) { 
0421       // hard wired: fiber Y's in the tiles
0422       int fiber = ich/2;
0423       float y   = fiberY[fiber];
0424       sumFU    += hlHelper->rawPeak[ich]; sumFC += tileCalPeak[ich];
0425       YU       += y*sumFU;      YC    += y*sumFC;
0426       if(sumFU >= muxFU) {muxUFiber = ich/2; muxFU = sumFU;}
0427       if(sumFC >= muxFC) {muxCFiber = ich/2; muxFC = sumFC;}
0428       fiberLY[ich/2]  = sumFC;
0429     } else {
0430       sumFU =  hlHelper->rawPeak[ich]; sumFC  = tileCalPeak[ich];
0431     }
0432   }
0433   YU /= uSum;   YC /= cSum;  YF = getYFit();
0434   evtsum.hitfiber   = muxCFiber;
0435   evtsum.hitfiberpc = muxFC;
0436   evtsum.yChi2      = yFit->GetChisquare()/yFit->GetNDF();
0437   evtsum.ySigma     = yFit->GetParameter(2);
0438   runsum.rmax  -> Fill(muxUCh, muxU);  runsum.cmax  -> Fill(muxCCh, muxC);
0439   runsum.rfsum -> Fill(muxUFiber, muxFU);
0440   runsum.cfsum -> Fill(muxCFiber, muxFC);
0441   runsum.rtsum -> Fill(muxUFiber, uSum);
0442   runsum.ctsum -> Fill(muxCFiber, cSum);
0443   ras   = ((hlHelper->rawPeak[muxUFiber*2+1]-hlHelper->rawPeak[muxUFiber*2])  /(hlHelper->rawPeak[muxUFiber*2]+hlHelper->rawPeak[muxUFiber*2+1]));
0444   cas   = ((tileCalPeak[muxCFiber*2+1]-tileCalPeak[muxCFiber*2])/(tileCalPeak[muxCFiber*2]+tileCalPeak[muxCFiber*2+1]));
0445   runsum.rfasym-> Fill(muxUFiber, ras);
0446   runsum.cfasym-> Fill(muxCFiber, cas);
0447   rtas  = (ruSum - luSum)  /(ruSum + luSum);
0448   ctas  = (rcSum - lcSum)  /(rcSum + lcSum);
0449   runsum.rtasym-> Fill(rtas);
0450   runsum.ctasym-> Fill(ctas);
0451   XU = 12.5+12.5*rtas; /* XC = 12.5+12.5*ctas; */   XC = AsymToX(ctas);
0452   runsum.XY    -> Fill(XC, YF);
0453   runsum.ry    -> Fill(YU);
0454   runsum.cy    -> Fill(YC);
0455   runsum.fy    -> Fill(YF);
0456   runsum.rx    -> Fill(XU);
0457   runsum.cx    -> Fill(XC);
0458   runsum.rimp  -> Fill(XU, YU);
0459   runsum.cimp  -> Fill(XC, YC);
0460   runsum.cx_pc -> Fill(XC, cSum);
0461   runsum.cy_pc -> Fill(YC, cSum);
0462   runsum.fy_pc -> Fill(YF, cSum);
0463   runsum.cimpW -> Fill(XC, YC, cSum);
0464   for(int ith=0; ith<CHANNELTHRESHOLDS; ith++) {if(hitcnt[ith]) runsum.treff->Fill(hitcnt[ith], ith, cSum) ; else break;}
0465   runsum.pc_imp ->Fill(XC, YC, cSum);
0466   runsum.pc_fimp->Fill(XC, YF, cSum);
0467   evtsum.evtpc      = cSum;
0468   evtsum.hitfibertd = evtsum.times[muxCFiber*2+1]-evtsum.times[muxCFiber*2];
0469   //  cout<<"<collectTileSummary> YF/tChi2/yChi2  "<<YF<<" * "<<evtsum.tChi2<<" * "<< evtsum.yChi2<<endl;
0470   runsum.y_yChi2 -> Fill(YF,           evtsum.yChi2);
0471   runsum.s_y_Chi2-> Fill(evtsum.tChi2, evtsum.yChi2);
0472   reject();
0473   runsum.trcode  -> Fill(evtsum.rejectCode);
0474   runsum.y_ys    -> Fill(YF,           evtsum.ySigma);
0475   if(evtsum.ySigma<4.5) {
0476     runsum.yCleaned -> Fill(YF);
0477     runsum.yc_rcode -> Fill(YF,evtsum.rejectCode); 
0478     if(!evtsum.rejectCode<100) runsum.yKept -> Fill(YF);
0479   }
0480   if(ACTIVECHANNELS>=TILECHANNELS+2){
0481     Double_t tower_pc = hlHelper->fitPeak[9]/0.05;     //  pixels in the tower
0482     runsum.td_r_l   -> Fill(evtsum.hitfibertd, evtsum.hitfiberpc);
0483     runsum.td_t_tw  -> Fill(evtsum.evttime-hlHelper->fitTime[9], evtsum.evtpc, tower_pc);
0484   }
0485   tileTrigger();
0486 }
0487 
0488 // **************************************************************************
0489 
0490 Double_t tileHelper::getYFit(){
0491   //  Let's do Y first
0492   yFit = new TF1("yFit", "[0]*exp(-abs(x-[1])/[2])", 0., tileSizeY);
0493   Double_t    par[]  =  {20., tileSizeY/2., 1.};
0494   yFit -> SetParameters(par);
0495   yFit -> SetParLimits(1, 0., tileSizeY);
0496   yFit -> SetParLimits(2, 0., 5.);
0497   fLY  = new TGraph(fibers, fiberY, fiberLY); 
0498   fLY -> Fit(yFit, "rq", "nq", 0., tileSizeY);
0499   yFit->GetParameters(par);
0500   runsum.Y0 -> Fill(par[0]);
0501   runsum.Y1 -> Fill(par[1]);
0502   runsum.Y2 -> Fill(par[2]);
0503   
0504   return par[1];
0505 }
0506 
0507 // **************************************************************************
0508 
0509 void tileHelper::tileImpact(){
0510   //  Let's do Y first
0511   yFit = new TF1("yFit", "[0]*exp(-abs(x-[1])/[2])", 0., tileSizeY);
0512   Double_t    par[]  =  {20., tileSizeY/2., 1.};
0513   yFit -> SetParameters(par);
0514   yFit -> SetParLimits(1, 0., tileSizeY);
0515   yFit -> SetParLimits(2, 0., 4.);
0516   fLY  = new TGraph(fibers, fiberY, fiberLY); 
0517   fLY -> Fit(yFit, "rq", "nq", 0., tileSizeY);
0518   yFit->GetParameters(par);
0519   runsum.Y0 -> Fill(par[0]);
0520   runsum.Y1 -> Fill(par[1]);
0521   runsum.Y2 -> Fill(par[2]);
0522   //  Now X
0523   Float_t slope(0.);
0524   Float_t dsl  = 0.5/20;
0525   for (int isl = 0;  isl< 20; isl++){
0526     slope = dsl/2.+dsl*isl;
0527     Double_t X = 12.5*(1 + ctas/slope); 
0528     runsum.XvsSl->Fill(X, slope);
0529     
0530   }
0531 }
0532 // **************************************************************************
0533 
0534 void tileHelper::fitTileSignal(){
0535   //  zero approximation to fit parameters
0536   //  use pulse data to find location of maximum
0537   //  hLabHelper  * hlHelper = hLabHelper ::getInstance();
0538   int      rTime = 0;
0539   Double_t rVal  = 0.;
0540   for (int iSample = RISETIME-1; iSample <= NSAMPLES-FALLTIME; iSample++) {
0541     if(evtsum.evtadcsum[iSample]>rVal) {
0542       rVal  = evtsum.evtadcsum[iSample];
0543       rTime = iSample;    
0544     }
0545   }
0546   rVal         -= PEDESTAL*detchannels;
0547   if(rVal<0.) rVal = 0.;
0548   Double_t  par0[6];
0549   par0[0] = rVal/3.;
0550   par0[1] = max(0.,(Double_t)(rTime-RISETIME));
0551   par0[2] = 4.;
0552   par0[3] = 1.6;
0553   par0[4] = PEDESTAL*detchannels;  //   we do not do pedestal subtrastion at this time
0554   par0[5] = 0;      //   slope of the pedestals is initially set to "0"
0555   evtsum.evtShape->SetParameters(par0);
0556   for(int parn=0; parn<6; parn++) evtsum.evtShape->SetParLimits(parn, par0Min[parn], par0Max[parn]);
0557   evtsum.evtShape->SetParLimits(4, PEDESTAL*(detchannels-1), PEDESTAL*(detchannels+1));
0558   //  fitFunc->cd(chan+1);  //  shapes[chan]->Draw();
0559   //  gPad->Modified();
0560   evtsum.evtGraph->Fit(evtsum.evtShape, "RQWM", "NQ", 0., (Double_t)NSAMPLES);
0561   evtsum.evtShape->GetParameters(evtsum.evtfitpar);
0562   evtsum.evtpedestal  = (evtsum.evtfitpar[4]+evtsum.evtfitpar[5]*rTime);   //  
0563   evtsum.evtSignal->SetParameters(evtsum.evtfitpar);
0564   evtsum.evtSignal->SetParameter(4, 0.);
0565   evtsum.evtSignal->SetParameter(5, 0.);
0566   evtsum.evtpeak  = evtsum.evtSignal->GetMaximum(par0Min[1], par0Max[1]);
0567   evtsum.evttime  = evtsum.evtSignal->GetMaximumX(par0Min[1], par0Max[1]);
0568   evtsum.tChi2    = evtsum.evtShape->GetChisquare()/evtsum.evtShape->GetNDF();
0569   //  cout<<"<fitTileSignal>  Peak  "<< evtsum.evtpeak<<"  Time "<<evtsum.evttime<<"  Chi2 "<<evtsum.tChi2 <<endl;
0570   //  ev_display->cd(chan+1);
0571   //  shapes[hlHelper->chan]->Draw("same");
0572 
0573 }
0574 
0575 // **************************************************************************
0576 
0577 tileHelper::runtilesummary::runtilesummary(){
0578   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0579   hlHelper-> fhcl-> cd();
0580   TString hn         = "Raw data: Events peaking in channel  RunId = ";  hn += runId;      hn += ";Channel #;ADC Peak (events with max in this channel) [counts]";
0581   rmax       = new TH2F("rmax",    hn,   TILECHANNELS, 0, TILECHANNELS,   2050, 0., 2050.);
0582   hn         = "Total Pixel Count with tighter triggering  RunId = "; hn += runId; hn += ";HitMultThr;HitAmplThr [trigger counts];Pixels";
0583   treff      = new TH3F("treff",   hn, HITMULTTHRESHOLDS, 0, HITMULTTHRESHOLDS, CHANNELTHRESHOLDS, 0, CHANNELTHRESHOLDS, 1000, 0., 200.); 
0584   hn         = "Total Pixel Count vs impact position  RunId = "; hn += runId;      hn += ";Calibated X [cm];Calibrated Y [cm];Pixels";
0585   pc_imp     = new TH3F("pc_imp",   hn,   250, 0., 25., 150, 0., 15., 1000, 0, 200.);  
0586   hn        += " Y Fitted";
0587   pc_fimp    = new TH3F("pc_fimp",  hn,   250, 0., 25., 150, 0., 15., 1000, 0, 200.);  
0588   hn         = "Fit to the Pixel Count vs impact position  RunId = "; hn += runId; hn += ";Calibated X [cm];Calibrated Y [cm];Pixels";
0589   pc_pat     = new TH2F("pc_pat",  hn,   25, 0., 25., 15, 0., 15.);  
0590   hn        += " Y Fitted";  
0591   pc_fpat    = new TH2F("pc_fpat", hn,   25, 0., 25., 15, 0., 15.);  
0592   hn         = "Raw Fiber sum (R+L), Peak in this fiber  RunId = ";   hn += runId; hn += ";Fiber #;Fiber ADC Sum (events with max in this fiber) [counts]";
0593   rfsum      = new TH2F("rfsum",   hn,   TILEFIBERS, 0, TILEFIBERS, 2050, 0., 2050.);  
0594   hn         = "Raw Total sum, Peak in this fiber  RunId = ";  hn += runId;        hn += ";Fiber #;Event ADC Sum (events with max in this fiber) [counts]";
0595   rtsum      = new TH2F("rtsum",   hn,   TILEFIBERS, 0, TILEFIBERS, 2050, 0., 2050.); 
0596   hn         = "Raw Fiber Asymmetry  RunId = ";  hn += runId;                      hn += ";Fiber #;Fiber R/L Asymmetry";
0597   rfasym     = new TH2F("rfasym",  hn,   TILEFIBERS, 0, TILEFIBERS, 200, -1., 1.);  
0598   hn         = "Raw Event Asymmetry  RunId = ";  hn += runId;                      hn += ";Event Raw R/L Event Asymmetry";
0599   rtasym     = new TH1F("rtasym",  hn,   200, -1., 1.);
0600   hn         = "Raw X  RunId = ";   hn += runId;                                   hn += ";Event Raw X-coordinate [cm]";
0601   rx         = new TH1F("rx",      hn,   250, 0., 25.);
0602   hn         = "Raw Y  RunId = ";   hn += runId;                                   hn += ";Event Raw Y-coordinate [cm]";
0603   ry         = new TH1F("ry",      hn,   150, 0., 15.);
0604   hn         = "Raw impact position  RunId = ";  hn += runId;                      hn += ";Raw X [cm];Raw Y [cm]";
0605   rimp       = new TH2F("rimp",    hn,   250,  0., 25., 150,  0., 15.); 
0606 
0607   hn         = "Calibrated data: Fitted Maxima  RunId = ";  hn += runId;           hn += ";Channel #; Calibrated Fit Peak [pixels]";
0608   cdata      = new TH2F("cdata",   hn,   TILECHANNELS, 0, TILECHANNELS,   1000, 0., 200.);
0609   hn         = "Calibrated data: Events peaking in channel  RunId = ";  hn += runId;  hn += ";Channel #;Calibrated Fit Peak (evens with max in this channel) [pixels]";
0610   cmax       = new TH2F("cmax",    hn,   TILECHANNELS, 0, TILECHANNELS,   1000, 0., 200.);
0611   hn         = "Calibrated Fiber sum (R+L), Peak in the fiber  RunId = ";  hn += runId;  hn += ";Fiber #;Calibrated Fiber Sum (events with max in this fiber) [pixels]";
0612   cfsum      = new TH2F("cfsum",   hn,   TILEFIBERS, 0, TILEFIBERS, 1000, 0., 200.);    
0613   hn         = "Calibrated Total sum, Peak in this fiber  RunId = ";  hn += runId;    hn += ";Fiber #;Calibrated Total Sum (events with max in this fiber) [pixels]";
0614   ctsum      = new TH2F("ctsum",   hn,   TILEFIBERS, 0, TILEFIBERS, 1000, 0., 200.);  
0615   hn         = "Calibrated Fiber Asymmetry  RunId = ";  hn += runId;               hn += ";Fiber #;Fiber R/L Asymmetry (calibrated)";
0616   cfasym     = new TH2F("cfasym",  hn,   TILEFIBERS, 0, TILEFIBERS, 200, -1., 1.);  
0617   hn         = "Calibrated Event Asymmetry  RunId = ";  hn += runId;               hn += ";Event R/L Event Asymmetry (calibrated)";
0618   ctasym     = new TH1F("ctasym",  hn,   200, -1., 1.);
0619   hn         = "Calibrated X  RunId = ";   hn += runId;                            hn += ";Calibrated X [cm]";
0620   cx         = new TH1F("cx",      hn,   250, 0., 25.);
0621   hn         = "Calibrated Y  RunId = ";   hn += runId;                            hn += ";Calibrated Y [cm]";
0622   cy         = new TH1F("cy",      hn,   150, 0., 15.);
0623   hn         = "Fitted Y  RunId = ";   hn += runId;                                hn += ";Fitted Y [cm]";
0624   fy         = new TH1F("fy",      hn,   150, 0., 15.);
0625   hn         = "Calibrated Impact Position  RunId = ";  hn += runId;               hn += ";Calibrated X [cm];Calibrated Y [cm]";
0626   cimp       = new TH2F("cimp",    hn,   250,  0., 25., 150,  0., 15.); 
0627   hn         = "Calibrated Pixel Count (Y) vs X-impact  RunId = ";  hn += runId;   hn +=";Calibrated X [cm];Calibrated Total Sum [pixels]";
0628   cx_pc      = new TH2F("cx_pc",   hn,   250,  0., 25., 1000, 0., 200.);
0629   hn         = "Calibrated Pixel Count (Y) vs Y-impact  RunId = ";  hn += runId;   hn +=";Calibrated Y [cm];Calibrated Total Sum [pixels]";
0630   cy_pc      = new TH2F("cy_pc",   hn,   150,  0., 15., 1000, 0., 200.);
0631   hn         = "Calibrated Pixel Count (Y) vs Fitted Y-impact  RunId = ";  hn += runId;   hn +=";Fitted Y [cm];Calibrated Total Sum [pixels]";
0632   fy_pc      = new TH2F("fy_pc",   hn,   150,  0., 15., 1000, 0., 200.);
0633   hn         = "Calibrated Weighted Impact Position  RunId = ";  hn += runId;      hn +=";Calibrated X [cm];Calibrated Y [cm]";
0634   cimpW      = new TH2F("cimpW",   hn,   250,  0., 25., 150,  0., 15.); 
0635   
0636   //  Trigger related data
0637   //  Pixel count above threshold vs total pixel count (based on calibrated data)
0638   hn         = "Uncalibrated trigger sum [trigger counts]  RunId = ";  hn += runId; hn += runId; hn += ";Trigger hits;";
0639   trhits     = new TH1F("trhits",  hn, 500, 0., 500.);
0640   hn         = "Number of hits in triggered events vs channel threshold(Y) and total trigger sum (X)  RunId = ";  hn += runId; hn += ";TotalSum [trigger counts];Trigger THr [trigger counts]";
0641   hits_tpc   = new TH2F("hits_tpc",hn, 500, 0., 100., CHANNELTHRESHOLDS, 0, CHANNELTHRESHOLDS); 
0642 
0643   //  
0644   //   debugging muon impact computations
0645   //  TH2        * XvsSl, * XY;     //  X value vs assumption about slope in the X vs Asymmetry 
0646   //  TH1        * y0, *y1, *y2;    //  Parameters from the tile-Y fit
0647   Y0       = new TH1F("Y0", "Y-coordinate fit P(0)", 200, 0., 200.);
0648   Y1       = new TH1F("Y1", "Y-coordinate fit P(1)", 100, 0., 15.);
0649   Y2       = new TH1F("Y2", "Y-coordinate fit P(2)", 100, 0., 10.);
0650   y_yChi2  = new TH2F("y_yChi2",  "yFit: Chi2 at different fitted Y values", 150,0.,15., 100,0.,50.);
0651   s_y_Chi2 = new TH2F("s_y_Chi2", "Chi2 of tile signal vs Chi2 of Y coordinate", 100, 0., 500., 100, 0., 50.);
0652   y_ys     = new TH2F("y_ys",     "yFit sigma vs fitted Y", 150,0.,15. ,100, 0., 5.);
0653   yCleaned = new TH1F("yCleaned", "Fitted Y after event rejection (fit sigma <4.5)", 150,0.,15.);
0654   yKept    = new TH1F("yKept",    "Fitted Y after event rejection (s<4.5, rejectcode=0)", 150,0.,15.);
0655   XvsSl    = new TH2F("XvsSl",   "X-Coordinate computed using different asymmentry slopes", 275, -15., 40., 20, 0., 0.5);
0656   XY       = new TH2F("XY",      "Muon Impact:  x  vs  Y. Fit Y vs Best X", 225, -15., 40., 75, 0., 15.);
0657   td_r_l   = new TH2F("td_r_l",  "R/L Timing difference in hit fiber vs PixelCount in Hit Fiber", 200, -1, 1., 200, 0., 100.);
0658   td_t_tw  = new TH3F("td_t_tw", "Tile-to-tower time difference(X) vs Pixel Counts in tile(Y) and tower(Z)", 200, -1. ,1., 100, 0., 100., 100, 0., 100.);
0659   trcode   = new TH1F("trcode",  "Tile reject codes ", 1000, 0., 1000.);
0660   yc_rcode = new TH2F("yc_rcode","Tile reject code vs Cleaned Y ", 150, 0., 15., 1000, 0., 1000.);
0661 
0662 }
0663 // **************************************************************************
0664 
0665 Int_t tileHelper::reject(){
0666   //  test fit results in individual channels, compare to fit of the tile sum
0667   //  hLabHelper  * hlHelper = hLabHelper ::getInstance();
0668   hLabHelper * hlHelper = hLabHelper::getInstance();
0669   evtsum.rejectCode = 0;
0670   for(int ich = 0; ich< TILECHANNELS; ich++) {if(evtsum.chi2[ich]>10.) evtsum.rejectCode ++;}
0671   for(int ich = 0; ich< TILECHANNELS; ich++) {if(abs(evtsum.times[ich]-evtsum.evttime)>2.) evtsum.rejectCode += 10;}
0672   for(int ich = 0; ich< TILECHANNELS; ich++) {if(hlHelper->rawPeak[ich]<2.) evtsum.rejectCode += 100;}
0673   return  evtsum.rejectCode;
0674 }
0675 
0676 
0677 // **************************************************************************
0678 
0679 tileLightYield::tileLightYield(){
0680   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0681   hlHelper->getFileLists();
0682   //  rootfile * rfile;
0683   list<rootfile>::iterator    next;
0684   for (next = hlHelper->roots.begin(); next != hlHelper->roots.end(); ++next){
0685     ;
0686   }
0687 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!