Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hcalHelper.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hcalHelper.C
1 
2 // *************************************************************************
3 hcalHelper::hcalHelper(){
4  if(ACTIVECHANNELS) {
5  cout<<"<hcalHelper::hcalHelper> Wrong instantiation sequence - hLabHelper already exists. EXITING"<<endl;
6  delete this;
7  }
8  ACTIVECHANNELS = 192; // just to insure everything works as it must
9  hLabHelper::getInstance();
10  eventseq = 0;
11  eventsread = 0;
12  displayMode = 3; // default - only final summary is displayed
13  // TODO: needs better understanding - how to handle HG Gains
14  // create calorimeter
15 
16  // we create all graphic objects in the directory of temporary file (rootTempFileName)
17 
18  t1044 = new hcal(); // it does stacks configuring at this time
19  // ONLY ACTIVE channels will be copyed from .prdf to adc
20  CHTOTAL = detchannels;
21  ACTIVECHANNELS = CHTOTAL;
22  channels = ACTIVECHANNELS;
23  samples = NSAMPLES;
24  activeCh = new Int_t [ACTIVECHANNELS];
25  adc = new Int_t * [ACTIVECHANNELS];
26  for (int ich=0; ich<ACTIVECHANNELS; ich++) {
27  adc[ich] = new Int_t [NSAMPLES];
28  }
29  // parse calibration information
30  t1044->setCalibration();
31  // make a look-up table for easy access to adc channels
32  Int_t ilive = 0;
33  for(int istk = 0; istk<CALSTACKS; istk++) {
34  if(t1044->alive[istk]) {
35  t1044->stacks[istk]->updateMap(ilive);
36  Int_t stackid = t1044->stacks[istk]->stackId;
37  Int_t chtocp = t1044->stacks[istk]->twrsInStack*t1044->stacks[istk]->gains;
38  for(int ich = 0; ich < chtocp; ich++) {
39  activeCh[ilive] = (stackid==EMC? emcCh[ich] :
40  (stackid==HINNER? hcalInnerCh[ich] :
41  (stackid==HOUTER? hcalOuterCh[ich] :
42  (stackid==HODO? hodoCh[ich] :
43  counters[ich]))));
44  ilive ++;
45  }
46  // t1044->stacks[istk]->print();
47  }
48  }
49  cout<<"<hcalHelper::hcalHelper>: "<<ACTIVECHANNELS<<" Channels Initialized"<<endl;
50  rdata = new TH2F("rdata", "Raw Peak Values for all active channels",ACTIVECHANNELS, 0, ACTIVECHANNELS, 1024, -2048., 2048.);
51 
52 }
53 
54 
55 // **************************************************************************
56 void hcalHelper::setRunKind(TString kind) {
57  hLabHelper * hlHelper = hLabHelper ::getInstance();
58  hcalHelper * hHelper = hcalHelper ::getInstance();
59  hlHelper -> setRunKind(kind);
60  if(kind=="cosmics") {
61  hHelper->t1044->alive[EMC] = kFALSE;
62  hHelper->t1044->alive[HODO] = kFALSE;
63  hHelper->t1044->alive[SCINT] = kFALSE;
64  delete hHelper->t1044->stacks[EMC];
65  delete hHelper->t1044->stacks[HODO];
66  delete hHelper->t1044->stacks[SCINT];
67  hHelper->t1044->stacks[EMC] = NULL;
68  hHelper->t1044->stacks[HODO] = NULL;
69  hHelper->t1044->stacks[SCINT] = NULL;
70  }
71 }
72 
73 // **************************************************************************
74 
75 int hcalHelper::evLoop(int run, int evToProcess, int fToProcess){
76  runnumber = run;
77  if(runnumber<=0) return 0;
78  hLabHelper * hlHelper = hLabHelper ::getInstance();
79  hlHelper->runnumber = runnumber;
80  hlHelper -> setPRDFRun(runnumber, kFALSE);
81  hlHelper -> makeCanvasDirectory();
82  hcalTree::getInstance();
83  // hlHelper-> fhcl->cd();
84 
85  int OK;
86 
87  Eventiterator *it = new fileEventiterator(hlHelper->prdfName, OK);
88  if (OK)
89  {
90  cout <<"<evLoop> Couldn't open input file " << hlHelper->prdfName << endl;
91  delete(it);
92  return 0;
93  }
94  Event *evt;
95  // Loop over events
96  while ( (evt = it->getNextEvent()) != 0 )
97  {
98  t1044 -> readyForEvent = kFALSE;
99  eventseq = evt->getEvtSequence();
100  eventReject = 0;
101  eventTrigger = 0;
102  if ( evt->getEvtType() != 1 ) {
103  cout << "eventseq " << eventseq << " event type = " << evt->getEvtType() << endl;
104  continue;
105  }
106 
107  // Get sPHENIX Packet
108  Packet_hbd_fpgashort *sp21101 =
109  dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101));
110  Packet_hbd_fpgashort *sp21102 =
111  dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21102));
112  if(sp21101&&sp21102){
113  sp21101->setNumSamples( NSAMPLES ); sp21102->setNumSamples( NSAMPLES );
114  for (int ich=0; ich<ACTIVECHANNELS; ich++) {
115  Int_t spCh = activeCh[ich]<144? activeCh[ich] :activeCh[ich]-144;
116  Packet_hbd_fpgashort* spacket = activeCh[ich]<144? sp21101 : sp21102;
117  for (int is=0; is<NSAMPLES; is++) {
118  adc[ich][is] = spacket->iValue(spCh,is);
119  // cout<<"<evLoop Event "<< hlHelper->eventseq<<" ich/is/adc "<<ich<<" - "<<is<<" - "<< hlHelper->adc[ich][is]<<" active "<<hlHelper->active[ich]<<endl;
120  }
121  }
122  delete sp21101; delete sp21102;
123  } else {
124  cout << "EventSeq " << eventseq << " type = " << evt->getEvtType() <<" Missing packets " <<(int)sp21101<<"/"<<(int)sp21102<<endl;
125  continue;
126  }
127  eventsread++;
128  if(eventsread%1000==0) cout <<"RUN "<<runnumber<<" EVENT "<<eventseq <<" - "<<eventsread<<endl;
129  // invert raw data in non-calorimeter stacks to make all signals negative
130  // subtract 2048 everywhere
131  for (int istk = 0; istk<CALSTACKS; istk++){
132  if(!t1044->alive[istk]) continue;
133  for (int itw=0; itw<t1044->stacks[istk]->twrsInStack; itw++){
134  for(int ig=0; ig<t1044->stacks[istk]->gains; ig++){
135  int iadc = t1044->stacks[istk]->towers[itw]->adcChannel[ig];
136  for(int is=0; is<NSAMPLES; is++){
137  if(t1044->stacks[istk]->stackKind==CALOR) {
138  adc[iadc][is] = adc[iadc][is]-2048.;
139  } else {
140  adc[iadc][is] = (activeCh[iadc]%2? -(adc[iadc][is]-2048.) : adc[iadc][is]-2048.);
141  }
142  // knowing where this value goes we may immediately set satFlags
143  t1044->stacks[istk]->towers[itw]->graph[ig]->SetPoint(is,(Double_t)is,adc[t1044->stacks[istk]->towers[itw]->adcChannel[ig]][is]);
144  if(adc[iadc][is]<undflow||adc[iadc][is]>ovrflow)
145  t1044->stacks[istk]->towers[itw]->satFlag[ig] += 1;
146  } // is
147  } // ig
148  } // itw
149  } // istk
150  // amount of data displayed in per event and per run display
151  // 0 - nothing displayed except Summary at the end of processing
152  // 1 - digitizations for everything
153  // 2 - lego plots of individual amplitudes in calorimeter and digitization for counter
154  Int_t rejectST = t1044->getStackTiming();
155  if(rejectST<3 || !triggerOn) { // HCal is not entirely empty or I want to see all events
156  Int_t rejectCR = collectRaw();
157  Int_t rejectCL = t1044->update(displayMode, kFALSE);
158  eventReject = rejectST*1000 + rejectCR*100 + rejectCL;
159  eventTrigger = t1044->collectTrPrimitives();
160  if(displayMode<3) cout<<"Event "<<eventseq<<" Trigger "<<eventTrigger<<" eventReject "<<eventReject<<" "<<rejectST<<"/"<<rejectCR<<"/"<<rejectCL<<endl;;
161  t1044->displayRaw(displayMode);
162  }
163 
164  delete evt;
165  if(displayMode<3){
166  // if(std::cin.get()==(Int_t)"q") break;
167  if(std::cin.get()==113) break;
168  // htree -> thcl -> Fill();
169  }
170  t1044->clean();
171  if ((evToProcess>0&&eventsread>=evToProcess)||eventsread>270000) break;
172  }
173  // display (print) un summary. Verbosity mode varies
174  t1044->displaySummary(displayMode);
175 
176 
177  cout<<"ALLDONE for RUN "<<runnumber<<" Events "<<eventsread<<endl;
178  delete it;
179  hlHelper -> thcl -> Write();
180  hlHelper -> fhcl -> Write();
181  // TODO - if multiple files are processed closing the root file should be done somewhere else
182  // TODO - at the moment
183  // hlHelper -> renameAndCloseRF();
184  // closeLoop();
185  return eventseq;
186 }
187 
188 
189 // **************************************************************************
190 // amount of per event and per run display
191 // 0 - nothing displayed except Summary at the end of processing
192 // 1 - digitizations for everything
193 // 2 - lego plots of individual amplitudes in calorimeter and digitization for counters
194 
195 void hcalHelper::setDisplayMode(Int_t mode){
196  displayMode = mode;
197  cout<<"<hcalHelper::setDisplayMode> Mode is set to "<<displayMode<<endl;
198 }
199 
200 // **************************************************************************
201 
202 void hcalHelper::setCTriggerOn(Bool_t ON){
203  triggerOn = ON;
204  cout<<"<hcalHelper::setCTriggerOn> triggerOn is set to "<<triggerOn<<endl;
205 }
206 
207 // **************************************************************************
208 
209 void hcalHelper::setCTriggerGRange(Int_t gain){
210  triggerGain = gain;
211  TRGAINRANGE = gain;
212  cout<<"<hcalHelper::setDCTriggerGRange> Gaine is set to "<<triggerGain<<endl;
213 }
214 
215 // **************************************************************************
216 void hcalHelper::updateCalibration(){
217  // hLabHelper * hlHelper = hLabHelper ::getInstance();
218  ;
219 }
220 
221 // **************************************************************************
222 
223 void hcalHelper::hcalPattern(Int_t nx, Int_t ny, Int_t run, Int_t mod){;}
224 
225 // **************************************************************************
226 
227 // 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
228 Int_t hcalHelper::collectRaw(){
229  // currently no rejection at this stage
230  rejectRaw = 0;
231  // stacks first
232  for(int istk = 0; istk<CALSTACKS; istk++){
233  if(!(t1044->stacks[istk])) continue;
234  for(int ig = 0; ig<t1044->stacks[istk]->gains; ig++){
235  // Double_t fPeak = t1044->stacks[istk]->fitPeak[ig];
236  Int_t fTime = (int)(t1044->stacks[istk]->fitTime[ig]+0.5);
237  for (int itw = 0; itw<t1044->stacks[istk]->twrsInStack; itw++){
238  Int_t ich = t1044->stacks[istk]->towers[itw]->adcChannel[ig];
239  // extract raw peak position and peak value for this channel
240  Int_t iss = max((fTime-2), 1);// Int_t ise = min((fTime+3), NSAMPLES);
241  Double_t psum(0.); float sUsed(0.);
242  for (int is = 0; is<iss; is++) { if(adc[ich][is]>=undflow&&adc[ich][is]<=ovrflow) {psum += adc[ich][is]; sUsed += 1.;}}
243 
244  if(sUsed) psum /= sUsed; else psum=0.;
245  Double_t rVal = adc[ich][fTime] - psum;
246  t1044->stacks[istk]->towers[itw]->rawPed[ig] = psum;
247  t1044->stacks[istk]->towers[itw]->rawAmpl[ig] = rVal;
248  t1044->stacks[istk]->towers[itw]->rawTime[ig] = fTime; // rTime;
249  rdata -> Fill(t1044->stacks[istk]->towers[itw]->adcChannel[ig], rVal);
250  }
251 
252  }
253  // all kinds of scintillators etc
254  }
255  return rejectRaw;
256 }
257 // **************************************************************************
258 void hcalHelper::hcalImpact(){;}
259 // **************************************************************************
260 void hcalHelper::fitHCalSignal(){;}
261 // **************************************************************************
262 void hcalHelper::hcalDisplay(){;}
263 // **************************************************************************
264 // FIXES to data
265 void hcalHelper::dofixes(){
266  // if(runnumber>=1125 && runnumber<1152){
267  // // invert pulse in channel 134 (number 9 in todays parlance
268  // for (int is = 0; is < NSAMPLES; is++) {adc[8][is] -= 2*PEDESTAL; adc[8][is] = abs(adc[8][is]);}
269  // }
270  ;
271 }
272 
273 // **************************************************************************
274 
275 Int_t hcalHelper::reject(){
276  return false;
277 }
278 
279 // **************************************************************************
280 hcal::hcal(){
281 
282  maxStacks = CALSTACKS; // CALSTACKS
283  stacks = new stack * [maxStacks];
284  alive = new Bool_t [maxStacks];
285  kinds = new Int_t [maxStacks];
286  kinds[EMC] = CALOR; alive[EMC] = kTRUE; // EMC
287  kinds[HINNER] = CALOR; alive[HINNER] = kTRUE; // HInner
288  kinds[HOUTER] = CALOR; alive[HOUTER] = kTRUE; // Outer HCal
289  kinds[HODO] = COUNTER; alive[HODO] = kTRUE; // HODOSCOPE
290  kinds[SCINT] = COUNTER; alive[SCINT] = kTRUE; // Beam counters
291  kinds[CHER] = COUNTER; alive[CHER] = kFALSE; // Cherenkov's
292  activeStacks = 0;
293  activeCalStacks = 0;
294 
295  for (int istk = 0; istk < maxStacks; istk++) stacks[istk] = NULL;
296 
297  // this is just a proyotype, it has well defined configuration
298  detchannels = 0;
299  hgDetChannels = 0;
300  lgDetChannels = 0;
301  if(alive[EMC]) {
302  stacks[EMC] = new stack(EMC, EMCTOWERS, EMCCOLUMNS, EMCROWS);
303  // cout<<"MMM "<<stacks[EMC]->mapUpdated<<endl;
304  stacks[EMC]->stackKind = kinds[EMC];
305  activeStacks ++;
306  activeCalStacks ++;
307  detchannels += EMCTOWERS*EMCGAINS; // only one gain (treat as LG)
308  lgDetChannels += EMCTOWERS*EMCGAINS;
309  // cout<<"MMM "<<stacks[EMC]->mapUpdated<<endl;
310  }
311  if(alive[HINNER]) {
312  stacks[HINNER] = new stack(HINNER, HCALTOWERS, HCALCOLUMNS, HCALROWS);
313  stacks[HINNER]->stackKind = kinds[HINNER];
314  activeStacks ++;
315  activeCalStacks ++;
316  detchannels += HCALTOWERS*HCALGAINS; // dual gain readout in HCal
317  hgDetChannels += HCALTOWERS;
318  lgDetChannels += HCALTOWERS;
319  }
320  if(alive[HOUTER]) {
321  stacks[HOUTER] = new stack(HOUTER, HCALTOWERS, HCALCOLUMNS, HCALROWS);
322  stacks[HOUTER]->stackKind = kinds[HOUTER];
323  activeStacks ++;
324  activeCalStacks ++;
325  detchannels += HCALTOWERS*HCALGAINS;
326  hgDetChannels += HCALTOWERS;
327  lgDetChannels += HCALTOWERS;
328  }
329  if(alive[HODO]) {
330  stacks[HODO] = new stack(HODO, T1044HODO, T1044HODO/2, 2);
331  stacks[HODO]->stackKind = kinds[HODO];
332  activeStacks ++;
333  detchannels += T1044HODO;
334  lgDetChannels += T1044HODO;
335  }
336  if(alive[SCINT]) {
337  stacks[SCINT] = new stack(SCINT, T1044TRIGGERCH, T1044TRIGGERCH/2, 2);
338  stacks[SCINT]->stackKind = kinds[SCINT];
339  activeStacks ++;
340  detchannels += T1044TRIGGERCH;
341  lgDetChannels += T1044TRIGGERCH;
342  }
343  for(int istk = 0; istk<CALSTACKS; istk++) {
344  if(alive[istk]&&stacks[istk]->mapUpdated) cout<<"Stack "<<stacks[istk]->stackId<<" mapUpdated already set "<<stacks[istk]->mapUpdated<<endl;
345  }
346  readyForEvent = kTRUE;
347  clean();
348 
349  // Summary histos for the whole calorimeter
350  TString totaN = "tota";
351  TString totaT = "Total Amplitude in Calorimeter System";
352  totalamp = new TH1F(totaN, totaT, 10000, 0., 10000.);
353  TString toteN = "tote";
354  TString toteT = "Total Energy in Calorimeter System";
355  totale = new TH1F(toteN, toteT, 1000, 0., 60.);
356  TString tlcgN = "tlcg";
357  TString tlcgT = "Shower Center of gravity (units of Stacks)";
358  totallcg = new TH1F(tlcgN, tlcgT, 200, 0., 2.);
359  TString sclcgN = "sclcg";
360  TString sclcgT = "Shower Center in scintillators (units of Stacks)";
361  scintlcg = new TH1F(sclcgN, sclcgT, 200, 0., 2.);
362 }
363 
364 // **************************************************************************
365 
366 void hcal::setCalibration(){
367  // set energy scales
368  for(int istk = 0; istk<CALSTACKS; istk++){
369  if(!alive[istk]) continue;
370  stacks[istk]->stackECalib = stECalib[stacks[istk]->stackId];
371  // TODO: set correction coefficients in towers (currently unity)
372  if(stacks[istk]->stackKind == CALOR) {
373  for (int itw=0; itw< stacks[istk]->twrsInStack; itw++) {
374  stacks[istk]->towers[itw]->twrECalib = stECalib[stacks[istk]->stackId];
375  // must be scaled by the Muon peak positions
376  stacks[istk]->towers[itw]->twrEScale = 1.;
377  }
378  }
379  }
380  //
381 }
382 
383 // **************************************************************************
384 hcal::~hcal(){
385  if(totalamp) delete totalamp;
386  if(totale) delete totale;
387  if(totallcg) delete totallcg;
388  if(scintlcg) delete scintlcg;
389 
390 }
391 
392 // **************************************************************************
393 
394 Int_t hcal::getStackTiming(){
395  for(int istk = 0; istk<CALSTACKS; istk++){
396  if(!stacks[istk]) continue;
397  Int_t stackReject = stacks[istk]->getStackTime();
398  if(kinds[istk] == CALOR && stackReject>10) rejectEvent ++ ;
399  }
400  // raise reject flag
401  return rejectEvent;
402 }
403 
404 // **************************************************************************
405 
406 Int_t hcal::collectTrPrimitives(){
407  int accept = 0;
408 
409  // CALOR primitives
410  for(int istk = 0; istk<CALSTACKS; istk++){
411  if(!stacks[istk]) continue;
412  if(!stacks[istk]->stackKind == CALOR) continue;
413  // kill LED events at this time
414  Int_t triggerFlag = stacks[istk]->getTrigger();
415  accept += triggerFlag*pow(100,stacks[istk]->stackId);
416  }
417  return accept;
418 }
419 
420 
421 // **************************************************************************
422 // done after collectRaw already collected and updated initial values for peaks and times for all objects (towers, counters)
423 // It also updates graphical objects in towers so if you need individual towers visible - take care
424 Int_t hcal::update(Int_t displayMode, Bool_t fitShapes){
425  // updates & compute unique amplitudes (in objects with multiple ranges), times and energies
426  for(int istk = 0; istk<CALSTACKS; istk++){
427  if(!stacks[istk]) continue;
428  // update adds 1000 to whatever is already stored in reject (saturation in towers, below zero threshold). Returns reject
429  Int_t sr = stacks[istk]->update(fitShapes); // stack reject
430  if(stacks[istk]->stackKind==CALOR&&sr) rejectEvent += pow(10,stacks[istk]->stackId);
431  }
432  // TODO: collect global event information for presentation
433 
434  for(int istk = 0; istk<CALSTACKS; istk++){
435  if(!alive[istk]||stacks[istk]->stackKind!=CALOR) continue;
436  // stacks[istk]->print();
437  amplTotal += (stacks[istk]->totCorAmpl)*stAScale[istk];
438  eTotal += stacks[istk]->E;
439  calLCG += (float)(stacks[istk]->stackId)*stacks[istk]->totCorAmpl*stAScale[istk];
440  }
441  calLCG /= amplTotal;
442  totalamp -> Fill(amplTotal);
443  totale -> Fill(eTotal);
444  totallcg -> Fill(calLCG);
445  // now scintillators
446  if(!alive[SCINT]) return 0;
447  Double_t tcampl(0.);
448  for (int sc = 0; sc<3; sc++){
449  scintLCG += (sc+1)*stacks[SCINT]->towers[sc+3]->campl;
450  tcampl += stacks[SCINT]->towers[sc+4]->campl;
451  }
452  scintLCG /= tcampl;
453  scintlcg -> Fill(scintLCG);
454  //
455  return 0;
456 }
457 // **************************************************************************
458 void hcal::displayRaw(Int_t mode){
459  if(!mode || mode>2) return;
460  for(int istk = 0; istk<CALSTACKS; istk++){
461  if(!stacks[istk]) continue;
462  // kill LED
463  // if(stacks[istk]->triggerFlag && stacks[istk]->triggerHits<12) {
464  if(mode==1||(mode==2&&stacks[istk]->stackKind==CALOR)) {
465  stacks[istk]->displayEvent(mode);
466  // if(mode==1)
467  stacks[istk]->displayADCSum();
468  }
469  // }
470  }
471 }
472 
473 // **************************************************************************
474 void hcal::displaySummary(Int_t mode){
475  hcalHelper * hHelper = hcalHelper ::getInstance();
476  hLabHelper * hlHelper = hLabHelper ::getInstance();
477  for(int istk = 0; istk<CALSTACKS; istk++){
478  if(!stacks[istk]) continue;
479  stacks[istk]->displayTowerSummary(mode);
480  stacks[istk]->displayStackSummary(mode);
481  }
482  TCanvas * calSum;
483  TString csN = "csSum";
484  TString csT = "Calorimeter Summary for Run "; csT += hHelper->runnumber;
485  calSum = (TCanvas *)(gROOT->FindObject(csN));
486  if(!calSum) calSum = new TCanvas(csN, csT, 50, 50, 500, 250);
487  else calSum->Clear();
488  calSum -> SetEditable(kTRUE);
489  calSum -> Divide(4,1);
490  calSum -> cd(1);
491  totalamp -> Draw();
492  calSum -> cd(2);
493  totale -> Draw();
494  calSum -> cd(3);
495  totallcg -> Draw();
496  calSum -> cd(4);
497  scintlcg -> Draw();
498  calSum -> Update();
499 
500  TCanvas * satSum;
501  TString satsN = "satsum";
502  TString satsT = "Saturation in Calorimeters Run "; satsT += hHelper->runnumber;
503  satSum = (TCanvas *)(gROOT->FindObject(satsN));
504  if(!satSum) satSum = new TCanvas(satsN, satsT, 50, 50, 600, 400);
505  else satSum->Clear();
506  satSum->Divide(2,3);
507  for(int istk=0; istk<CALSTACKS; istk++) {
508  if(!alive[istk]||stacks[istk]->stackKind != CALOR) continue;
509  if(stacks[istk]->stackId==EMC) {
510  satSum -> cd(1);
511  satSum -> GetPad(1)->SetLogy();
512  stacks[istk]->satProb[0]->Draw();
513  satSum -> cd(2);
514  totale -> Draw();
515  }
516  if(stacks[istk]->stackId==HINNER) {
517  satSum -> cd(3);
518  satSum -> GetPad(3)->SetLogy();
519  stacks[istk]->satProb[0]->Draw();
520  satSum -> cd(4);
521  satSum -> GetPad(4)->SetLogy();
522  stacks[istk]->satProb[1]->Draw();
523  }
524  if(stacks[istk]->stackId==HOUTER) {
525  satSum -> cd(5);
526  satSum -> GetPad(5)->SetLogy();
527  stacks[istk]->satProb[0]->Draw();
528  satSum -> cd(6);
529  satSum -> GetPad(6)->SetLogy();
530  stacks[istk]->satProb[1]->Draw();
531  }
532  }
533  satSum->Update();
534  if(!hlHelper->cDirectory.IsNull()){
535  TString cName = hlHelper->cDirectory;
536  cName += "calSum_"; cName += hHelper->runnumber; cName += ".gif";
537  calSum->SaveAs(cName);
538  cName = hlHelper->cDirectory; cName += "satSum_"; cName += hHelper->runnumber; cName += ".gif";
539  satSum->SaveAs(cName);
540  }
541 }
542 // **************************************************************************
543 
544 void hcal::clean(){
545  for(int istk = 0; istk<CALSTACKS; istk++){
546  if(!stacks[istk]) continue;
547  stacks[istk]->clean();
548  }
549  rejectEvent = 0;
550  readyForEvent = true;
551  amplTotal = 0.;
552  eTotal = 0.;
553  calLCG = 0.;
554  scintLCG = 0.;
555 }
556 
557 // **************************************************************************
558 
559 stack::stack(const Int_t sId, const Int_t twrs, const Int_t xch, const Int_t ych){
560  reject = 0;
561  mapUpdated = kFALSE;
562  eventsSeen = 0;
563  // cout<<"SSS STACK "<<mapUpdated<<endl;
564  triggerFlag = 0;
565  triggerHits = 0;
566  triggerSum = 0;
567  stackECalib = 0.;
568  trAmp = NULL;
569  stackId = sId;
570  twrsInStack = twrs;
571  stackLoc = -1; // location of this stack in ADC still not defined
572  nTwrsX = xch;
573  nTwrsY = ych;
574  gains = (stackId==HINNER || stackId==HOUTER)? 2 : 1;
575  towers = new tower * [twrsInStack];
576  for (int itw = 0; itw<twrsInStack; itw++){
577  int xSt = itw%nTwrsX;
578  int ySt = itw/nTwrsX;
579  towers[itw] = new tower(stackId, itw, xSt, ySt);
580  towers[itw]->adcChannel[0] = itw;
581  towers[itw]->adcChannel[1] = (stackId==HINNER || stackId==HOUTER)? itw+HCALTOWERS : -1;
582  }
583  // cout<<"SSS0 "<<mapUpdated<<endl;
584  for (int ig = 0; ig<gains; ig++){
585  adcsum[ig] = new Double_t [NSAMPLES];
586  fitPar[ig] = new Double_t [NPARAMETERS];
587  TString gT = "Stack "; gT += stackId; gT += (ig==0? " HIGH " : " LOW "); gT += "gain";
588  graph[ig] = new TGraph(NSAMPLES);
589  graph[ig] -> SetTitle(gT);
590  TString sum = "adcsum_"; sum += (ig==HIGH)? "HG" : "LG";
591  TString sig = "signal_"; sig += (ig==HIGH)? "HG" : "LG";
592  shape[ig] = new TF1(sum, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
593  signal[ig] = new TF1(sig, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
594  // Double_t times[ACTIVECHANNELS];
595  // Double_t chi2 [ACTIVECHANNELS];
596  }
597 
598  // cout<<"SSS1 "<<mapUpdated<<endl;
599 
600  E = 0.;
601  // summaries from FIT
602  // Float_t maxAmpSum = (stackId==EMC&&emcGainSelection==0)? 2000. : 1000.;
603  // Float_t maxPedSum = (stackId==EMC&&emcGainSelection==0)? 250. : 25.;
604  Float_t maxPedSum = 400.;
605  Float_t maxAmpSum = 16000.;
606  TString stfpN = "st_fped"; stfpN += stackId;
607  TString stfpT = "Stack "; stfpT += stackId; stfpT += " Pedestals - Fit [LG]";
608  stFPed = new TH1F(stfpN, stfpT, 100, -maxPedSum, maxPedSum);
609  TString stfaN = "st_famp"; stfaN += stackId;
610  TString stfaT = "Stack "; stfaT += stackId; stfaT += " Amplitudes - Fit [LG]";
611  stFAmpl = new TH1F(stfaN, stfaT, 550, -0.1*maxAmpSum, maxAmpSum);
612  TString stftN = "st_ftime"; stftN += stackId;
613  TString stftT = "Stack "; stftT += stackId; stftT += " Time - Fit [ADC ticks]";
614  stFTime = new TH1F(stftN, stftT, NSAMPLES, 0., NSAMPLES);
615  TString stc2N = "st_c2"; stc2N += stackId;
616  TString stc2T = "Stack "; stc2T += stackId; stc2T += " Fit Chi2 ";
617  stChi2 = new TH1F(stc2N, stc2T, 100, 0., 200.);
618  TString stgrN = "st_gr"; stgrN += stackId;
619  TString stgrT = "Stack "; stgrT += stackId; stgrT += " Ratio HG/LG ";
620  stFGR = new TH1F(stgrN, stgrT, 100, -20., 20.);
621 
622  // cout<<"SSS2 "<<mapUpdated<<endl;
623  // summaries from contributing towers
624  TString stspN = "st_sped"; stspN += stackId;
625  TString stspT = "Stack "; stspT += stackId; stspT += " Stack Pedestals Sum [LG]";
626  stSPed = new TH1F(stspN, stspT, 100, -maxPedSum, maxPedSum);
627  TString stsaN = "st_samp"; stsaN += stackId;
628  TString stsaT = "Stack "; stsaT += stackId; stsaT += " Stack Amplitude Sum [LG]";
629  stSAmpl = new TH1F(stsaN, stsaT, 550, -0.1*maxAmpSum, maxAmpSum);
630  TString stseN = "st_sen"; stseN += stackId;
631  TString stseT = "Stack "; stseT += stackId; stseT += " Stack Energy Sum [MeV]";
632  stSE = new TH1F(stseN, stseT, 1100, -5., 50.);
633  TString ststN = "st_stime"; ststN += stackId;
634  TString ststT = "Stack "; ststT += stackId; ststT += " <Stack Time> [clock ticks]";
635  stAvT = new TH1F(ststN, ststT, NSAMPLES, 0., NSAMPLES);
636  TString ststsN = "st_stsig"; ststsN += stackId;
637  TString ststsT = "Stack "; ststsT += stackId; ststsT += " Time - RMS [clock ticks]";
638  stTRMS = new TH1F(ststsN, ststsT, NSAMPLES, 0., NSAMPLES);
639  // Hit impacts in the stack
640  TString hitCGN = "hitSt"; hitCGN += stackId;
641  TString hitCGT = "Stack "; hitCGT += stackId; hitCGT += " Impact X vs Y [towers]";
642  hitCG = new TH2F(hitCGN, hitCGT, nTwrsX, 0, nTwrsX, nTwrsY, 0, nTwrsY);
643 
644  // SATURATION
645  for(int ig=0; ig<gains; ig++) {
646  TString satN = "stsat"; satN += stackId; satN += ig;
647  TString satT = "Stack "; satT += stackId; satT += ig==0? " HG " : " LG "; satT += "saturation probability vs # of sat. samples";
648  satProb[ig] = new TH1F(satN, satT, 25, 0., 25.);
649  }
650 
651 
652  if(stackKind != CALOR/*||hlHelper->runKind!="cosmics"*/) return;
653  // cosmic related summaries
654  TString trAmpN = "trAmp"; trAmpN += stackId;
655  TString trAmpT = "Stack "; trAmpT += stackId; trAmpT += " Track Amplitude";
656  trAmp = new TH1F(trAmpN, trAmpT, 300, 0., 300.);
657  TString trCh2N = "trCh2"; trCh2N += stackId;
658  TString trCh2T = "Stack "; trCh2T += stackId; trCh2T += " Track Chi2";
659  trChi2 = new TH1F(trCh2N, trCh2T, 100, 0., 10.);
660  TString trCrN = "trCr"; trCrN += stackId;
661  TString trCrT = "Stack "; trCrT += stackId; trAmpT += " Track Crossing";
662  trCr = new TH1F(trCrN, trCrT, 40, 0., 4.);
663  TString trSlN = "trSl"; trSlN += stackId;
664  TString trSlT = "Stack "; trSlT += stackId; trAmpT += " Track Slope";
665  trSl = new TH1F(trSlN, trSlT, 50, -1., 1.);
666  TString trACh2N = "trACh2"; trACh2N += stackId;
667  TString trACh2T = "Stack "; trACh2T += stackId; trACh2T += " Track Amplitude vs Chi2";
668  trAmpCh2 = new TH2F(trACh2N, trACh2T, 300, 0, 300, 100, 0., 2.);
669  TString trAstHN = "trAstH"; trAstHN += stackId;
670  TString trAstHT = "Stack "; trAstHT += stackId; trAstHT += " Track Amplitude vs # of Hits in Stack";
671  trAstH = new TH2F(trAstHN, trAstHT,300, 0, 300, 100, 0., 16.);
672  TString trASlN = "trASl"; trASlN += stackId;
673  TString trASlT = "Stack "; trASlT += stackId; trASlT += " Track Amplitude vs # Track Slope";
674  trASl = new TH2F(trASlN, trASlT, 300, 0, 300, 100, 0., 1.);
675  // Next scatterplot is the basis of calibration
676  // Its data are HG data divided by H/L gain ratio and so to make this
677  // picture representative the number of bins is chose equal to
678  // amplimit*hgratio[stackId]
679  // cout<<"SSS4 "<<mapUpdated<<endl;
680  Float_t amplimit = 100; // LG ADC counts
681  Int_t nx = amplimit*hlgratios[stackId];
682  TString trATwrN = "trATwrN"; trATwrN += stackId;
683  TString trATwrT = "Stack "; trATwrT += stackId; trATwrT += " Tagged Twr Amplitude vs Twr # ";
684  trATwr = new TH2F(trATwrN, trATwrT, nx, 0, 100, twrsInStack, 0., twrsInStack);
685  // cout<<"SSSExit "<<mapUpdated<<endl;
686 
687 }
688 
689 // **************************************************************************
690 stack::~stack(){
691  if(stackId < 0 || !towers) return;
692  for (int itw = 0; itw < twrsInStack; itw++){
693  if(towers[itw]) delete towers[itw];
694  }
695  delete stFGR;
696  delete stFPed;
697  delete stFAmpl;
698  delete stFTime;
699  delete stChi2;
700  delete stSPed;
701  delete stSAmpl;
702  delete stSE;
703  delete stAvT;
704  delete stTRMS;
705  delete hitCG;
706  if(trAmp) {
707  delete trAmp;
708  delete trChi2;
709  delete trCr;
710  delete trSl;
711  delete trAmpCh2;
712  delete trAstH;
713  delete trASl;
714  delete trATwr;
715  }
716 
717 }
718 
719 // **************************************************************************
720 
721 void stack::displayEvent(Int_t mode){
722  if(!mode || mode>2) return;
723  hcalHelper * hHelper = hcalHelper::getInstance();
724  TCanvas * pattern[2];
725  for(int ig = 0; ig < gains; ig++){
726  TString stckN = "StPat"; stckN += stackId; stckN += "_G_"; stckN += ig;
727  TString stckT = "Stack "; stckT += stackId; stckT += (ig==HIGH)? " HGain" : " LGain"; stckT += " Event Pattern Display Run "; stckT += hHelper->runnumber;
728  pattern[ig] = (TCanvas *)(gROOT->FindObject(stckN));
729  if(!pattern[ig]) pattern[ig] = new TCanvas(stckN, stckT, 50*(stackId+ig), 50*(stackId+ig), 100*nTwrsX, 100*nTwrsY);
730  else pattern[ig]->Clear();
731  pattern[ig] -> SetEditable(kTRUE);
732  // display digitizations
733  if(mode == 1 || stackKind != CALOR) {
734  pattern[ig] -> Divide(nTwrsX,nTwrsY);
735  for (int itw = 0; itw<twrsInStack; itw++){
736  Int_t ipad = (nTwrsY-(towers[itw]->y)-1)*nTwrsX + (towers[itw]->x)+1;
737  pattern[ig] -> cd(ipad);
738  // cout<<"<stack::displayPattern "<<nTwrsX<<" "<<nTwrsY<<" "<<x<<" "<<y<<" "<<ipad<<endl;
739  towers[itw]->graph[ig]->SetMarkerStyle(20);
740  towers[itw]->graph[ig]->SetMarkerSize(0.4);
741  towers[itw]->graph[ig]->Draw("acp");
742  // towers[itw]->shape[ig]->Draw("same");
743  }
744  pattern[ig]->Update();
745  } else {
746  TH2 * stackLego;
747  // display stack lego plot - only for calorimeter stacks
748  TString slN = "st"; slN += stackId; slN += "_G_"; slN += ig;
749  TString slT = "Mixed Amplitudes [LG] in Stack "; slT += stackId;
750  if(!(stackLego = (TH2F*) (gROOT->FindObject(slN)))) stackLego = new TH2F(slN, slT, nTwrsX, 0, nTwrsX, nTwrsY, 0, nTwrsY);
751  for (int itw = 0; itw<twrsInStack; itw++){
752  // towers[itw]->print();
753  // cout<<"stack "<<stackId<<" itw/x/y/E "<<itw<<" / "<<towers[itw]->x<<" / "<< towers[itw]->y<<" / "<< towers[itw]->campl<<endl;
754  stackLego->SetBinContent(1+towers[itw]->x, 1+towers[itw]->y, towers[itw]->campl);
755  }
756  stackLego -> Draw("lego");
757  pattern[ig]->Update();
758  if(ig==0) break;
759  }
760  }
761 }
762 
763 // **************************************************************************
764 
765 void stack::displayTowerSummary(Int_t mode){
766  if(mode<=2) return;
767  hLabHelper * hlHelper = hLabHelper ::getInstance();
768  hcalHelper * hHelper = hcalHelper ::getInstance();
769  TCanvas * corTwrsA;
770  TString cTAN = "rtA_st"; cTAN += stackId;
771  TString cTAT = "Stack "; cTAT += stackId; cTAT += " Hit Amplitudes in Towers [LG units] Run "; cTAT += hHelper->runnumber;
772  corTwrsA = (TCanvas *)(gROOT->FindObject(cTAN));
773  if(!corTwrsA) corTwrsA = new TCanvas(cTAN, cTAT, 50*stackId, 50*stackId, 100*nTwrsX, 100*nTwrsY);
774  else corTwrsA->Clear();
775  corTwrsA -> SetEditable(kTRUE);
776  corTwrsA -> Divide(nTwrsX,nTwrsY);
777  for (int itw = 0; itw<twrsInStack; itw++){
778  Int_t ipad = (nTwrsY-(towers[itw]->y)-1)*nTwrsX + (towers[itw]->x)+1;
779  corTwrsA -> cd(ipad);
780  corTwrsA -> GetPad(ipad)->SetLogy();
781  towers[itw]->twrAmpl->SetMarkerStyle(20);
782  towers[itw]->twrAmpl->SetMarkerSize(0.4);
783  towers[itw]->twrAmpl->Draw();
784  // rawTwrsT -> cd(ipad);
785  // // rawTwrsA -> GetPad(ipad)->SetLogy();
786  // towers[itw]->twrTime->SetMarkerStyle(20);
787  // towers[itw]->twrTime->SetMarkerSize(0.4);
788  // towers[itw]->twrTime->Draw();
789  }
790  corTwrsA->Update();
791  if(!hlHelper->cDirectory.IsNull()){
792  TString cName = hlHelper->cDirectory;
793  cName += "corTwrsA_"; cName += hHelper->runnumber; cName += ".gif";
794  corTwrsA->SaveAs(cName);
795  }
796 
797  // rawTwrsT->Update();
798 }
799 // **************************************************************************
800 
801 void stack::displayStackSummary(Int_t mode){
802  if(mode != 3 || stackKind != CALOR) return;
803  hcalHelper * hHelper = hcalHelper ::getInstance();
804  hLabHelper * hlHelper = hLabHelper ::getInstance();
805  TCanvas * stackSum(NULL), * muTracks(NULL);
806  TString stN = "stSum"; stN += stackId;
807  TString stT = "Stack "; stT += stackId; stT += " Summary for Run "; stT += hHelper->runnumber;
808  stackSum = (TCanvas *)(gROOT->FindObject(stN));
809  if(!stackSum) stackSum = new TCanvas(stN, stT, 50*stackId, 50*stackId, 400, 400);
810  else stackSum->Clear();
811  stackSum -> SetEditable(kTRUE);
812  stackSum -> Divide(4,2);
813  // TH1 * stFPed, * stFAmpl, * stFTime, * stChi2;
814  // TH1 * stSPed, * stSAmpl, * stSE, * stAvT, * stTRMS;
815  stackSum ->cd(5);
816  stSPed ->Draw();
817  stackSum ->cd(6);
818  stackSum -> GetPad(6)->SetLogy();
819  stSAmpl ->Draw();
820  stackSum ->cd(7);
821  stAvT ->Draw();
822  stackSum ->cd(8);
823  stackSum -> GetPad(8)->SetLogy();
824  stSE ->Draw();
825 
826  stackSum ->cd(1);
827  stFPed ->Draw();
828  stackSum ->cd(2);
829  stackSum -> GetPad(2)->SetLogy();
830  stFAmpl ->Draw();
831  stackSum ->cd(3);
832  stFGR ->Draw();
833  // stFTime ->Draw();
834  stackSum ->cd(4);
835  stChi2 ->Draw();
836  stackSum ->Update();
837  // saturation probability
838  if(eventsSeen) {
839  for (int ig=0; ig<gains; ig++) {
840  satProb[ig]->Scale(1./(float)eventsSeen);
841  satProb[ig]->Scale(1./(float)twrsInStack);
842  }
843  }
844 
845  if(!hlHelper->cDirectory.IsNull()){
846  TString cName = hlHelper->cDirectory;
847  cName += "stackSum_"; cName += stackId; cName += "_"; cName += hHelper->runnumber; cName += ".gif";
848  stackSum->SaveAs(cName);
849  }
850  if(hlHelper->runKind == "cosmics"){
851  TString muN = "muTr"; muN += stackId;
852  TString muT = "Stack "; muT += stackId; muT += " Cosmic Muons in Run "; muT += hHelper->runnumber;
853  muTracks = (TCanvas *)(gROOT->FindObject(muN));
854  if(!muTracks) muTracks = new TCanvas(muN, muT, 50*stackId, 50*stackId, 400, 400);
855  else muTracks->Clear();
856  muTracks -> SetEditable(kTRUE);
857  muTracks -> Divide(3,2);
858  muTracks -> cd(1);
859  trAmp -> Draw();
860  muTracks -> cd(2);
861  trCr -> Draw();
862  muTracks -> cd(3);
863  trSl -> Draw();
864  muTracks -> cd(4);
865  trASl -> Draw();
866  // trChi2 -> Draw();
867  muTracks -> cd(5);
868  trAmpCh2 -> Draw();
869  muTracks -> cd(6);
870  trAstH -> Draw();
871  muTracks -> Update();
872 
873  if(!hlHelper->cDirectory.IsNull()){
874  TString cName = hlHelper->cDirectory;
875  cName += "muTracks_"; cName += stackId; cName += "_"; cName += hHelper->runnumber; cName += ".gif";
876  muTracks->SaveAs(cName);
877  }
878  }
879 }
880 
881 // **************************************************************************
882 
883 void stack::displayADCSum(){
884  hcalHelper * hHelper = hcalHelper ::getInstance();
885  hLabHelper * hlHelper = hLabHelper ::getInstance();
886  TString stckN = "StEvS"; stckN += stackId;
887  TString stckT = "Stack "; stckT += stackId; stckT += " Stack Event Sum Display Run "; stckT += hHelper->runnumber;
888  TCanvas * stckEvSumDisplay = (TCanvas *)(gROOT->FindObject(stckN));
889  if(!stckEvSumDisplay) stckEvSumDisplay = new TCanvas(stckN, stckT, 25*stackId, 25*stackId, 600, 400);
890  else stckEvSumDisplay->Clear();
891  stckEvSumDisplay -> Divide(1,gains);
892  for(int ig = 0; ig < gains; ig++){
893  stckEvSumDisplay->cd(ig+1);
894  // for (int is = 0; is < NSAMPLES; is++) evtGraph[ig]->SetPoint(is,(Double_t)is,evtadcsum[ig][is]);
895  // cout<<"<tileDisplay> call to fitTileSignal"<<endl;
896  // eventSum->Divide(1,2);
897  // eventSum->cd(1);
898  // eventSum->SetEditable(kTRUE);
899  graph[ig]->SetMarkerStyle(20);
900  graph[ig]->SetMarkerSize(0.4);
901  graph[ig]->Draw("acp");
902  // eventSum->cd(2);
903  shape[ig]->Draw("same");
904  // evtShape[ig]->Draw();
905  // evtGraph[ig]->Draw("same");
906  // stckDisplay->Update();
907  }
908  stckEvSumDisplay->Update();
909  if(!hlHelper->cDirectory.IsNull()){
910  TString cName = hlHelper->cDirectory;
911  cName += "stackEvSum_"; cName += stackId; cName += "_"; cName += hHelper->runnumber; cName += ".gif";
912  stckEvSumDisplay->SaveAs(cName);
913  }
914 }
915 // **************************************************************************
916 
917 // clean stack for a new event
918 void stack::print(){
919  cout<<"<stack::print> stackId "<<stackId<<" Gains "<<gains<<" twrsInStack "<<twrsInStack<<" stackLoc "<<stackLoc<<" TwrsXY "<<nTwrsX<<"/"<<nTwrsY<<" mapUpdated "<<mapUpdated<<endl;
920  for (int itw = 0; itw < twrsInStack; itw++){
921  if(towers[itw]) {
922  towers[itw]->print();
923  }
924  }
925 }
926 
927 
928 
929 // **************************************************************************
930 void stack::updateMap(Int_t stckloc){
931  // print();
932  if(mapUpdated){
933  cout<<"<stack::updateMap> Already updated - EXIT"<<endl;
934  return;
935  }
936  stackLoc = stckloc;
937  for (int itw = 0; itw<twrsInStack; itw++){
938  towers[itw]->adcChannel[0] += stckloc;
939  if(gains==2) towers[itw]->adcChannel[1] += stckloc;
940  }
941  mapUpdated = kTRUE;
942 }
943 
944 
945 // **************************************************************************
946 
947 // collectdata from towers in stack for a new event
948 Int_t stack::update(Bool_t fitShape){
949  eventsSeen ++;
950  Int_t used = 0;
951  for (int itw = 0; itw < twrsInStack; itw++){
952  if(towers[itw]) {
953  Int_t rejected = towers[itw]->update(fitShape);
954  if(rejected<=2){
955  used++;
956  totPed += towers[itw]->cped;
957  totAmpl += towers[itw]->rampl;
958  totCorAmpl += towers[itw]->campl;
959  // with fit time used in every tower to identify event time
960  // all towers in the stack will now have identical times
961  avTwrTime += towers[itw]->ctime;
962  rmsTwrTime += towers[itw]->ctime*towers[itw]->ctime;
963  E += towers[itw]->E;
964  } // else if(rejected>5) towers[itw]->display();
965  for(int ig = 0; ig<gains; ig++){
966  satProb[ig] -> Fill(towers[itw]->satFlag[ig]);
967  }
968  }
969  }
970  if(!used) {
971  // everything in saturation
972  reject += 1000;
973  return reject;
974  }
975  avTwrTime /= used;
976  rmsTwrTime = 0.; // (twrsUsed==0? 0. : sqrt(rmsTwrTime/twrsUsed-avTwrTime*avTwrTime));
977  stSE -> Fill(E);
978  stSPed -> Fill(totPed);
979  stSAmpl -> Fill(totCorAmpl);
980  stAvT -> Fill(avTwrTime);
981  stTRMS -> Fill(rmsTwrTime);
982 
983  // Compute impact CG and find cosmic track
984  if(stackKind == CALOR) {
985  getStackImpact();
986  hLabHelper * hlHelper = hLabHelper ::getInstance();
987  if(hlHelper->runKind=="cosmics") getStackTrack();
988  }
989  return reject;
990 }
991 
992 
993 // **************************************************************************
994 
995 // getStackTime creates summary event (one per gain)
996 Int_t stack::getStackTime(){
997  hcalHelper * hHelper = hcalHelper::getInstance();
998  for(int ig = 0; ig < gains; ig++){
999  // 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