Prototype2/hcalLab/hcalHelper.C

0001
0002 // *************************************************************************
0003 hcalHelper::hcalHelper(){
0004   if(ACTIVECHANNELS) {
0005     cout<<"<hcalHelper::hcalHelper>  Wrong instantiation  sequence  - hLabHelper already exists.  EXITING"<<endl;
0006     delete this;
0007   }
0008   ACTIVECHANNELS = 192;   // just to insure everything works as it must
0009   hLabHelper::getInstance();
0010   eventseq    = 0;
0011   eventsread  = 0;
0012   displayMode = 3;   //   default - only final summary is displayed
0013   //  TODO:  needs better understanding - how to handle HG Gains
0014   //  create calorimeter
0015
0016   //  we create all graphic objects in the directory of temporary file (rootTempFileName)
0017
0018   t1044 = new hcal();    //  it does stacks configuring at this time
0019   //  ONLY ACTIVE channels will be copyed from .prdf to adc
0020   CHTOTAL        = detchannels;
0021   ACTIVECHANNELS = CHTOTAL;
0022   channels       = ACTIVECHANNELS;
0023   samples        = NSAMPLES;
0024   activeCh    = new Int_t   [ACTIVECHANNELS];
0025   adc         = new Int_t * [ACTIVECHANNELS];
0026   for (int ich=0; ich<ACTIVECHANNELS; ich++) {
0027     adc[ich]        = new Int_t  [NSAMPLES];
0028   }
0029   //  parse calibration information
0030   t1044->setCalibration();
0031   //  make a look-up table for easy access to adc channels
0032   Int_t  ilive = 0;
0033   for(int istk = 0; istk<CALSTACKS; istk++) {
0034     if(t1044->alive[istk]) {
0035       t1044->stacks[istk]->updateMap(ilive);
0036       Int_t stackid   = t1044->stacks[istk]->stackId;
0037       Int_t chtocp    = t1044->stacks[istk]->twrsInStack*t1044->stacks[istk]->gains;
0038       for(int ich = 0; ich < chtocp; ich++)  {
0039     activeCh[ilive] = (stackid==EMC?     emcCh[ich] :
0040                (stackid==HINNER?  hcalInnerCh[ich]  :
0041                 (stackid==HOUTER?  hcalOuterCh[ich]  :
0042                  (stackid==HODO?    hodoCh[ich]      :
0043                   counters[ich]))));
0044     ilive ++;
0045       }
0046       //      t1044->stacks[istk]->print();
0047     }
0048   }
0049   cout<<"<hcalHelper::hcalHelper>:  "<<ACTIVECHANNELS<<" Channels Initialized"<<endl;
0050   rdata        = new TH2F("rdata", "Raw Peak Values for all active channels",ACTIVECHANNELS, 0, ACTIVECHANNELS, 1024, -2048., 2048.);
0051
0052 }
0053
0054
0055 // **************************************************************************
0056 void hcalHelper::setRunKind(TString kind) {
0057   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0058   hcalHelper  * hHelper  = hcalHelper ::getInstance();
0059   hlHelper -> setRunKind(kind);
0060   if(kind=="cosmics") {
0061     hHelper->t1044->alive[EMC]   = kFALSE;
0062     hHelper->t1044->alive[HODO]  = kFALSE;
0063     hHelper->t1044->alive[SCINT] = kFALSE;
0064     delete hHelper->t1044->stacks[EMC];
0065     delete hHelper->t1044->stacks[HODO];
0066     delete hHelper->t1044->stacks[SCINT];
0067     hHelper->t1044->stacks[EMC]   = NULL;
0068     hHelper->t1044->stacks[HODO]  = NULL;
0069     hHelper->t1044->stacks[SCINT] = NULL;
0070   }
0071 }
0072
0073 // **************************************************************************
0074
0075 int hcalHelper::evLoop(int run, int evToProcess, int fToProcess){
0076   runnumber      = run;
0077   if(runnumber<=0) return 0;
0078   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0079   hlHelper->runnumber = runnumber;
0080   hlHelper -> setPRDFRun(runnumber, kFALSE);
0081   hlHelper -> makeCanvasDirectory();
0082   hcalTree::getInstance();
0083   //  hlHelper-> fhcl->cd();
0084
0085   int OK;
0086
0087   Eventiterator *it =  new fileEventiterator(hlHelper->prdfName, OK);
0088   if (OK)
0089     {
0090       cout <<"<evLoop> Couldn't open input file " << hlHelper->prdfName << endl;
0091       delete(it);
0092       return 0;
0093     }
0094   Event *evt;
0095   // Loop over events
0096   while ( (evt = it->getNextEvent()) != 0 )
0097     {
0098       t1044 -> readyForEvent = kFALSE;
0099       eventseq = evt->getEvtSequence();
0100       eventReject  = 0;
0101       eventTrigger = 0;
0102       if ( evt->getEvtType() != 1 ) {
0103     cout << "eventseq " << eventseq << " event type = " << evt->getEvtType() << endl;
0104     continue;
0105       }
0106
0107       // Get sPHENIX Packet
0108       Packet_hbd_fpgashort *sp21101 =
0109     dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101));
0110       Packet_hbd_fpgashort *sp21102 =
0111     dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21102));
0112       if(sp21101&&sp21102){
0113     sp21101->setNumSamples( NSAMPLES );  sp21102->setNumSamples( NSAMPLES );
0114     for (int ich=0; ich<ACTIVECHANNELS; ich++)    {
0115       Int_t spCh  = activeCh[ich]<144?  activeCh[ich] :activeCh[ich]-144;
0116       Packet_hbd_fpgashort* spacket = activeCh[ich]<144? sp21101  :  sp21102;
0117       for (int is=0; is<NSAMPLES; is++)     {
0118         adc[ich][is] =  spacket->iValue(spCh,is);
0119         //        cout<<"<evLoop  Event "<< hlHelper->eventseq<<"  ich/is/adc  "<<ich<<" - "<<is<<" - "<< hlHelper->adc[ich][is]<<" active "<<hlHelper->active[ich]<<endl;
0120       }
0121     }
0122     delete sp21101; delete sp21102;
0123       } else {
0124     cout << "EventSeq " << eventseq << " type = " << evt->getEvtType() <<" Missing packets " <<(int)sp21101<<"/"<<(int)sp21102<<endl;
0125     continue;
0126       }
0127       eventsread++;
0128       if(eventsread%1000==0)  cout <<"RUN  "<<runnumber<<"  EVENT "<<eventseq <<" - "<<eventsread<<endl;
0129       //  invert raw data in non-calorimeter stacks to make all signals negative
0130       //  subtract 2048 everywhere
0131       for (int istk = 0; istk<CALSTACKS; istk++){
0132     if(!t1044->alive[istk]) continue;
0133     for (int itw=0; itw<t1044->stacks[istk]->twrsInStack; itw++){
0134       for(int ig=0; ig<t1044->stacks[istk]->gains; ig++){
0135         int iadc = t1044->stacks[istk]->towers[itw]->adcChannel[ig];
0136         for(int is=0; is<NSAMPLES; is++){
0137           if(t1044->stacks[istk]->stackKind==CALOR) {
0138         adc[iadc][is] =  adc[iadc][is]-2048.;
0139           } else {
0140         adc[iadc][is] = (activeCh[iadc]%2? -(adc[iadc][is]-2048.)   :  adc[iadc][is]-2048.);
0141           }
0142           //  knowing where this value goes we may immediately set satFlags
0143           t1044->stacks[istk]->towers[itw]->graph[ig]->SetPoint(is,(Double_t)is,adc[t1044->stacks[istk]->towers[itw]->adcChannel[ig]][is]);
0144           if(adc[iadc][is]<undflow||adc[iadc][is]>ovrflow)
0145         t1044->stacks[istk]->towers[itw]->satFlag[ig] += 1;
0146         } //  is
0147       } //  ig
0148     } //  itw
0149       } //  istk
0150       //  amount of data displayed in per event and per run display
0151       //  0  - nothing displayed except Summary at the end of processing
0152       //  1  - digitizations for everything
0153       //  2  - lego plots of individual amplitudes in calorimeter and digitization for counter
0154       Int_t rejectST = t1044->getStackTiming();
0155       if(rejectST<3 || !triggerOn)  {  //  HCal is not entirely empty or I want to see all events
0156     Int_t rejectCR = collectRaw();
0157     Int_t rejectCL = t1044->update(displayMode, kFALSE);
0158     eventReject    = rejectST*1000 + rejectCR*100 + rejectCL;
0159     eventTrigger   = t1044->collectTrPrimitives();
0160     if(displayMode<3) cout<<"Event "<<eventseq<<"  Trigger "<<eventTrigger<<"  eventReject "<<eventReject<<"  "<<rejectST<<"/"<<rejectCR<<"/"<<rejectCL<<endl;;
0161     t1044->displayRaw(displayMode);
0162       }
0163
0164       delete evt;
0165       if(displayMode<3){
0166     //  if(std::cin.get()==(Int_t)"q") break;
0167     if(std::cin.get()==113) break;
0168     //      htree -> thcl -> Fill();
0169       }
0170       t1044->clean();
0171       if ((evToProcess>0&&eventsread>=evToProcess)||eventsread>270000) break;
0172     }
0173   //  display (print) un summary. Verbosity mode varies
0174   t1044->displaySummary(displayMode);
0175
0176
0177   cout<<"ALLDONE for RUN  "<<runnumber<<"  Events "<<eventsread<<endl;
0178   delete it;
0179   hlHelper -> thcl -> Write();
0180   hlHelper -> fhcl -> Write();
0181   //  TODO - if multiple files are processed closing the root file should be done somewhere else
0182   //  TODO - at the moment
0183   //  hlHelper -> renameAndCloseRF();
0184   //  closeLoop();
0185   return eventseq;
0186 }
0187
0188
0189 // **************************************************************************
0190 //  amount of per event and per run display
0191 //  0  - nothing displayed except Summary at the end of processing
0192 //  1  - digitizations for everything
0193 //  2  - lego plots of individual amplitudes in calorimeter and digitization for counters
0194
0195 void hcalHelper::setDisplayMode(Int_t mode){
0196   displayMode = mode;
0197   cout<<"<hcalHelper::setDisplayMode>  Mode is set to "<<displayMode<<endl;
0198 }
0199
0200 // **************************************************************************
0201
0202 void hcalHelper::setCTriggerOn(Bool_t  ON){
0203   triggerOn = ON;
0204   cout<<"<hcalHelper::setCTriggerOn>  triggerOn is set to "<<triggerOn<<endl;
0205 }
0206
0207 // **************************************************************************
0208
0209 void hcalHelper::setCTriggerGRange(Int_t gain){
0210   triggerGain = gain;
0211   TRGAINRANGE = gain;
0212   cout<<"<hcalHelper::setDCTriggerGRange>  Gaine is set to "<<triggerGain<<endl;
0213 }
0214
0215 // **************************************************************************
0216 void hcalHelper::updateCalibration(){
0217   //  hLabHelper  * hlHelper = hLabHelper ::getInstance();
0218   ;
0219 }
0220
0221 // **************************************************************************
0222
0223 void  hcalHelper::hcalPattern(Int_t nx, Int_t ny, Int_t run, Int_t mod){;}
0224
0225 // **************************************************************************
0226
0227 //  To call after global timing is extracted for individual stacks (note that scintillating counters are considered "stacks" - they are probably individually timed. Scintillating hodoscope is stack in it own rights
0228 Int_t  hcalHelper::collectRaw(){
0229   // currently no rejection at this stage
0230   rejectRaw = 0;
0231   //  stacks first
0232   for(int istk = 0; istk<CALSTACKS; istk++){
0233     if(!(t1044->stacks[istk])) continue;
0234     for(int ig = 0; ig<t1044->stacks[istk]->gains; ig++){
0235       //     Double_t fPeak = t1044->stacks[istk]->fitPeak[ig];
0236       Int_t  fTime = (int)(t1044->stacks[istk]->fitTime[ig]+0.5);
0237       for (int itw = 0; itw<t1044->stacks[istk]->twrsInStack; itw++){
0238     Int_t ich = t1044->stacks[istk]->towers[itw]->adcChannel[ig];
0239     //  extract raw peak position and peak value for this channel
0240     Int_t iss = max((fTime-2), 1);// Int_t ise = min((fTime+3), NSAMPLES);
0241     Double_t psum(0.); float sUsed(0.);
0242     for (int is = 0;  is<iss; is++) { if(adc[ich][is]>=undflow&&adc[ich][is]<=ovrflow) {psum += adc[ich][is]; sUsed += 1.;}}
0243
0244     if(sUsed) psum /= sUsed; else psum=0.;
0245     Double_t rVal = adc[ich][fTime] - psum;
0246     t1044->stacks[istk]->towers[itw]->rawPed[ig]  = psum;
0247     t1044->stacks[istk]->towers[itw]->rawAmpl[ig] = rVal;
0248     t1044->stacks[istk]->towers[itw]->rawTime[ig] = fTime;    // rTime;
0249     rdata -> Fill(t1044->stacks[istk]->towers[itw]->adcChannel[ig], rVal);
0250       }
0251
0252     }
0253   //  all kinds of scintillators etc
0254   }
0255   return rejectRaw;
0256 }
0257 // **************************************************************************
0258 void  hcalHelper::hcalImpact(){;}
0259 // **************************************************************************
0260 void  hcalHelper::fitHCalSignal(){;}
0261 // **************************************************************************
0262 void  hcalHelper::hcalDisplay(){;}
0263 // **************************************************************************
0264 //  FIXES to data
0265 void hcalHelper::dofixes(){
0266   // if(runnumber>=1125 && runnumber<1152){
0267   //   //  invert pulse in channel 134 (number 9 in todays parlance
0268   //   for (int is = 0; is < NSAMPLES; is++) {adc[8][is] -= 2*PEDESTAL; adc[8][is] = abs(adc[8][is]);}
0269   // }
0270   ;
0271 }
0272
0273 // **************************************************************************
0274
0275 Int_t hcalHelper::reject(){
0276   return false;
0277 }
0278
0279 // **************************************************************************
0280 hcal::hcal(){
0281
0282   maxStacks    = CALSTACKS;                                //  CALSTACKS
0283   stacks       = new stack *  [maxStacks];
0284   alive        = new Bool_t   [maxStacks];
0285   kinds        = new Int_t    [maxStacks];
0286   kinds[EMC]    = CALOR;      alive[EMC]     = kTRUE;    //  EMC
0287   kinds[HINNER] = CALOR;      alive[HINNER]  = kTRUE;    //  HInner
0288   kinds[HOUTER] = CALOR;      alive[HOUTER]  = kTRUE;    //  Outer HCal
0289   kinds[HODO]   = COUNTER;    alive[HODO]    = kTRUE;    //  HODOSCOPE
0290   kinds[SCINT]  = COUNTER;    alive[SCINT]   = kTRUE;    //  Beam counters
0291   kinds[CHER]   = COUNTER;    alive[CHER]    = kFALSE;   //  Cherenkov's
0292   activeStacks    = 0;
0293   activeCalStacks = 0;
0294
0295   for (int istk = 0; istk < maxStacks; istk++) stacks[istk] = NULL;
0296
0297   //  this is just a proyotype, it has well defined configuration
0298   detchannels      = 0;
0299   hgDetChannels    = 0;
0300   lgDetChannels    = 0;
0301   if(alive[EMC])      {
0302     stacks[EMC]     = new stack(EMC,    EMCTOWERS,   EMCCOLUMNS, EMCROWS);
0303     //    cout<<"MMM "<<stacks[EMC]->mapUpdated<<endl;
0304     stacks[EMC]->stackKind = kinds[EMC];
0305     activeStacks ++;
0306     activeCalStacks ++;
0307     detchannels    += EMCTOWERS*EMCGAINS;    //   only one gain (treat as LG)
0308     lgDetChannels  += EMCTOWERS*EMCGAINS;
0309     //    cout<<"MMM "<<stacks[EMC]->mapUpdated<<endl;
0310   }
0311   if(alive[HINNER])   {
0312     stacks[HINNER]  = new stack(HINNER, HCALTOWERS, HCALCOLUMNS, HCALROWS);
0313     stacks[HINNER]->stackKind = kinds[HINNER];
0314     activeStacks ++;
0315     activeCalStacks ++;
0316     detchannels    += HCALTOWERS*HCALGAINS;  //   dual gain readout in HCal
0317     hgDetChannels  += HCALTOWERS;
0318     lgDetChannels  += HCALTOWERS;
0319   }
0320   if(alive[HOUTER])   {
0321     stacks[HOUTER]  = new stack(HOUTER, HCALTOWERS, HCALCOLUMNS, HCALROWS);
0322     stacks[HOUTER]->stackKind = kinds[HOUTER];
0323     activeStacks ++;
0324     activeCalStacks ++;
0325     detchannels    += HCALTOWERS*HCALGAINS;
0326     hgDetChannels  += HCALTOWERS;
0327     lgDetChannels  += HCALTOWERS;
0328   }
0329   if(alive[HODO])   {
0330     stacks[HODO]  = new stack(HODO, T1044HODO, T1044HODO/2, 2);
0331     stacks[HODO]->stackKind = kinds[HODO];
0332     activeStacks ++;
0333     detchannels    += T1044HODO;
0334     lgDetChannels  += T1044HODO;
0335   }
0336   if(alive[SCINT])   {
0337     stacks[SCINT]  = new stack(SCINT, T1044TRIGGERCH, T1044TRIGGERCH/2, 2);
0338     stacks[SCINT]->stackKind = kinds[SCINT];
0339     activeStacks ++;
0340     detchannels    += T1044TRIGGERCH;
0341     lgDetChannels  += T1044TRIGGERCH;
0342   }
0343   for(int istk = 0; istk<CALSTACKS; istk++) {
0344     if(alive[istk]&&stacks[istk]->mapUpdated) cout<<"Stack "<<stacks[istk]->stackId<<" mapUpdated already set "<<stacks[istk]->mapUpdated<<endl;
0345   }
0346   readyForEvent = kTRUE;
0347   clean();
0348
0349   //  Summary histos for the whole calorimeter
0350   TString  totaN  = "tota";
0351   TString  totaT  = "Total Amplitude in Calorimeter System";
0352   totalamp = new TH1F(totaN,  totaT, 10000, 0., 10000.);
0353   TString  toteN  = "tote";
0354   TString  toteT  = "Total Energy in Calorimeter System";
0355   totale   = new TH1F(toteN,  toteT, 1000, 0., 60.);
0356   TString  tlcgN  = "tlcg";
0357   TString  tlcgT  = "Shower Center of gravity (units of Stacks)";
0358   totallcg = new TH1F(tlcgN,  tlcgT,  200, 0., 2.);
0359   TString  sclcgN  = "sclcg";
0360   TString  sclcgT  = "Shower Center in scintillators (units of Stacks)";
0361   scintlcg = new TH1F(sclcgN, sclcgT, 200, 0., 2.);
0362 }
0363
0364 // **************************************************************************
0365
0366 void hcal::setCalibration(){
0367   //  set energy scales
0368   for(int istk = 0; istk<CALSTACKS; istk++){
0369     if(!alive[istk]) continue;
0370     stacks[istk]->stackECalib = stECalib[stacks[istk]->stackId];
0371   //  TODO:   set correction coefficients in towers (currently unity)
0372     if(stacks[istk]->stackKind == CALOR) {
0373       for (int itw=0; itw< stacks[istk]->twrsInStack; itw++) {
0374     stacks[istk]->towers[itw]->twrECalib = stECalib[stacks[istk]->stackId];
0375     //  must be scaled by the Muon peak positions
0376     stacks[istk]->towers[itw]->twrEScale  = 1.;
0377       }
0378     }
0379   }
0380   //
0381 }
0382
0383 // **************************************************************************
0384 hcal::~hcal(){
0385   if(totalamp)  delete totalamp;
0386   if(totale)    delete totale;
0387   if(totallcg)  delete totallcg;
0388   if(scintlcg)  delete scintlcg;
0389
0390 }
0391
0392 // **************************************************************************
0393
0394 Int_t hcal::getStackTiming(){
0395   for(int istk = 0; istk<CALSTACKS; istk++){
0396     if(!stacks[istk]) continue;
0397     Int_t stackReject = stacks[istk]->getStackTime();
0398     if(kinds[istk] == CALOR && stackReject>10) rejectEvent ++ ;
0399   }
0400   //  raise reject flag
0401   return rejectEvent;
0402 }
0403
0404 // **************************************************************************
0405
0406 Int_t hcal::collectTrPrimitives(){
0407   int accept = 0;
0408
0409   //  CALOR primitives
0410   for(int istk = 0; istk<CALSTACKS; istk++){
0411     if(!stacks[istk]) continue;
0412     if(!stacks[istk]->stackKind == CALOR) continue;
0413     //  kill LED events at this time
0414     Int_t triggerFlag = stacks[istk]->getTrigger();
0415     accept += triggerFlag*pow(100,stacks[istk]->stackId);
0416   }
0417   return accept;
0418 }
0419
0420
0421 // **************************************************************************
0422 //  done after collectRaw already collected and updated initial values for peaks and times for all objects (towers, counters)
0423 //  It also updates graphical objects in towers so if you need individual towers visible - take care
0424 Int_t hcal::update(Int_t displayMode, Bool_t fitShapes){
0425   //  updates & compute unique amplitudes (in objects with multiple ranges), times and energies
0426   for(int istk = 0; istk<CALSTACKS; istk++){
0427     if(!stacks[istk]) continue;
0428     //  update adds 1000 to whatever is already stored in reject (saturation in towers, below zero threshold). Returns reject
0429     Int_t sr  = stacks[istk]->update(fitShapes);     //  stack reject
0430     if(stacks[istk]->stackKind==CALOR&&sr) rejectEvent += pow(10,stacks[istk]->stackId);
0431   }
0432   //  TODO:  collect global event information for presentation
0433
0434   for(int istk = 0; istk<CALSTACKS; istk++){
0435     if(!alive[istk]||stacks[istk]->stackKind!=CALOR) continue;
0436     //    stacks[istk]->print();
0437     amplTotal += (stacks[istk]->totCorAmpl)*stAScale[istk];
0438     eTotal    += stacks[istk]->E;
0439     calLCG    += (float)(stacks[istk]->stackId)*stacks[istk]->totCorAmpl*stAScale[istk];
0440   }
0441   calLCG /= amplTotal;
0442   totalamp -> Fill(amplTotal);
0443   totale   -> Fill(eTotal);
0444   totallcg -> Fill(calLCG);
0445   //  now scintillators
0446   if(!alive[SCINT]) return 0;
0447   Double_t tcampl(0.);
0448   for (int sc = 0; sc<3; sc++){
0449     scintLCG += (sc+1)*stacks[SCINT]->towers[sc+3]->campl;
0450     tcampl   += stacks[SCINT]->towers[sc+4]->campl;
0451   }
0452   scintLCG /= tcampl;
0453   scintlcg -> Fill(scintLCG);
0454   //
0455   return 0;
0456 }
0457 // **************************************************************************
0458 void hcal::displayRaw(Int_t mode){
0459   if(!mode || mode>2) return;
0460   for(int istk = 0; istk<CALSTACKS; istk++){
0461     if(!stacks[istk]) continue;
0462     //  kill LED
0463     //    if(stacks[istk]->triggerFlag && stacks[istk]->triggerHits<12) {
0464     if(mode==1||(mode==2&&stacks[istk]->stackKind==CALOR)) {
0465       stacks[istk]->displayEvent(mode);
0466     //      if(mode==1)
0467       stacks[istk]->displayADCSum();
0468     }
0469     //  }
0470   }
0471 }
0472
0473 // **************************************************************************
0474 void hcal::displaySummary(Int_t mode){
0475   hcalHelper  * hHelper  = hcalHelper ::getInstance();
0476   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0477   for(int istk = 0; istk<CALSTACKS; istk++){
0478     if(!stacks[istk]) continue;
0479     stacks[istk]->displayTowerSummary(mode);
0480     stacks[istk]->displayStackSummary(mode);
0481   }
0482   TCanvas * calSum;
0483   TString csN = "csSum";
0484   TString csT = "Calorimeter Summary for Run "; csT += hHelper->runnumber;
0485   calSum = (TCanvas *)(gROOT->FindObject(csN));
0486   if(!calSum) calSum  = new TCanvas(csN, csT, 50, 50, 500, 250);
0487   else calSum->Clear();
0488   calSum   -> SetEditable(kTRUE);
0489   calSum   -> Divide(4,1);
0490   calSum   -> cd(1);
0491   totalamp -> Draw();
0492   calSum   -> cd(2);
0493   totale   -> Draw();
0494   calSum   -> cd(3);
0495   totallcg -> Draw();
0496   calSum   -> cd(4);
0497   scintlcg -> Draw();
0498   calSum   -> Update();
0499
0500   TCanvas * satSum;
0501   TString satsN = "satsum";
0502   TString satsT = "Saturation in Calorimeters   Run "; satsT += hHelper->runnumber;
0503   satSum = (TCanvas *)(gROOT->FindObject(satsN));
0504   if(!satSum) satSum  = new TCanvas(satsN, satsT, 50, 50, 600, 400);
0505   else satSum->Clear();
0506   satSum->Divide(2,3);
0507   for(int istk=0; istk<CALSTACKS; istk++)  {
0508     if(!alive[istk]||stacks[istk]->stackKind != CALOR)  continue;
0509     if(stacks[istk]->stackId==EMC) {
0510       satSum -> cd(1);
0511       satSum -> GetPad(1)->SetLogy();
0512       stacks[istk]->satProb[0]->Draw();
0513       satSum -> cd(2);
0514       totale -> Draw();
0515     }
0516     if(stacks[istk]->stackId==HINNER) {
0517       satSum -> cd(3);
0518       satSum -> GetPad(3)->SetLogy();
0519       stacks[istk]->satProb[0]->Draw();
0520       satSum -> cd(4);
0521       satSum -> GetPad(4)->SetLogy();
0522       stacks[istk]->satProb[1]->Draw();
0523     }
0524     if(stacks[istk]->stackId==HOUTER) {
0525       satSum -> cd(5);
0526       satSum -> GetPad(5)->SetLogy();
0527       stacks[istk]->satProb[0]->Draw();
0528       satSum -> cd(6);
0529       satSum -> GetPad(6)->SetLogy();
0530       stacks[istk]->satProb[1]->Draw();
0531     }
0532   }
0533   satSum->Update();
0534   if(!hlHelper->cDirectory.IsNull()){
0535     TString cName  = hlHelper->cDirectory;
0536     cName += "calSum_";  cName += hHelper->runnumber; cName += ".gif";
0537     calSum->SaveAs(cName);
0538     cName  = hlHelper->cDirectory; cName += "satSum_"; cName += hHelper->runnumber; cName += ".gif";
0539     satSum->SaveAs(cName);
0540   }
0541 }
0542 // **************************************************************************
0543
0544 void hcal::clean(){
0545   for(int istk = 0; istk<CALSTACKS; istk++){
0546     if(!stacks[istk]) continue;
0547     stacks[istk]->clean();
0548   }
0549   rejectEvent   = 0;
0550   readyForEvent = true;
0551   amplTotal = 0.;
0552   eTotal    = 0.;
0553   calLCG    = 0.;
0554   scintLCG  = 0.;
0555 }
0556
0557 // **************************************************************************
0558
0559 stack::stack(const Int_t sId, const Int_t twrs, const Int_t xch, const Int_t ych){
0560   reject      =  0;
0561   mapUpdated  =  kFALSE;
0562   eventsSeen  =  0;
0563   //  cout<<"SSS STACK "<<mapUpdated<<endl;
0564   triggerFlag =  0;
0565   triggerHits =  0;
0566   triggerSum  =  0;
0567   stackECalib =  0.;
0568   trAmp       =  NULL;
0569   stackId      = sId;
0570   twrsInStack  = twrs;
0571   stackLoc     = -1;      //  location of this stack in ADC still not defined
0572   nTwrsX       = xch;
0573   nTwrsY       = ych;
0574   gains        = (stackId==HINNER || stackId==HOUTER)?  2  :  1;
0575   towers       = new tower * [twrsInStack];
0576   for (int itw = 0; itw<twrsInStack; itw++){
0577     int xSt = itw%nTwrsX;
0578     int ySt = itw/nTwrsX;
0579     towers[itw] = new tower(stackId, itw, xSt, ySt);
0580     towers[itw]->adcChannel[0]  = itw;
0581     towers[itw]->adcChannel[1]  = (stackId==HINNER || stackId==HOUTER)? itw+HCALTOWERS : -1;
0582   }
0583   //  cout<<"SSS0 "<<mapUpdated<<endl;
0584   for (int ig = 0; ig<gains; ig++){
0585     adcsum[ig]  = new Double_t [NSAMPLES];
0586     fitPar[ig]  = new Double_t [NPARAMETERS];
0587     TString gT  = "Stack ";  gT += stackId; gT += (ig==0? " HIGH " : " LOW "); gT += "gain";
0588     graph[ig]   = new TGraph(NSAMPLES);
0589     graph[ig] -> SetTitle(gT);
0590     TString sum = "adcsum_"; sum += (ig==HIGH)? "HG"  :  "LG";
0591     TString sig = "signal_"; sig += (ig==HIGH)? "HG"  :  "LG";
0592     shape[ig]   = new TF1(sum, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
0593     signal[ig]  = new TF1(sig, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
0594     //    Double_t    times[ACTIVECHANNELS];
0595     //    Double_t    chi2 [ACTIVECHANNELS];
0596   }
0597
0598   //  cout<<"SSS1 "<<mapUpdated<<endl;
0599
0600   E            = 0.;
0601   //  summaries from FIT
0602   // Float_t  maxAmpSum = (stackId==EMC&&emcGainSelection==0)? 2000.  :  1000.;
0603   // Float_t  maxPedSum = (stackId==EMC&&emcGainSelection==0)? 250.   :  25.;
0604   Float_t maxPedSum = 400.;
0605   Float_t maxAmpSum = 16000.;
0606   TString  stfpN = "st_fped"; stfpN += stackId;
0607   TString  stfpT = "Stack ";  stfpT += stackId;  stfpT += " Pedestals - Fit [LG]";
0608   stFPed  = new TH1F(stfpN, stfpT, 100, -maxPedSum, maxPedSum);
0609   TString  stfaN = "st_famp"; stfaN += stackId;
0610   TString  stfaT = "Stack ";  stfaT += stackId;  stfaT += " Amplitudes - Fit [LG]";
0611   stFAmpl = new TH1F(stfaN, stfaT, 550, -0.1*maxAmpSum, maxAmpSum);
0612   TString  stftN = "st_ftime"; stftN += stackId;
0613   TString  stftT = "Stack ";   stftT += stackId; stftT += " Time - Fit [ADC ticks]";
0614   stFTime = new TH1F(stftN, stftT, NSAMPLES, 0., NSAMPLES);
0615   TString  stc2N = "st_c2";   stc2N += stackId;
0616   TString  stc2T = "Stack ";  stc2T += stackId;  stc2T += " Fit Chi2 ";
0617   stChi2  = new TH1F(stc2N, stc2T, 100, 0., 200.);
0618   TString  stgrN = "st_gr";   stgrN += stackId;
0619   TString  stgrT = "Stack ";  stgrT += stackId;  stgrT += " Ratio HG/LG ";
0620   stFGR   = new TH1F(stgrN, stgrT, 100, -20., 20.);
0621
0622   //  cout<<"SSS2 "<<mapUpdated<<endl;
0623   //  summaries from contributing towers
0624   TString  stspN = "st_sped"; stspN += stackId;
0625   TString  stspT = "Stack ";  stspT += stackId;  stspT += " Stack Pedestals Sum [LG]";
0626   stSPed  = new TH1F(stspN, stspT, 100, -maxPedSum, maxPedSum);
0627   TString  stsaN = "st_samp"; stsaN += stackId;
0628   TString  stsaT = "Stack ";  stsaT += stackId;  stsaT += " Stack Amplitude Sum  [LG]";
0629   stSAmpl = new TH1F(stsaN, stsaT, 550, -0.1*maxAmpSum, maxAmpSum);
0630   TString  stseN = "st_sen";  stseN += stackId;
0631   TString  stseT = "Stack ";  stseT += stackId; stseT += " Stack Energy Sum [MeV]";
0632   stSE    = new TH1F(stseN, stseT, 1100, -5., 50.);
0633   TString  ststN = "st_stime"; ststN += stackId;
0634   TString  ststT = "Stack ";   ststT += stackId; ststT += " <Stack Time> [clock ticks]";
0635   stAvT   = new TH1F(ststN, ststT, NSAMPLES, 0., NSAMPLES);
0636   TString  ststsN = "st_stsig"; ststsN += stackId;
0637   TString  ststsT = "Stack ";   ststsT += stackId; ststsT += " Time - RMS [clock ticks]";
0638   stTRMS  = new TH1F(ststsN, ststsT, NSAMPLES, 0., NSAMPLES);
0639   //  Hit impacts in the stack
0640   TString hitCGN   = "hitSt";  hitCGN += stackId;
0641   TString hitCGT   = "Stack "; hitCGT += stackId; hitCGT += "  Impact X vs Y [towers]";
0642   hitCG    = new TH2F(hitCGN, hitCGT, nTwrsX, 0, nTwrsX, nTwrsY, 0, nTwrsY);
0643
0644   //  SATURATION
0645   for(int ig=0; ig<gains; ig++) {
0646     TString satN = "stsat";  satN += stackId; satN += ig;
0647     TString satT = "Stack "; satT += stackId; satT += ig==0? " HG " : " LG "; satT += "saturation probability vs # of sat. samples";
0648     satProb[ig] = new TH1F(satN, satT, 25, 0., 25.);
0649   }
0650
0651
0652   if(stackKind != CALOR/*||hlHelper->runKind!="cosmics"*/) return;
0653   //  cosmic related summaries
0654   TString trAmpN = "trAmp";    trAmpN += stackId;
0655   TString trAmpT = "Stack ";   trAmpT += stackId; trAmpT += "  Track Amplitude";
0656   trAmp    = new TH1F(trAmpN,    trAmpT,   300,  0., 300.);
0657   TString trCh2N = "trCh2";    trCh2N += stackId;
0658   TString trCh2T = "Stack ";   trCh2T += stackId; trCh2T += "  Track Chi2";
0659   trChi2   = new TH1F(trCh2N,    trCh2T,   100,  0., 10.);
0660   TString trCrN = "trCr";     trCrN   += stackId;
0661   TString trCrT = "Stack ";   trCrT   += stackId; trAmpT += "  Track Crossing";
0662   trCr     = new TH1F(trCrN,     trCrT,     40,  0., 4.);
0663   TString trSlN = "trSl";     trSlN   += stackId;
0664   TString trSlT = "Stack ";   trSlT   += stackId; trAmpT += "  Track Slope";
0665   trSl     = new TH1F(trSlN,   trSlT,     50, -1., 1.);
0666   TString trACh2N = "trACh2";  trACh2N += stackId;
0667   TString trACh2T = "Stack ";  trACh2T += stackId; trACh2T += "  Track Amplitude vs Chi2";
0668   trAmpCh2 = new TH2F(trACh2N, trACh2T, 300, 0,  300, 100, 0., 2.);
0669   TString trAstHN = "trAstH";  trAstHN += stackId;
0670   TString trAstHT = "Stack ";  trAstHT += stackId; trAstHT += "  Track Amplitude vs # of Hits in Stack";
0671   trAstH   = new TH2F(trAstHN, trAstHT,300, 0,  300, 100, 0., 16.);
0672   TString trASlN = "trASl";    trASlN += stackId;
0673   TString trASlT = "Stack ";   trASlT += stackId; trASlT += "  Track Amplitude vs # Track Slope";
0674   trASl    = new TH2F(trASlN,  trASlT, 300, 0,  300, 100, 0., 1.);
0675   //  Next scatterplot is the basis of calibration
0676   //  Its data are HG data divided by H/L gain ratio and so to make this
0677   //  picture representative the number of bins is chose equal to
0678   //  amplimit*hgratio[stackId]
0679   //  cout<<"SSS4 "<<mapUpdated<<endl;
0680   Float_t amplimit = 100;  //  LG ADC counts
0681   Int_t   nx       = amplimit*hlgratios[stackId];
0682   TString trATwrN = "trATwrN"; trATwrN += stackId;
0683   TString trATwrT = "Stack ";  trATwrT += stackId; trATwrT += "  Tagged Twr Amplitude vs Twr # ";
0684   trATwr   = new TH2F(trATwrN,  trATwrT, nx, 0, 100, twrsInStack, 0., twrsInStack);
0685   //  cout<<"SSSExit "<<mapUpdated<<endl;
0686
0687 }
0688
0689 // **************************************************************************
0690 stack::~stack(){
0691   if(stackId < 0 || !towers) return;
0692   for (int itw = 0; itw < twrsInStack; itw++){
0693     if(towers[itw])  delete towers[itw];
0694   }
0695   delete stFGR;
0696   delete stFPed;
0697   delete stFAmpl;
0698   delete stFTime;
0699   delete stChi2;
0700   delete stSPed;
0701   delete stSAmpl;
0702   delete stSE;
0703   delete stAvT;
0704   delete stTRMS;
0705   delete hitCG;
0706   if(trAmp) {
0707     delete trAmp;
0708     delete trChi2;
0709     delete trCr;
0710     delete trSl;
0711     delete trAmpCh2;
0712     delete trAstH;
0713     delete trASl;
0714     delete trATwr;
0715   }
0716
0717 }
0718
0719 // **************************************************************************
0720
0721 void stack::displayEvent(Int_t mode){
0722   if(!mode || mode>2)  return;
0723   hcalHelper * hHelper = hcalHelper::getInstance();
0724   TCanvas * pattern[2];
0725   for(int ig = 0; ig < gains; ig++){
0726     TString stckN = "StPat";   stckN += stackId; stckN += "_G_";    stckN += ig;
0727     TString stckT = "Stack ";  stckT += stackId; stckT += (ig==HIGH)? " HGain" : " LGain"; stckT += "  Event Pattern Display  Run "; stckT += hHelper->runnumber;
0728     pattern[ig] = (TCanvas *)(gROOT->FindObject(stckN));
0729     if(!pattern[ig]) pattern[ig]  = new TCanvas(stckN, stckT, 50*(stackId+ig), 50*(stackId+ig), 100*nTwrsX, 100*nTwrsY);
0730     else pattern[ig]->Clear();
0731     pattern[ig] -> SetEditable(kTRUE);
0732     //  display digitizations
0733     if(mode == 1 || stackKind != CALOR) {
0734       pattern[ig] -> Divide(nTwrsX,nTwrsY);
0735       for (int itw = 0; itw<twrsInStack; itw++){
0736     Int_t ipad = (nTwrsY-(towers[itw]->y)-1)*nTwrsX + (towers[itw]->x)+1;
0737     pattern[ig] -> cd(ipad);
0738     //      cout<<"<stack::displayPattern "<<nTwrsX<<" "<<nTwrsY<<"  "<<x<<" "<<y<<"  "<<ipad<<endl;
0739     towers[itw]->graph[ig]->SetMarkerStyle(20);
0740     towers[itw]->graph[ig]->SetMarkerSize(0.4);
0741     towers[itw]->graph[ig]->Draw("acp");
0742     //  towers[itw]->shape[ig]->Draw("same");
0743       }
0744       pattern[ig]->Update();
0745     } else {
0746       TH2 * stackLego;
0747       //  display stack lego plot - only for calorimeter stacks
0748       TString slN  = "st"; slN += stackId; slN += "_G_"; slN += ig;
0749       TString slT  = "Mixed Amplitudes [LG] in Stack "; slT += stackId;
0750       if(!(stackLego = (TH2F*) (gROOT->FindObject(slN)))) stackLego = new TH2F(slN, slT, nTwrsX, 0, nTwrsX, nTwrsY, 0, nTwrsY);
0751       for (int itw = 0; itw<twrsInStack; itw++){
0752     //  towers[itw]->print();
0753     //  cout<<"stack "<<stackId<<"  itw/x/y/E  "<<itw<<" / "<<towers[itw]->x<<" / "<< towers[itw]->y<<" / "<< towers[itw]->campl<<endl;
0754     stackLego->SetBinContent(1+towers[itw]->x, 1+towers[itw]->y, towers[itw]->campl);
0755       }
0756       stackLego  -> Draw("lego");
0757       pattern[ig]->Update();
0758       if(ig==0) break;
0759    }
0760   }
0761 }
0762
0763 // **************************************************************************
0764
0765 void stack::displayTowerSummary(Int_t mode){
0766   if(mode<=2)  return;
0767   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0768   hcalHelper  * hHelper = hcalHelper ::getInstance();
0769   TCanvas * corTwrsA;
0770   TString cTAN = "rtA_st";  cTAN += stackId;
0771   TString cTAT = "Stack ";  cTAT += stackId;  cTAT += " Hit Amplitudes in Towers [LG units]  Run ";  cTAT += hHelper->runnumber;
0772   corTwrsA = (TCanvas *)(gROOT->FindObject(cTAN));
0773   if(!corTwrsA) corTwrsA  = new TCanvas(cTAN, cTAT, 50*stackId, 50*stackId, 100*nTwrsX, 100*nTwrsY);
0774     else corTwrsA->Clear();
0775   corTwrsA -> SetEditable(kTRUE);
0776   corTwrsA -> Divide(nTwrsX,nTwrsY);
0777   for (int itw = 0; itw<twrsInStack; itw++){
0778     Int_t ipad = (nTwrsY-(towers[itw]->y)-1)*nTwrsX + (towers[itw]->x)+1;
0779     corTwrsA -> cd(ipad);
0780     corTwrsA -> GetPad(ipad)->SetLogy();
0781     towers[itw]->twrAmpl->SetMarkerStyle(20);
0782     towers[itw]->twrAmpl->SetMarkerSize(0.4);
0783     towers[itw]->twrAmpl->Draw();
0784     // rawTwrsT -> cd(ipad);
0785     // //    rawTwrsA -> GetPad(ipad)->SetLogy();
0786     // towers[itw]->twrTime->SetMarkerStyle(20);
0787     // towers[itw]->twrTime->SetMarkerSize(0.4);
0788     // towers[itw]->twrTime->Draw();
0789   }
0790   corTwrsA->Update();
0791   if(!hlHelper->cDirectory.IsNull()){
0792     TString cName  = hlHelper->cDirectory;
0793     cName += "corTwrsA_";  cName += hHelper->runnumber;   cName += ".gif";
0794     corTwrsA->SaveAs(cName);
0795   }
0796
0797   //  rawTwrsT->Update();
0798 }
0799 // **************************************************************************
0800
0801 void stack::displayStackSummary(Int_t mode){
0802   if(mode != 3 || stackKind != CALOR) return;
0803   hcalHelper  * hHelper  = hcalHelper ::getInstance();
0804   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0805   TCanvas * stackSum(NULL), * muTracks(NULL);
0806   TString stN = "stSum";  stN += stackId;
0807   TString stT = "Stack "; stT += stackId;  stT += "  Summary for Run "; stT += hHelper->runnumber;
0808   stackSum = (TCanvas *)(gROOT->FindObject(stN));
0809   if(!stackSum) stackSum  = new TCanvas(stN, stT, 50*stackId, 50*stackId, 400, 400);
0810     else stackSum->Clear();
0811   stackSum -> SetEditable(kTRUE);
0812   stackSum -> Divide(4,2);
0813   //  TH1        * stFPed, * stFAmpl, * stFTime, * stChi2;
0814   //  TH1        * stSPed, * stSAmpl, * stSE,    * stAvT,   * stTRMS;
0815   stackSum ->cd(5);
0816   stSPed   ->Draw();
0817   stackSum ->cd(6);
0818   stackSum -> GetPad(6)->SetLogy();
0819   stSAmpl  ->Draw();
0820   stackSum ->cd(7);
0821   stAvT    ->Draw();
0822   stackSum ->cd(8);
0823   stackSum -> GetPad(8)->SetLogy();
0824   stSE     ->Draw();
0825
0826   stackSum ->cd(1);
0827   stFPed   ->Draw();
0828   stackSum ->cd(2);
0829   stackSum -> GetPad(2)->SetLogy();
0830   stFAmpl  ->Draw();
0831   stackSum ->cd(3);
0832   stFGR    ->Draw();
0833   //  stFTime  ->Draw();
0834   stackSum ->cd(4);
0835   stChi2   ->Draw();
0836   stackSum ->Update();
0837   //  saturation probability
0838   if(eventsSeen) {
0839     for (int ig=0; ig<gains; ig++) {
0840       satProb[ig]->Scale(1./(float)eventsSeen);
0841       satProb[ig]->Scale(1./(float)twrsInStack);
0842     }
0843   }
0844
0845   if(!hlHelper->cDirectory.IsNull()){
0846     TString cName  = hlHelper->cDirectory;
0847     cName += "stackSum_"; cName += stackId; cName += "_"; cName += hHelper->runnumber;  cName += ".gif";
0848     stackSum->SaveAs(cName);
0849   }
0850   if(hlHelper->runKind == "cosmics"){
0851     TString muN = "muTr";   muN += stackId;
0852     TString muT = "Stack "; muT += stackId;  muT += "  Cosmic Muons in Run "; muT += hHelper->runnumber;
0853     muTracks = (TCanvas *)(gROOT->FindObject(muN));
0854     if(!muTracks) muTracks  = new TCanvas(muN, muT, 50*stackId, 50*stackId, 400, 400);
0855     else muTracks->Clear();
0856     muTracks -> SetEditable(kTRUE);
0857     muTracks -> Divide(3,2);
0858     muTracks -> cd(1);
0859     trAmp    -> Draw();
0860     muTracks -> cd(2);
0861     trCr     -> Draw();
0862     muTracks -> cd(3);
0863     trSl     -> Draw();
0864     muTracks -> cd(4);
0865     trASl    -> Draw();
0866     //   trChi2   -> Draw();
0867     muTracks -> cd(5);
0868     trAmpCh2 -> Draw();
0869     muTracks -> cd(6);
0870     trAstH   -> Draw();
0871     muTracks -> Update();
0872
0873     if(!hlHelper->cDirectory.IsNull()){
0874       TString cName  = hlHelper->cDirectory;
0875       cName += "muTracks_"; cName += stackId; cName += "_"; cName += hHelper->runnumber;  cName += ".gif";
0876       muTracks->SaveAs(cName);
0877     }
0878   }
0879 }
0880
0881 // **************************************************************************
0882
0883 void stack::displayADCSum(){
0884   hcalHelper  * hHelper  = hcalHelper ::getInstance();
0885   hLabHelper  * hlHelper = hLabHelper ::getInstance();
0886   TString stckN = "StEvS";    stckN += stackId;
0887   TString stckT = "Stack "; stckT += stackId; stckT += "  Stack Event Sum Display  Run "; stckT += hHelper->runnumber;
0888   TCanvas * stckEvSumDisplay = (TCanvas *)(gROOT->FindObject(stckN));
0889   if(!stckEvSumDisplay) stckEvSumDisplay  = new TCanvas(stckN, stckT, 25*stackId, 25*stackId, 600, 400);
0890   else stckEvSumDisplay->Clear();
0891   stckEvSumDisplay -> Divide(1,gains);
0892   for(int ig = 0; ig < gains; ig++){
0893     stckEvSumDisplay->cd(ig+1);
0894     //    for (int is = 0; is < NSAMPLES; is++)  evtGraph[ig]->SetPoint(is,(Double_t)is,evtadcsum[ig][is]);
0895     //  cout<<"<tileDisplay>   call to fitTileSignal"<<endl;
0896     //  eventSum->Divide(1,2);
0897     //  eventSum->cd(1);
0898     //  eventSum->SetEditable(kTRUE);
0899     graph[ig]->SetMarkerStyle(20);
0900     graph[ig]->SetMarkerSize(0.4);
0901     graph[ig]->Draw("acp");
0902     //  eventSum->cd(2);
0903     shape[ig]->Draw("same");
0904     // evtShape[ig]->Draw();
0905     // evtGraph[ig]->Draw("same");
0906     //    stckDisplay->Update();
0907   }
0908   stckEvSumDisplay->Update();
0909   if(!hlHelper->cDirectory.IsNull()){
0910     TString cName  = hlHelper->cDirectory;
0911     cName += "stackEvSum_"; cName += stackId; cName += "_"; cName += hHelper->runnumber;  cName += ".gif";
0912     stckEvSumDisplay->SaveAs(cName);
0913   }
0914 }
0915 // **************************************************************************
0916
0917 //   clean stack for a new event
0918 void stack::print(){
0919   cout<<"<stack::print>  stackId "<<stackId<<"  Gains "<<gains<<"  twrsInStack "<<twrsInStack<<"  stackLoc "<<stackLoc<<"  TwrsXY "<<nTwrsX<<"/"<<nTwrsY<<"  mapUpdated "<<mapUpdated<<endl;
0920   for (int itw = 0; itw < twrsInStack; itw++){
0921     if(towers[itw])  {
0922       towers[itw]->print();
0923     }
0924   }
0925 }
0926
0927
0928
0929 // **************************************************************************
0930 void stack::updateMap(Int_t stckloc){
0931   //  print();
0932   if(mapUpdated){
0933     cout<<"<stack::updateMap>  Already updated - EXIT"<<endl;
0934     return;
0935   }
0936   stackLoc     = stckloc;
0937   for (int itw = 0; itw<twrsInStack; itw++){
0938     towers[itw]->adcChannel[0] += stckloc;
0939     if(gains==2) towers[itw]->adcChannel[1] += stckloc;
0940   }
0941   mapUpdated = kTRUE;
0942 }
0943
0944
0945 // **************************************************************************
0946
0947 //   collectdata from towers in stack for a new event
0948 Int_t stack::update(Bool_t fitShape){
0949   eventsSeen ++;
0950   Int_t  used = 0;
0951   for (int itw = 0; itw < twrsInStack; itw++){
0952     if(towers[itw])  {
0953       Int_t rejected = towers[itw]->update(fitShape);
0954       if(rejected<=2){
0955     used++;
0956     totPed     += towers[itw]->cped;
0957     totAmpl    += towers[itw]->rampl;
0958     totCorAmpl += towers[itw]->campl;
0959     //  with fit time used in every tower to identify event time
0960     //  all towers in the stack will now have identical times
0961     avTwrTime  += towers[itw]->ctime;
0962     rmsTwrTime += towers[itw]->ctime*towers[itw]->ctime;
0963     E          += towers[itw]->E;
0964       } //  else  if(rejected>5) towers[itw]->display();
0965       for(int ig = 0; ig<gains; ig++){
0966     satProb[ig] -> Fill(towers[itw]->satFlag[ig]);
0967       }
0968     }
0969   }
0970   if(!used) {
0971     //  everything in saturation
0972     reject += 1000;
0973     return reject;
0974   }
0975   avTwrTime /= used;
0976   rmsTwrTime = 0.;    // (twrsUsed==0? 0.  :  sqrt(rmsTwrTime/twrsUsed-avTwrTime*avTwrTime));
0977   stSE    -> Fill(E);
0978   stSPed  -> Fill(totPed);
0979   stSAmpl -> Fill(totCorAmpl);
0980   stAvT   -> Fill(avTwrTime);
0981   stTRMS  -> Fill(rmsTwrTime);
0982
0983   //  Compute impact CG and find cosmic track
0984   if(stackKind == CALOR) {
0985     getStackImpact();
0986     hLabHelper  * hlHelper = hLabHelper ::getInstance();
0987     if(hlHelper->runKind=="cosmics") getStackTrack();
0988   }
0989   return reject;
0990 }
0991
0992
0993 // **************************************************************************
0994
0995 //   getStackTime  creates summary event (one per gain)
0996 Int_t stack::getStackTime(){
0997   hcalHelper * hHelper = hcalHelper::getInstance();
0998   for(int ig = 0; ig < gains; ig++){
0999     //  build eventsum (adc sum without pedestal subtruction)
1000     //  do fitting etc (HG/LG separately)
1001     twrsUsed[ig] = 0.;
1002     for (int itw = 0; itw < twrsInStack; itw++) {
1003       if(towers[itw]->satFlag[ig]>2) continue;
1004       twrsUsed[ig] ++;
1005       for(int is = 0; is < NSAMPLES; is++)
1006     adcsum[ig][is] += hHelper->adc[towers[itw]->adcChannel[ig]][is];
1007     }
1008     if(twrsUsed[ig]>=1) {
1009       for (int is = 0; is < NSAMPLES; is++)
1010     graph[ig]->SetPoint(is,(Double_t)is,adcsum[ig][is]);
1011     } else {
1012       reject += 10;
1013       continue;
1014     }
1015
1016     Double_t rVal, rTime;
1017     Int_t iss(RISETIME), ise(NSAMPLES-FALLTIME);
1018     Int_t rMaxVal  = adcsum[ig][iss]; Int_t maxsmpl = iss;
1019     Int_t rMinVal  = adcsum[ig][iss]; Int_t minsmpl = iss;
1020
1021     for (int is = iss; is<=ise; is++) {
1022       if(adcsum[ig][is]>rMaxVal) {
1023     rMaxVal  = adcsum[ig][is];
1024     maxsmpl  = is;
1025       }
1026       if(adcsum[ig][is]<rMinVal) {
1027     rMinVal  = adcsum[ig][is];
1028     minsmpl  = is;
1029       }
1030     }
1031     if(abs(rMinVal) > abs(rMaxVal)) {
1032       rVal = rMinVal;  rTime = minsmpl;
1033     } else {
1034       rVal = rMaxVal;  rTime = maxsmpl;
1035     }
1036     //    cout<<"<stack::getStackTime>:  Raw Signal at "<<rTime<<"  Raw Amplitude "<<rVal<<endl;
1037
1038     Double_t  par0[6];
1039     par0[0] = rMinVal;   //rVal;      //3.;
1040     par0[1] = minsmpl;   //rTime;  //max(0.,(Double_t)(rTime-RISETIME));
1041     par0[2] = 4.;
1042     par0[3] = 1.6;
1043     par0[4] = 0.; //PEDESTAL*twrsInStack;  //   we do not do pedestal subtrastion at this time
1044     par0[5] = 0.;                     //   slope of the pedestals is initially set to "0"
1045     shape[ig]->SetParameters(par0);
1046     for(int parn=0; parn<4; parn++) shape[ig]->SetParLimits(parn, par0Min[parn], par0Max[parn]);
1047     shape[ig]->SetParLimits(1, minsmpl-RISETIME, minsmpl+FALLTIME);// rTime-RISETIME, rTime+FALLTIME/2.);
1048     shape[ig]->SetParLimits(4, -256., 256.);
1049     shape[ig]->SetParLimits(5, -0.5,  0.5);
1050     graph[ig]->Fit(shape[ig], "RQWM", "NQ", 0., (Double_t)NSAMPLES);
1051     shape[ig]->GetParameters(fitPar[ig]);
1052     //   cout<<"<stack::getStackTime>: Gain "<<ig<<" minsmpl "<<rTime<<"  Parameters (0) "<<fitPar[ig][0]<<" (1) "<<fitPar[ig][1]<<" (2) "<<fitPar[ig][2]<<" (3) "<<fitPar[ig][3]<<" (4) "<<fitPar[ig][4]<<" (5) "<<fitPar[ig][5]<<endl;
1053     //  Signal is whatever left after fitted pedestal subtraction (removal)
1054     signal[ig]->SetParameters(fitPar[ig]);
1055     signal[ig]->SetParameter(4, 0.);
1056     signal[ig]->SetParameter(5, 0.);
1057     // Double_t sgp  = signal[ig]->GetMaximum(par0Min[1], par0Max[1]);
1058     // Double_t tmp  = signal[ig]->GetMaximumX(par0Min[1], par0Max[1]);
1059     Double_t sgn  = signal[ig]->GetMinimum(par0Min[1], par0Max[1]);
1060     Double_t tmn  = signal[ig]->GetMinimumX(par0Min[1], par0Max[1]);
1061     fitPeak[ig]  = sgn;
1062     fitTime[ig]  = tmn;
1063     fitPed[ig]  = (fitPar[ig][4]+fitPar[ig][5]*fitTime[ig]);
1064     fitChi2[ig] = shape[ig]->GetChisquare()/shape[ig]->GetNDF();
1065     //    cout<<"<stack::getStackTime>: Gain "<<ig<<"  Peak "<<  fitPeak[ig]<<"  Time "<< fitTime[ig]<<"  Chi2 "<<fitChi2[ig]<<endl;
1066   }  //  end of loop over gain ranges
1067
1068   //  select gain range to use for Summaries
1069   //  TODO - study chi2 vs ????? - can it really signal overflow
1070   gainToUse = 0;
1071   Float_t  gratio(0.);
1072   Float_t  scale = 1.;
1073   if(stackKind == CALOR) {
1074     if(reject==20||(stackId==EMC&&reject==10)) return reject;  // saturation
1075     if(stackId==EMC) {
1076       gainToUse      = 0;
1077       if(emcGainSelection == HIGH) scale = HLGRATIO;
1078     } else {
1079       if(!reject) {
1080     gainToUse = 0;
1081     gratio = fitPeak[0]/fitPeak[1];
1082     scale  = hlgratios[stackId];
1083       } else {
1084     gainToUse = 1;
1085       }
1086     }
1087     if(twrsUsed[gainToUse]==twrsInStack){
1088       Double_t amp    = -fitPeak[gainToUse]/scale;
1089       if(amp < STZEROSUPTHR) reject += 100;  //  stack empty
1090       if(stackId!=EMC&&gainToUse==0)  stFGR  -> Fill(gratio);
1091       stFPed -> Fill(fitPed[gainToUse]/scale);
1092       stFAmpl-> Fill(amp);
1093       stFTime-> Fill(fitTime[gainToUse]);
1094       stChi2 -> Fill(fitChi2[gainToUse]/scale);
1095     }
1096   }
1097   return reject;
1098 }
1099
1100 // **************************************************************************
1101
1102 //   compute impact positions in units of towers
1103 //   return
1104 Int_t  stack::getStackImpact(){
1105   if(stackKind!=CALOR) return twrsInStack;
1106   xImpact = -1.;
1107   yImpact = -1.;
1108   Int_t used(0);
1109   double sxw(0.), syw(0.), sw(0.);
1110   for (int itw=0; itw < twrsInStack; itw++) {
1111     if(towers[itw]->satFlag[towers[itw]->gainToUse]<=2){
1112       used++;
1113       sw   += towers[itw]->campl;
1114       sxw  += towers[itw]->campl*(Float_t)towers[itw]->x;
1115       syw  += towers[itw]->campl*(Float_t)towers[itw]->y;
1116     }
1117   }
1118   if(used&&sw>0.) {
1119     xImpact = sxw/sw;
1120     yImpact = syw/sw;
1121   }
1122   //  cout<<"<stack::getStackImpact> Stack "<<stackId<<"  sw "<<sw<<" X/Y "<<xImpact<<"/"<<yImpact<<endl;
1123   hitCG -> Fill(xImpact, yImpact);
1124   return twrsInStack - used;
1125 }
1126
1127 // **************************************************************************
1128
1129 //   clean stack for a new eventcompute impact positions in units of towers
1130 //   return kFALSE if
1131 Int_t  stack::getStackTrack(){
1132   if(stackKind!=CALOR) return 0;
1133   double ixs(0.), iys(0.), ix2s(0.), iy2s(0.), ixys(0.);
1134   double crossing(0.), slope(0.);
1135   Int_t  hits[nTwrsY];
1136   //  tag maximum amplitudes in each row
1137   //  it is possibly better to compute CG's in every row and use those
1138   //  for tracking muons
1139   for(int r =0; r<nTwrsY; r++){   //  loop over raws
1140     hits[r]  = 0;   Double_t hitAmpl = 0.;
1141     for (int c = 0; c<nTwrsX; c++) {
1142       int itw = nTwrsX*r+c;
1143       if(towers[itw]->campl>hitAmpl) {hits[r] = itw;  hitAmpl = towers[itw]->campl;}
1144       if(towers[itw]->campl>TWRAMPTHR&&towers[itw]->satFlag[towers[itw]->gainToUse]<=2) stackHits ++;
1145     }
1146   }
1147   for (int r =0; r<nTwrsY; r++){    //  loop over towers peaking in rows
1148     Double_t amp  = towers[hits[r]]->campl;
1149     if(amp>TWRAMPTHR) {
1150       trackHits ++;
1151       trackAmpl +=  amp;
1152       Float_t  x = (Float_t)towers[hits[r]]->x;
1153       Float_t  y = (Float_t)towers[hits[r]]->y;
1154       ixs       += (amp*x);
1155       iys       += (amp*y);
1156       ix2s      += (amp*x*x);
1157       iy2s      += (amp*y*y);
1158       ixys      += (amp*x*y);
1159     }
1160   }
1161   if(trackHits<=2) return trackHits;
1162   //  we solve for x=F(y) - cosmics is vertical
1163   if(trackAmpl>0.) {
1164     //      icr    = (ixs*ixys - ix2s*iys)/(ixs*ixs-is*ix2s);
1165     //      islope = (-is*ixys+ixs*iys)   /(ixs*ixs-is*ix2s);
1166     crossing  = (iys*ixys - iy2s*ixs)/(iys*iys-trackAmpl*iy2s);
1167     slope     = (-trackAmpl*ixys+iys*ixs)   /(iys*iys-trackAmpl*iy2s);
1168   }
1169
1170   //  convert back to y = F(x)
1171
1172   //  now estimator
1173   for (int itw = 0; itw < twrsInStack; itw++){
1174     Double_t amp  = towers[itw]->campl;;
1175     Float_t  x = (Float_t)towers[itw]->x;
1176     Float_t  y = (Float_t)towers[itw]->y;
1177     float  est(0.), estx(0.), esty(0.);
1178     if(slope!=0.) {
1179       est     = abs((y - x/slope+crossing/slope)/sqrt(pow(1./slope,2)+1.));
1180       //      estx    = ((x-(islope*y+icr))*(x-(islope*y+icr)));
1181       //      esty    = islope!=0? (y-(x-icr)/islope)*(y-(x-icr)/islope)  : 0.;
1182     } else {
1183       //  crossing and slope are just tower coordinate in this case
1184       estx = (x-crossing);
1185       esty = (y-slope);
1186       est  = sqrt(estx*estx+esty*esty)/2;
1187     }
1188     trchi2 += est*amp;
1189   }
1190   trchi2 /= totAmpl;
1191   trAmp   -> Fill(trackAmpl);
1192   trChi2  -> Fill(trchi2);
1193   trCr    -> Fill(crossing);
1194   trSl    -> Fill(slope);
1195   trAmpCh2-> Fill(trackAmpl, trchi2);
1196   trAstH  -> Fill(trackAmpl, stackHits);
1197   trASl   -> Fill(trackAmpl, abs(slope));
1198   if(abs(slope)<0.45||abs(slope)>0.55) {
1199     for (int itw = 0; itw < twrsInStack; itw++){
1200       if(towers[itw]->campl>TWRAMPTHR)
1201     trATwr->Fill(towers[itw]->campl, (float)itw);
1202     }
1203   }
1204   //  cout<<"<stack::getStackTrack> Stack "<<stackId<<"  Hits "<<stackHits<<"/"<<trackHits<<"  Cr/Sl "<<crossing<<"/"<<slope<<"  Chi2 "<<trchi2<<"  Tr Ampl "<<trackAmpl<<endl;
1205   return trackHits;
1206 }
1207
1208 // **************************************************************************
1209
1210 //   clean stack for a new event
1211 Int_t  stack::getTrigger(){
1212
1213   for (int itw = 0; itw < twrsInStack; itw++) {
1214     Float_t twrAmp = (gains==1)?  -(towers[itw]->rawAmpl[0])  :  -(towers[itw]->rawAmpl[TRGAINRANGE]);
1215     if((int)(twrAmp)/8>TWRAMPTHR) triggerHits ++;
1216     triggerSum += twrAmp;
1217   }
1218   if((int)triggerSum/8>STTOTAMPTHR)                             triggerFlag += 1;
1219   if(triggerHits >= STHITMIN && triggerHits <=STHITMAX)         triggerFlag += 10;
1220   return triggerFlag;
1221 }
1222 // **************************************************************************
1223
1224 //   clean stack for a new event
1225 void stack::clean(){
1226   for (int itw = 0; itw < twrsInStack; itw++) {  if(towers[itw])  towers[itw]->clean();}
1227   for (int ig = 0; ig<gains; ig++){
1228     for (int is=0; is<NSAMPLES; is++) {
1229       adcsum[ig][is] = 0.;
1230     }
1231     for(int ip=0; ip<NPARAMETERS; ip++){
1232       fitPar[ig][ip] = 0.;
1233     }
1234     fitPed[ig]  = 0.;
1235     fitPeak[ig] = 0.;
1236     fitChi2[ig] = 0.;
1237   }
1238   reject      = 0;
1239   totPed      = 0.;
1240   totAmpl     = 0.;
1241   totCorAmpl  = 0.;
1242   avTwrTime   = 0.;
1243   rmsTwrTime  = 0.;
1244   twrsUsed[0] = 0;
1245   twrsUsed[1] = 0;
1246   triggerHits = 0;
1247   triggerSum  = 0.;
1248   triggerFlag = 0;
1249   stackHits   = 0;
1250   trackHits   = 0;
1251   crossing    = 0.;
1252   slope       = 0.;
1253   trchi2      = 0.;
1254   trackAmpl   = 0.;
1255   E           = 0.;
1256 }
1257
1258
1259 // ******************************** TOWER ***********************************
1260
1261 tower::tower(Int_t stk, Int_t ind, Int_t xSt, Int_t ySt){
1262   gainToUse =  0;
1263   reject    =  kFALSE;
1264   attached  =  kFALSE;
1265   stackId   = stk;
1266   index     = ind;
1267   x         = xSt;
1268   y         = ySt;
1269   twrECalib = 1.;
1270   twrEScale = 1.;
1271   gains        = (stackId==HINNER || stackId==HOUTER)?  2  :  1;
1272   if(gains==1) gainToUse = 0;
1273   for(int ig = 0; ig < gains; ig++) {
1274     graph[ig]   = new TGraph(NSAMPLES);
1275     TString title = "Stack ";     title += stackId; title += " Tower "; title += index;
1276     title += "(";     title += x; title += "/";     title += y;         title += ") ";
1277     title += (ig==HIGH)? " HGain"  :  " LGain";
1278     TString sig = "sig_"; sig += (ig==HIGH)? "HG_"  :  "LG_";
1279     sig += stackId;   sig   += "_";                sig   += index;
1280     TString shp = "shp_"; shp += (ig==HIGH)? "HG_"  :  "LG_";
1281     shp += stackId;   shp   += "_";                shp   += index;
1282     graph[ig]->SetTitle(title);
1283     shape[ig]   = new TF1(shp, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
1284     signal[ig]  = new TF1(sig, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
1285   }
1286   //  run summaries
1287   TString  cpN = "twr_ped"; cpN += stackId; cpN += index;
1288   TString  cpT = "Stack "; cpT += stackId; cpT += " Tower "; cpT += index; cpT += " Mixed pedestals [LG]";
1289   twrPed  = new TH1F(cpN, cpT, 100, -5., 5.);
1290   TString  caN = "twr_amp"; caN += stackId; caN += index;
1291   TString  caT = "Stack "; caT += stackId; caT += " Tower "; caT += index; caT += " Mixed Amplitudes [LG]";
1292   twrAmpl = new TH1F(caN, caT, 250, -50., 200.);
1293   TString  ctN = "twr_time"; ctN += stackId; ctN += index;
1294   TString  ctT = "Stack "; ctT += stackId; ctT += " Tower "; ctT += index; ctT += " Time [ADC ticks]";
1295   twrTime = new TH1F(ctN, ctT, NSAMPLES, 0., NSAMPLES);
1296   TString  ceN = "twr_en"; ceN += stackId; ceN += index;
1297   TString  ceT = "Stack "; ceT += stackId; ceT += " Tower "; ceT += index; ceT += " Energy [MeV]";
1298   twrE    = new TH1F(ceN, ceT, 250, -500., 20000.);
1299
1300   clean();
1301 }
1302
1303 // **************************************************************************
1304
1305 //   update event information from hcalHelper
1306 //   must be called after stackTiming (Bool_t fit) and collectRaw()
1307 Int_t tower::update(Bool_t fitShape){
1308   //  hcalHelper * hHelper = hcalHelper::getInstance();
1309
1310   //  decide on gainToUse for this tower
1311   gainToUse = 0;
1312   if(gains==2&&satFlag[0]>0)  gainToUse = 1;     //  all gain ranges saturated
1313
1314   //  compute unique values (bridge the gain ranges). Keep raw amplitudes as they came from collectRaw (signed and unscaled - we still need those for trigger emulator
1315   rped  = rawPed[gainToUse];
1316   rampl = -rawAmpl[gainToUse];
1317   if(gains>1&&gainToUse==0) {
1318     rped  /= hlgratios[stackId];
1319     rampl /= hlgratios[stackId];
1320   }
1321   rtime = rawTime[gainToUse];
1322   cped  = rped; campl = rampl; ctime = rtime;   //  no corrections applied yet
1323   //  if this is EMC and it has emcGainSelection = 0 correct it for HLGRATIO=16
1324   if(stackId == EMC && emcGainSelection==0) {
1325     rped  /= HLGRATIO;
1326     rampl /= HLGRATIO;
1327     cped  /= HLGRATIO;
1328     campl /= HLGRATIO;
1329   }
1330   if(campl<TWRAMPTHR) return NSAMPLES;
1331   E     = campl*twrECalib;
1332   //  store event data
1333   twrPed  -> Fill(cped);
1334   twrAmpl -> Fill(campl);
1335   twrTime -> Fill(ctime);
1336   twrE    -> Fill(E);
1337   //  trigger primitive (thresholding)
1338   // if(stackId==HOUTER) {
1339   //    display();
1340   // }
1341   //  if(!fitShape) return reject;
1342   return satFlag[gainToUse];
1343 }
1344 // **************************************************************************
1345
1346 //   clean tower for a new event
1347 void tower::clean(){
1348   for(int ig = 0; ig<gains; ig++)    {rawAmpl[ig] = 0.; rawTime[ig] = 0.; rawPed[ig] = 0.; satFlag[ig]=0;}
1349   reject   =  kFALSE;
1350   rped     = 0.;
1351   rampl    = 0.;
1352   rtime    = 0.;
1353   cped     = 0.;
1354   campl    = 0.;
1355   ctime    = 0.;
1356   E        = 0.;
1357   gainToUse  = 0;
1358   attached   =  kFALSE;
1359 }
1360
1361 // **************************************************************************
1362
1363 //   clean stack for a new event
1364 void tower::print(){
1365   cout<<"<tower::print>: St/twr "<<stackId<<"/"<<index<<"  X/Y "<<x<<"/"<<y<<" gainToUse "<<gainToUse<<"  Channels "<<adcChannel[0]<<"/"<<adcChannel[1]<<"  SatFlag "<<satFlag[0]<<"/"<<satFlag[1]<<"  Raw ampl "<<rawAmpl[0]<<" / "<<rawAmpl[1]<<"  Raw time "<<rawTime[0]<<" / "<<rawTime[1]<<"  Corrected ampl/time "<<campl<<"/"<<ctime<<"  Reject "<<reject<<endl;
1366 }
1367
1368
1369 // **************************************************************************
1370
1371 void tower::display(){
1372   hcalHelper * hHelper = hcalHelper::getInstance();
1373   TString twrN = "twr";    twrN += stackId; twrN += index;
1374   TString twrT = "Stack "; twrT += stackId; twrT += " Tower "; twrT += index; twrT += "Run "; twrT += hHelper->runnumber;
1375   TCanvas * twrDisplay = (TCanvas *)(gROOT->FindObject(twrN));
1376   if(!twrDisplay) twrDisplay  = new TCanvas(twrN, twrT, 0, 0, 400, 400);
1377   else twrDisplay->Clear();
1378   twrDisplay -> Divide(2,2);
1379   for(int ig=0; ig<gains;ig++){
1380     twrDisplay -> cd(ig*2+1);
1381     graph[ig]->SetMarkerStyle(20);
1382     graph[ig]->SetMarkerSize(0.4);
1383     graph[ig]->Draw("acp");
1384     twrDisplay -> cd(ig*2+2);
1385     twrAmpl  ->Draw();
1386   }
1387   twrDisplay->Update();
1388   print();
1389 }
1390