EMCAL/SMcalib/QA/LEDRef_evtdis.C

   1 ///////////////////////////////////////////////
   2 // Simple event display for channel spectra
   3 // Input file should be raw DATE root tree file
   4 ///////////////////////////////////////////////
   5
   6 // for the readout: 1 FEE reads out LED ref. info for the whole SuperModule
   7 // 1 CSP per StripModule
   8 const int NSTRIPS = 24; // number of StripModules in SuperModule
   9 const int TOTCHAN = NSTRIPS * 2;        // *2 since we have high gain and low gain
  10 // we group the calibrations in 3 sets of 8 strips
  11 const int NSETS = 3;
  12 const int NSTRIPS_IN_SET = 8;
  13
  14 // gamma2 fit function
  15 double fitfun(double *x, double *par) {
  16   double Amp    = par[0];
  17   double Tmax   = par[1];
  18   double Tau    = par[2];
  19   double Ped    = par[3];
  20   double gammaN   = par[4];
  21   double t = 0;
  22   if(Tau) t = (x[0] - Tmax + Tau)/Tau;
  23   if(t<0) t = 0;
  24
  25   // Since we have now set gammaN to 2, we replace the pow(t, 2) call in
  26   // double f = Amp * pow(t,gammaN) * exp(gammaN*(1-t)) + Ped;
  27   // with just t*t
  28   double f = Amp * t*t * exp(gammaN*(1-t)) + Ped;
  29   return f;
  30 }
  31
  32 // main method
  33 void LEDRef_evtdis(const int runno = 615,
  34                    const int gainv = 0,  /*0=low, 1=high*/
  35                    const int evtnum= -10,
  36                    int ymax=1023, // set the scale of plots
  37                    const int delay = 1)  // -1=no delay, wait for input, X>=0 => sleep aprox. X sec. after making plot
  38 {
  39   // set ranges to plot
  40   const int strip_f = 0; // first
  41   const int strip_l = NSTRIPS - 1;
  42
  43   const int nsamples = 65; // number of ADC time samples per channel and event
  44
  45   const int saveplot = 0;
  46   const int numbering      = 1; // 0: no numbering, 1: nubering on each plot
  47   const int dofit = 0; // 0: no fit, 1: try to fit the spectra
  48   const int debug    = 0;
  49   const float gammaN = 2;
  50   // end of setup
  51
  52   // Assume we are just interested in the 1st segment, _0.root below for fname*
  53   Char_t fname[256];
  54   sprintf(fname, "/local/data/Run_%09d.Seq_1A.Stream_0.root",runno);
  55   cout << "TOTCHAN " << TOTCHAN << endl;
  56
  57   // set up a raw reader of the data
  58   AliRawReader *rawReader = NULL;
  59   rawReader = new AliRawReaderRoot(fname);
  60   AliCaloRawStream *in = NULL;
  61   in = new AliCaloRawStream(rawReader,"EMCAL");
  62
  63   // set up histograms
  64   TH1F *hfit[TOTCHAN];
  65   TF1 *f1[TOTCHAN];
  66   char ch_label[TOTCHAN][100];
  67   char buff1[100];
  68   char name[80];
  69   for(int i=0; i<TOTCHAN; i++) {
  70     sprintf(buff1,"hfit_%d",i);
  71     hfit[i] = new TH1F(buff1,"hfit", nsamples , -0.5, nsamples - 0.5);
  72     hfit[i]->SetDirectory(0);
  73     sprintf(name,"f1_%d",i);
  74     f1[i] = new TF1(name,fitfun,0,70,5);
  75     f1[i]->SetLineWidth(2);
  76     f1[i]->SetLineColor(2);
  77
  78     //  int idx = istrip + NSTRIPS * gain; // encoding used later
  79     int gain = i / (NSTRIPS);
  80     int istrip = i % NSTRIPS;
  81     sprintf(ch_label[i], "Strip%02d", istrip);
  82   }
  83
  84   TCanvas *cc1 = new TCanvas("cc1","3 columns of 8 strips each",600,800);
  85   int numcol = NSETS;
  86   int numrow = NSTRIPS_IN_SET;
  87   cc1->Divide(numcol, numrow);
  88
  89   TText *t = new TText;
  90   t->SetTextSize(0.17);
  91   int clr[2] = {4,2}; // colors
  92
  93   // figure out which events we should look at
  94   int firstevent = evtnum;
  95   int lastevent = evtnum;
  96   if (evtnum < 0) { // get a bunch of events
  97     firstevent = 0;
  98     lastevent = - evtnum;
  99   }
 100   if (evtnum == 0) { // get all events
 101     firstevent = 0;
 102     lastevent = 1000000;
 103   }
 104
 105   Int_t iev =0;
 106   AliRawEventHeaderBase *aliHeader=NULL;
 107   while ( rawReader->NextEvent() && iev < firstevent) {
 108     aliHeader = (AliRawEventHeaderBase*) rawReader->GetEventHeader();
 109     iev++;
 110   }
 111
 112   // loop over selected events
 113   while ( rawReader->NextEvent() && iev <= lastevent) {
 114     aliHeader = (AliRawEventHeaderBase*) rawReader->GetEventHeader();
 115     int runNumber = aliHeader->Get("RunNb");
 116
 117     cout << "Found run number " << runNumber << endl;
 118
 119     // reset histograms
 120     for(int i=0; i<TOTCHAN; i++) {
 121       hfit[i]->Reset();
 122     }
 123
 124     // get events (the "1" ensures that we actually select all events for now)
 125     if ( 1 || aliHeader->Get("Type") == AliRawEventHeaderBase::kPhysicsEvent ) {
 126       const UInt_t * evtId = aliHeader->GetP("Id");
 127       int evno_raw = (int) evtId[0];
 128       int timestamp = aliHeader->Get("Timestamp");
 129
 130       cout << " evno " << evno_raw
 131            << " size " << aliHeader->GetEventSize()
 132            << " type " << aliHeader->Get("Type")
 133            << " type name " << aliHeader->GetTypeName()
 134            << " timestamp " << timestamp
 135            << endl;
 136
 137       /// process_event stream
 138       while ( in->Next() ) {
 139
 140         int strip = in->GetColumn();
 141         int gain = in->GetRow();
 142
 143         if (in->IsLEDMonData()) {
 144
 145           int idx = strip + NSTRIPS*gain;
 146           //cout << "hist idx " << idx << endl;
 147
 148           if (idx < 0 || idx > TOTCHAN) {
 149             cout << "Hist idx out of range: " << idx << endl;
 150           }
 151           else { // reasonable range of idx
 152             hfit[idx]->SetBinContent(in->GetTime(), in->GetSignal());
 153           }
 154
 155         } // LED Ref data only
 156
 157       } // Raw data read
 158
 159       // Next: let's actually plot the data..
 160       for (Int_t strip = strip_f; strip <= strip_l; strip++) {
 161
 162         int idx = strip + NSTRIPS*gainv;
 163
 164         // which set/column does the strip belong in
 165         int iset = strip / NSTRIPS_IN_SET;
 166         int within_set = strip % NSTRIPS_IN_SET;
 167         // on which pad should we plot it?
 168         int pad_id = (NSTRIPS_IN_SET-1-within_set)*NSETS + iset + 1;
 169
 170         cout << "strip " << strip
 171              << ". set="<< iset << ", within_set=" << within_set
 172              << ", pad=" << pad_id << endl;
 173         cc1->cd(pad_id);
 174         hfit[idx]->SetTitle("");
 175         hfit[idx]->SetFillColor(5);
 176         hfit[idx]->SetMaximum(ymax);
 177         hfit[idx]->SetMinimum(0);
 178         // we may or may not decide to fit the data
 179         if (dofit) {
 180           f1[i]->SetParameter(0, 0); // initial guess; zero amplitude :=)
 181           hfit[idx]->Fit(f1[i]);
 182         }
 183         hfit[idx]->Draw();
 184         if( numbering ) {
 185           t->SetTextColor(clr[gainv]);
 186           t->DrawTextNDC(0.65,0.65,ch_label[idx]);
 187         }
 188       }
 189
 190       // add some extra text on the canvas
 191       // print a box showing run #, evt #, and timestamp
 192       cc1->cd();
 193       // first draw transparent pad
 194       TPad *trans = new TPad("trans","",0,0,1,1);
 195       trans->SetFillStyle(4000);
 196       trans->Draw();
 197       trans->cd();
 198       // then draw text
 199       TPaveText *label = new TPaveText(.2,.11,.8,.14,"NDC");
 200       //  label->Clear();
 201       label->SetBorderSize(1);
 202       label->SetFillColor(0);
 203       label->SetLineColor(clr[gainv]);
 204       label->SetTextColor(clr[gainv]);
 205       //label->SetFillStyle(0);
 206       TDatime d;
 207       d.Set(timestamp);
 208       sprintf(name,"Run %d, Event %d, Hist Max %d, %s",runno,iev,ymax,d.AsString());
 209       label->AddText(name);
 210       label->Draw();
 211       cc1->Update();
 212       cout << "Done" << endl;
 213
 214       // some shenanigans to hold the plotting, if requested
 215       if (firstevent != lastevent) {
 216         if (delay == -1) {
 217           // wait for character input before proceeding
 218           cout << " enter y to proceed " << endl;
 219           char dummy[2];
 220           cin >> dummy;
 221           cout << " read " << dummy << endl;
 222           if (strcmp(dummy, "y")==0) {
 223             cout << " ok, continuing with event " << iev+1 << endl;
 224           }
 225           else {
 226             cout << " ok, exiting " << endl;
 227             //exit(1);
 228           }
 229         }
 230         else {
 231           cout << "Sleeping for " << delay * 500 << endl;
 232           gSystem->Sleep(delay * 500);
 233         }
 234       }
 235
 236       // save plot, if setup/requested to do so
 237       char plotname[100];
 238       if (saveplot==1) {
 239         sprintf(plotname,"Run_%d_LEDRef_Ev%d_Gain%d_MaxHist%d.gif",
 240                 runno,iev,gainv,ymax);
 241         cout <<"SAVING plot:"<< plotname << endl;
 242         cc1->SaveAs(plotname);
 243       }
 244
 245     } // event selection
 246
 247   iev ++;
 248   } // event loop
 249
 250 }