ZDC/ZDCLASERda.cxx

   1 /*
   2
   3 This program reads the DAQ data files passed as argument using the monitoring library.
   4
   5 It computes the average event size and populates local "./result.txt" file with the
   6 result.
   7
   8 The program reports about its processing progress.
   9
  10 Messages on stdout are exported to DAQ log system.
  11
  12 DA for ZDC standalone pedestal runs
  13
  14 Contact: Chiara.Oppedisano@to.infn.it
  15 Link:
  16 Run Type: STANDALONE_LASER_RUN
  17 DA Type: LDC
  18 Number of events needed: no constraint (tipically ~10^3)
  19 Input Files:
  20 Output Files: ZDCLaser.dat
  21 Trigger Types Used: Standalone Trigger
  22
  23 */
  24 #define PEDDATA_FILE  "ZDCPedestal.dat"
  25 #define MAPDATA_FILE  "ZDCChMapping.dat"
  26 #define LASHISTO_FILE "ZDCLaserHisto.root"
  27 #define LASDATA_FILE  "ZDCLaserCalib.dat"
  28
  29 #include <stdio.h>
  30 #include <stdlib.h>
  31 #include <Riostream.h>
  32
  33 // DATE
  34 #include <event.h>
  35 #include <monitor.h>
  36 #include <daqDA.h>
  37
  38 //ROOT
  39 #include <TROOT.h>
  40 #include <TPluginManager.h>
  41 #include <TH1F.h>
  42 #include <TF1.h>
  43 #include <TFile.h>
  44 #include <TFitter.h>
  45 #include "TMinuitMinimizer.h"
  46
  47 //AliRoot
  48 #include <AliRawReaderDate.h>
  49 #include <AliRawEventHeaderBase.h>
  50 #include <AliZDCRawStream.h>
  51
  52
  53 /* Main routine
  54       Arguments: list of DATE raw data files
  55 */
  56 int main(int argc, char **argv) {
  57
  58   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
  59                                         "*",
  60                                         "TStreamerInfo",
  61                                         "RIO",
  62                                         "TStreamerInfo()");
  63
  64   TMinuitMinimizer m;
  65   gROOT->GetPluginManager()->AddHandler("ROOT::Math::Minimizer", "Minuit","TMinuitMinimizer",
  66       "Minuit", "TMinuitMinimizer(const char *)");
  67   TVirtualFitter::SetDefaultFitter("Minuit");
  68
  69   int status = 0;
  70   int const kNChannels = 24;
  71   int const kNScChannels = 32;
  72
  73   /* log start of process */
  74   printf("\n ZDC LASER program started\n");
  75
  76   /* check that we got some arguments = list of files */
  77   if (argc<2) {
  78     printf("Wrong number of arguments\n");
  79     return -1;
  80   }
  81
  82   // --- Histograms for LASER runs
  83   //     20 signal channels + 2 reference PTMs
  84   //
  85   TH1F::AddDirectory(0);
  86   // --- Histos for reference PMTs (high gain chains)
  87   TH1F *hPMRefChg = new TH1F("hPMRefChg","hPMRefChg", 100,0.,1400.);
  88   TH1F *hPMRefAhg = new TH1F("hPMRefAhg","hPMRefAhg", 100,0.,1400.);
  89   TH1F *hPMRefClg = new TH1F("hPMRefClg","hPMRefClg", 100,0.,4000.);
  90   TH1F *hPMRefAlg = new TH1F("hPMRefAlg","hPMRefAlg", 100,0.,4000.);
  91   //
  92   // --- Histos for detector PMTs
  93   TH1F *hZNChg[5], *hZPChg[5], *hZNAhg[5], *hZPAhg[5], *hZEMhg[2];
  94   TH1F *hZNClg[5], *hZPClg[5], *hZNAlg[5], *hZPAlg[5], *hZEMlg[2];
  95   char hnamZNChg[20], hnamZPChg[20], hnamZNAhg[20], hnamZPAhg[20];
  96   char hnamZNClg[20], hnamZPClg[20], hnamZNAlg[20], hnamZPAlg[20];
  97   char hnamZEMhg[20], hnamZEMlg[20];
  98   for(Int_t j=0; j<5; j++){
  99     sprintf(hnamZNChg,"ZNChg-tow%d",j);
 100     sprintf(hnamZPChg,"ZPChg-tow%d",j);
 101     sprintf(hnamZNAhg,"ZNAhg-tow%d",j);
 102     sprintf(hnamZPAhg,"ZPAhg-tow%d",j);
 103     //
 104     hZNChg[j] = new TH1F(hnamZNChg, hnamZNChg, 100, 0., 1400.);
 105     hZPChg[j] = new TH1F(hnamZPChg, hnamZPChg, 100, 0., 1400.);
 106     hZNAhg[j] = new TH1F(hnamZNAhg, hnamZNAhg, 100, 0., 1400.);
 107     hZPAhg[j] = new TH1F(hnamZPAhg, hnamZPAhg, 100, 0., 1400.);
 108     //
 109     sprintf(hnamZNClg,"ZNClg-tow%d",j);
 110     sprintf(hnamZPClg,"ZPClg-tow%d",j);
 111     sprintf(hnamZNAlg,"ZNAlg-tow%d",j);
 112     sprintf(hnamZPAlg,"ZPAlg-tow%d",j);
 113     //
 114     hZNClg[j] = new TH1F(hnamZNClg, hnamZNClg, 100, 0., 4000.);
 115     hZPClg[j] = new TH1F(hnamZPClg, hnamZPClg, 100, 0., 4000.);
 116     hZNAlg[j] = new TH1F(hnamZNAlg, hnamZNAlg, 100, 0., 4000.);
 117     hZPAlg[j] = new TH1F(hnamZPAlg, hnamZPAlg, 100, 0., 4000.);
 118     //
 119     if(j<2){
 120       sprintf(hnamZEMhg,"ZEM%dhg",j);
 121       sprintf(hnamZEMlg,"ZEM%dlg",j);
 122       //
 123       hZEMhg[j] = new TH1F(hnamZEMhg, hnamZEMhg, 100, 0., 1400.);
 124       hZEMlg[j] = new TH1F(hnamZEMlg, hnamZEMlg, 100, 0., 4000.);
 125     }
 126   }
 127
 128   /* open result file */
 129   FILE *fp=NULL;
 130   fp=fopen("./result.txt","a");
 131   if (fp==NULL) {
 132     printf("Failed to open file\n");
 133     return -1;
 134   }
 135   /* report progress */
 136   daqDA_progressReport(10);
 137
 138   // *** To analyze LASER events you MUST have a pedestal data file!!!
 139   // *** -> check if a pedestal run has been analyzed
 140   int read = 0;
 141   read = daqDA_DB_getFile(PEDDATA_FILE, PEDDATA_FILE);
 142   if(read){
 143     printf("\t ERROR!!! ZDCPedestal.dat file NOT FOUND in DAQ db!!!\n");
 144     return -1;
 145   }
 146   else printf("\t ZDCPedestal.dat file retrieved from DAQ db\n");
 147
 148   FILE *filePed = fopen(PEDDATA_FILE,"r");
 149   if (filePed==NULL) {
 150     printf("\t ERROR!!! Can't open ZDCPedestal.dat file!!!\n");
 151     return -1;
 152   }
 153
 154   // 144 = 48 in-time + 48 out-of-time + 48 correlations
 155   Float_t readValues[2][6*kNChannels];
 156   Float_t MeanPedhg[kNChannels], MeanPedlg[kNChannels];
 157   Float_t CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
 158   // ***************************************************
 159   //   Unless we have a narrow correlation to fit we
 160   //    don't fit and store in-time vs. out-of-time
 161   //    histograms -> mean pedstal subtracted!!!!!!
 162   // ***************************************************
 163   //
 164   for(int jj=0; jj<6*kNChannels; jj++){
 165     for(int ii=0; ii<2; ii++){
 166        fscanf(filePed,"%f",&readValues[ii][jj]);
 167     }
 168     if(jj<kNChannels){
 169       MeanPedhg[jj] = readValues[0][jj];
 170       //printf("\t MeanPedhg[%d] = %1.1f\n",jj, MeanPedhg[jj]);
 171     }
 172     else if(jj>=kNChannels && jj<2*kNChannels){
 173       MeanPedlg[jj-kNChannels] = readValues[0][jj];
 174       //printf("\t MeanPedlg[%d] = %1.1f\n",jj-kNChannels, MeanPedlg[jj-kNChannels]);
 175     }
 176     else if(jj>4*kNChannels){
 177       CorrCoeff0[jj-4*kNChannels] = readValues[0][jj];
 178       CorrCoeff1[jj-4*kNChannels] = readValues[1][jj];;
 179     }
 180   }
 181
 182   FILE *mapFile4Shuttle;
 183
 184   /* report progress */
 185   daqDA_progressReport(20);
 186
 187
 188   /* init some counters */
 189   int nevents_physics=0;
 190   int nevents_total=0;
 191
 192   /* read the data files */
 193   int n;
 194   for(n=1;n<argc;n++) {
 195
 196     status=monitorSetDataSource( argv[n] );
 197     if (status!=0) {
 198       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
 199       return -1;
 200     }
 201
 202     /* report progress */
 203     /* in this example, indexed on the number of files */
 204     daqDA_progressReport(20+70*n/argc);
 205
 206     /* read the file */
 207     for(;;) {
 208       struct eventHeaderStruct *event;
 209       eventTypeType eventT;
 210
 211       /* get next event */
 212       status=monitorGetEventDynamic((void **)&event);
 213       if (status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
 214       if (status!=0) {
 215         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 216         return -1;
 217       }
 218
 219       /* retry if got no event */
 220       if (event==NULL) {
 221         break;
 222       }
 223
 224       // Initalize raw-data reading and decoding
 225       AliRawReader *reader = new AliRawReaderDate((void*)event);
 226       reader->Select("ZDC");
 227       // --- Reading event header
 228       //UInt_t evtype = reader->GetType();
 229       //printf("\n\t ZDCLASERda -> ev. type %d\n",evtype);
 230       //printf("\t ZDCLASERda -> run # %d\n",reader->GetRunNumber());
 231       //
 232       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
 233
 234
 235       /* use event - here, just write event id to result file */
 236       eventT=event->eventType;
 237
 238
 239       Int_t ich=0;
 240       Int_t adcMod[2*kNChannels], adcCh[2*kNChannels], sigCode[2*kNChannels];
 241       Int_t det[2*kNChannels], sec[2*kNChannels];
 242       for(Int_t y=0; y<2*kNChannels; y++){
 243         adcMod[y]=adcCh[y]=sigCode[y]=det[y]=sec[y]=0;
 244       }
 245
 246       Int_t iScCh=0;
 247       Int_t scMod[kNScChannels], scCh[kNScChannels], scSigCode[kNScChannels];
 248       Int_t scDet[kNScChannels], scSec[kNScChannels];
 249       for(Int_t y=0; y<kNScChannels; y++){
 250         scMod[y]=scCh[y]=scSigCode[y]=scDet[y]=scSec[y]=0;
 251       }
 252       //
 253       Int_t modNum=-1, modType=-1;
 254
 255       if(eventT==START_OF_DATA){
 256
 257         rawStreamZDC->SetSODReading(kTRUE);
 258
 259         // --------------------------------------------------------
 260         // --- Writing ascii data file for the Shuttle preprocessor
 261         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
 262         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 263         else{
 264           while((rawStreamZDC->Next())){
 265             if(rawStreamZDC->IsHeaderMapping()){ // mapping header
 266                modNum = rawStreamZDC->GetADCModule();
 267                modType = rawStreamZDC->GetModType();
 268             }
 269             if(rawStreamZDC->IsChMapping()){
 270               if(modType==1){ // ADC mapping ----------------------
 271                 adcMod[ich]  = rawStreamZDC->GetADCModFromMap(ich);
 272                 adcCh[ich]   = rawStreamZDC->GetADCChFromMap(ich);
 273                 sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
 274                 det[ich]     = rawStreamZDC->GetDetectorFromMap(ich);
 275                 sec[ich]     = rawStreamZDC->GetTowerFromMap(ich);
 276                 //
 277                 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",
 278                   ich,adcMod[ich],adcCh[ich],sigCode[ich],det[ich],sec[ich]);
 279                 //
 280                 //printf("  Mapping in DA -> %d ADC: mod %d ch %d, code %d det %d, sec %d\n",
 281                 //  ich,adcMod[ich],adcCh[ich],sigCode[ich],det[ich],sec[ich]);
 282                 //
 283                 ich++;
 284               }
 285               else if(modType==2){ //VME scaler mapping --------------------
 286                 scMod[iScCh]     = rawStreamZDC->GetScalerModFromMap(iScCh);
 287                 scCh[iScCh]      = rawStreamZDC->GetScalerChFromMap(iScCh);
 288                 scSigCode[iScCh] = rawStreamZDC->GetScalerSignFromMap(iScCh);
 289                 scDet[iScCh]     = rawStreamZDC->GetScDetectorFromMap(iScCh);
 290                 scSec[iScCh]    = rawStreamZDC->GetScTowerFromMap(iScCh);
 291                 //
 292                 fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",
 293                   iScCh,scMod[iScCh],scCh[iScCh],scSigCode[iScCh],scDet[iScCh],scSec[iScCh]);
 294                 //
 295                 //printf("  Mapping in DA -> %d Scaler: mod %d ch %d, code %d det %d, sec %d\n",
 296                 //  iScCh,scMod[iScCh],scCh[iScCh],scSigCode[iScCh],scDet[iScCh],scSec[iScCh]);
 297                 //
 298                 iScCh++;
 299               }
 300             }
 301           }
 302         }
 303         fclose(mapFile4Shuttle);
 304       }// SOD event
 305
 306       if(eventT==PHYSICS_EVENT){
 307         // --- Reading data header
 308         reader->ReadHeader();
 309         const AliRawDataHeader* header = reader->GetDataHeader();
 310         if(header) {
 311          UChar_t message = header->GetAttributes();
 312          if(message & 0x30){ // DEDICATED LASER RUN
 313             //printf("\t STANDALONE_LASER_RUN raw data found\n");
 314             continue;
 315          }
 316          else{
 317             printf("ZDCLASERda.cxx -> NO STANDALONE_LASER_RUN raw data found\n");
 318             return -1;
 319          }
 320         }
 321         else{
 322            printf("\t ATTENTION! No Raw Data Header found!!!\n");
 323            return -1;
 324         }
 325
 326         rawStreamZDC->SetSODReading(kTRUE);
 327
 328         if (!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 329         //
 330         // ----- Setting ch. mapping -----
 331         for(Int_t jk=0; jk<2*kNChannels; jk++){
 332           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
 333           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
 334           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
 335           rawStreamZDC->SetMapDet(jk, det[jk]);
 336           rawStreamZDC->SetMapTow(jk, sec[jk]);
 337         }
 338         //
 339         while(rawStreamZDC->Next()){
 340           Int_t index=-1;
 341           Int_t detector = rawStreamZDC->GetSector(0);
 342           Int_t sector = rawStreamZDC->GetSector(1);
 343
 344           if(rawStreamZDC->IsADCDataWord() && !(rawStreamZDC->IsUnderflow())
 345              && !(rawStreamZDC->IsOverflow()) && detector!=-1){
 346
 347             //printf("  IsADCWord %d, IsUnderflow %d, IsOverflow %d\n",
 348             //  rawStreamZDC->IsADCDataWord(),rawStreamZDC->IsUnderflow(),rawStreamZDC->IsOverflow());
 349
 350             if(sector!=5){ // Physics signals
 351               if(detector==1) index = sector;        // *** ZNC
 352               else if(detector==2) index = sector+5; // *** ZPC
 353               else if(detector==3) index = sector+9; // *** ZEM
 354               else if(detector==4) index = sector+12;// *** ZNA
 355               else if(detector==5) index = sector+17;// *** ZPA
 356             }
 357             else{ // Reference PMs
 358               index = (detector-1)/3+22;
 359             }
 360             //
 361             if(index==-1) printf("ERROR in ZDCLASERda.cxx -> det %d quad %d res %d index %d ADC %d\n",
 362               detector, sector, rawStreamZDC->GetADCGain(), index, rawStreamZDC->GetADCValue());
 363
 364             Float_t Pedestal=0.;
 365             if(rawStreamZDC->GetADCGain()==0) Pedestal = MeanPedhg[index];
 366             else if(rawStreamZDC->GetADCGain()==1) Pedestal = MeanPedlg[index];
 367             //
 368             Float_t CorrADC = rawStreamZDC->GetADCValue() - Pedestal;
 369             //
 370             //printf("\tdet %d sec %d res %d index %d ped %1.0f ADCcorr %1.0f\n",
 371             //  detector, sector, rawStreamZDC->GetADCGain(), index, Pedestal,CorrADC);
 372
 373             // **** Detector PMs
 374             if(sector!=5){
 375               if(rawStreamZDC->GetADCGain()==0){ // --- High gain chain ---
 376                 // ---- side C
 377                 if(detector==1) hZNChg[sector]->Fill(CorrADC);
 378                 else if(detector==2) hZPChg[sector]->Fill(CorrADC);
 379                 // ---- side A
 380                 else if(detector==4) hZNAhg[sector]->Fill(CorrADC);
 381                 else if(detector==5) hZPAhg[sector]->Fill(CorrADC);
 382                 // ---- ZEM
 383                 else if(detector==3) hZEMhg[sector-1]->Fill(CorrADC);
 384               }
 385               else if(rawStreamZDC->GetADCGain()==1){ // --- Low gain chain ---
 386                 // ---- side C
 387                 if(detector==1) hZNClg[sector]->Fill(CorrADC);
 388                 else if(detector==2) hZPClg[sector]->Fill(CorrADC);
 389                 // ---- side A
 390                 else if(detector==4) hZNAlg[sector]->Fill(CorrADC);
 391                 else if(detector==5) hZPAlg[sector]->Fill(CorrADC);
 392                 // ---- ZEM
 393                 else if(detector==3) hZEMlg[sector-1]->Fill(CorrADC);
 394               }
 395             }
 396             // **** Reference PMs
 397             else if(sector==5){
 398               if(rawStreamZDC->GetADCGain()==0){ // --- High gain chain ---
 399                 // ---- PMRef chain side C
 400                 if(detector==1) hPMRefChg->Fill(CorrADC);
 401                 // ---- PMRef side A
 402                 else if(detector==4) hPMRefAhg->Fill(CorrADC);
 403               }
 404               else if(rawStreamZDC->GetADCGain()==1){ // --- Low gain chain ---
 405                 // ---- PMRef chain side C
 406                 if(detector==1) hPMRefClg->Fill(CorrADC);
 407                 // ---- PMRef side A
 408                 else if(detector==4) hPMRefAlg->Fill(CorrADC);
 409               }
 410             }
 411           }//IsADCDataWord()+NOunderflow+NOoverflow
 412           //
 413          }
 414          //
 415          nevents_physics++;
 416          //
 417          delete reader;
 418          delete rawStreamZDC;
 419
 420       }//(if PHYSICS_EVENT)
 421       nevents_total++;
 422
 423       /* free resources */
 424       free(event);
 425
 426     }
 427   }
 428
 429   /* Analysis of the histograms */
 430   //
 431   Int_t det[2*kNChannels], quad[2*kNChannels];
 432   Int_t maxBin[2*kNChannels], nBin[2*kNChannels];
 433   Float_t xMax[2*kNChannels], maxXval[2*kNChannels], xlow[2*kNChannels];
 434   Float_t mean[2*kNChannels], sigma[2*kNChannels];
 435   TF1 *fun[2*kNChannels];
 436
 437   // ******** High gain chain ********
 438   for(Int_t k=0; k<5; k++){
 439     // --- ZNC
 440     det[k] = 1;
 441     quad[k] = k;
 442     maxBin[k] = hZNChg[k]->GetMaximumBin();
 443     nBin[k] = (hZNChg[k]->GetXaxis())->GetNbins();
 444     xMax[k] = (hZNChg[k]->GetXaxis())->GetXmax();
 445     if(nBin[k]!=0) maxXval[k] = maxBin[k]*xMax[k]/nBin[k];
 446     if(maxXval[k]-150.<0.) xlow[k]=0.;
 447     else xlow[k] = maxXval[k]-150.;
 448     // checking if histos are empty
 449     if(hZNChg[k]->GetEntries() == 0){
 450       printf("\n WARNING! Empty LASER histos -> ending DA WITHOUT writing output\n\n");
 451       return -1;
 452     }
 453     //
 454     hZNChg[k]->Fit("gaus","Q","",xlow[k],maxXval[k]+150.);
 455     fun[k] = hZNChg[k]->GetFunction("gaus");
 456     mean[k]  = (Float_t) (fun[k]->GetParameter(1));
 457     sigma[k] = (Float_t) (fun[k]->GetParameter(2));
 458     // --- ZPC
 459     det[k+5] = 2;
 460     quad[k+5] = k;
 461     maxBin[k+5] = hZPChg[k]->GetMaximumBin();
 462     nBin[k+5] = (hZPChg[k]->GetXaxis())->GetNbins();
 463     xMax[k+5] = (hZPChg[k]->GetXaxis())->GetXmax();
 464     if(nBin[k+5]!=0) maxXval[k+5] = maxBin[k+5]*xMax[k+5]/nBin[k+5];
 465     if(maxXval[k+5]-150.<0.) xlow[k+5]=0.;
 466     else xlow[k+5] = maxXval[k+5]-150.;
 467     hZPChg[k]->Fit("gaus","Q","",xlow[k+5],maxXval[k+5]+150.);
 468     fun[k+5] = hZPChg[k]->GetFunction("gaus");
 469     mean[k+5]  = (Float_t) (fun[k+5]->GetParameter(1));
 470     sigma[k+5] = (Float_t) (fun[k+5]->GetParameter(2));
 471     // --- ZEM1
 472     if(k<2){
 473       det[k+10] = 3;
 474       quad[k+10] = k+1;
 475       maxBin[k+10] = hZEMhg[k]->GetMaximumBin();
 476       nBin[k+10] = (hZEMhg[k]->GetXaxis())->GetNbins();
 477       xMax[k+10] = (hZEMhg[k]->GetXaxis())->GetXmax();
 478       if(nBin[k+10]!=0) maxXval[k+10] = maxBin[k+10]*xMax[k+10]/nBin[k+10];
 479       if(maxXval[k+10]-150.<0.) xlow[k+10]=0.;
 480       else xlow[k+10] = maxXval[k+10]-150.;
 481       hZEMhg[k]->Fit("gaus","Q","",xlow[k+10],maxXval[k+10]+150.);
 482       fun[k+10] = hZEMhg[k]->GetFunction("gaus");
 483       mean[k+10]  = (Float_t) (fun[k+10]->GetParameter(1));
 484       sigma[k+10] = (Float_t) (fun[k+10]->GetParameter(2));
 485     }
 486     // --- ZNA
 487     det[k+12] = 4;
 488     quad[k+12] = k;
 489     maxBin[k+12] = hZNAhg[k]->GetMaximumBin();
 490     nBin[k+12] = (hZNAhg[k]->GetXaxis())->GetNbins();
 491     xMax[k+12] = (hZNAhg[k]->GetXaxis())->GetXmax();
 492     if(nBin[k+12]!=0) maxXval[k+12] = maxBin[k+12]*xMax[k+12]/nBin[k+12];
 493     if(maxXval[k+12]-150.<0.) xlow[k+12]=0.;
 494     else xlow[k+12] = maxXval[k+12]-150.;
 495     hZNAhg[k]->Fit("gaus","Q","",xlow[k+12],maxXval[k+12]+150.);
 496     fun[k+12] = hZNAhg[k]->GetFunction("gaus");
 497     mean[k+12]  = (Float_t) (fun[k+12]->GetParameter(1));
 498     sigma[k+12] = (Float_t) (fun[k+12]->GetParameter(2));
 499     // --- ZPA
 500     det[k+17] = 4;
 501     quad[k+17] = 5;
 502     maxBin[k+17] = hZPAhg[k]->GetMaximumBin();
 503     nBin[k+17] = (hZPAhg[k]->GetXaxis())->GetNbins();
 504     xMax[k+17] = (hZPAhg[k]->GetXaxis())->GetXmax();
 505     if(nBin[k+17]!=0) maxXval[k+17] = maxBin[k+17]*xMax[k+17]/nBin[k+17];
 506     if(maxXval[k+17]-150.<0.) xlow[k+17]=0.;
 507     else xlow[k+17] = maxXval[k+17]-150.;
 508     hZPAhg[k]->Fit("gaus","Q","",xlow[k+17],maxXval[k+17]+150.);
 509     fun[k+17] = hZPAhg[k]->GetFunction("gaus");
 510     mean[k+17]  = (Float_t) (fun[k+17]->GetParameter(1));
 511     sigma[k+17] = (Float_t) (fun[k+17]->GetParameter(2));
 512   }
 513   // ~~~~~~~~ PM Ref side C ~~~~~~~~
 514   det[22] = 1;
 515   quad[22] = 5;
 516   maxBin[22] = hPMRefChg->GetMaximumBin();
 517   nBin[22] = (hPMRefChg->GetXaxis())->GetNbins();
 518   xMax[22] = (hPMRefChg->GetXaxis())->GetXmax();
 519   if(nBin[22]!=0) maxXval[22] = maxBin[22]*xMax[22]/nBin[22];
 520   if(maxXval[22]-150.<0.) xlow[22]=0.;
 521   else xlow[22] = maxXval[22];
 522   hPMRefChg->Fit("gaus","Q","",xlow[22],maxXval[22]+150.);
 523   fun[22] = hPMRefChg->GetFunction("gaus");
 524   mean[22]  = (Float_t) (fun[22]->GetParameter(1));
 525   sigma[22] = (Float_t) (fun[22]->GetParameter(2));
 526   // ~~~~~~~~ PM Ref side A ~~~~~~~~
 527   det[23] = 4;
 528   quad[23] = 5;
 529   maxBin[23] = hPMRefAhg->GetMaximumBin();
 530   nBin[23] = (hPMRefAhg->GetXaxis())->GetNbins();
 531   xMax[23] = (hPMRefAhg->GetXaxis())->GetXmax();
 532   if(nBin[23]!=0) maxXval[23] = maxBin[23]*xMax[23]/nBin[23];
 533   if(maxXval[23]-100.<0.) xlow[23]=0.;
 534   else xlow[23] = maxXval[23];
 535   hPMRefAhg->Fit("gaus","Q","",xlow[23],maxXval[23]+100.);
 536   fun[23] = hPMRefAhg->GetFunction("gaus");
 537   mean[23]  = (Float_t) (fun[23]->GetParameter(1));
 538   sigma[23] = (Float_t) (fun[23]->GetParameter(2));
 539
 540   // ******** Low gain chain ********
 541   Int_t kOffset = 24;
 542   for(Int_t k=0; k<5; k++){
 543     // --- ZNC
 544     det[k+kOffset] = 1;
 545     quad[k+kOffset] = k;
 546     maxBin[k+kOffset] = hZNClg[k]->GetMaximumBin();
 547     nBin[k+kOffset] = (hZNClg[k]->GetXaxis())->GetNbins();
 548     xMax[k+kOffset] = (hZNClg[k]->GetXaxis())->GetXmax();
 549     if(nBin[k+kOffset]!=0) maxXval[k+kOffset] = maxBin[k+kOffset]*xMax[k+kOffset]/nBin[k+kOffset];
 550     if(maxXval[k+kOffset]-150.<0.) xlow[k+kOffset]=0.;
 551     else xlow[k+kOffset] = maxXval[k+kOffset]-150.;
 552     hZNClg[k]->Fit("gaus","Q","",xlow[k+kOffset],maxXval[k+kOffset]+150.);
 553     fun[k+kOffset] = hZNClg[k]->GetFunction("gaus");
 554     mean[k+kOffset]  = (Float_t) (fun[k+kOffset]->GetParameter(1));
 555     sigma[k+kOffset] = (Float_t) (fun[k+kOffset]->GetParameter(2));
 556     // --- ZPC
 557     det[k+kOffset+5] = 2;
 558     quad[k+kOffset+5] = k;
 559     maxBin[k+kOffset+5] = hZPClg[k]->GetMaximumBin();
 560     nBin[k+kOffset+5] = (hZPClg[k]->GetXaxis())->GetNbins();
 561     xMax[k+kOffset+5] = (hZPClg[k]->GetXaxis())->GetXmax();
 562     if(nBin[k+kOffset+5]!=0) maxXval[k+kOffset+5] = maxBin[k+kOffset+5]*xMax[k+kOffset+5]/nBin[k+kOffset+5];
 563     if(maxXval[k+kOffset+5]-150.<0.) xlow[k+kOffset+5]=0.;
 564     else xlow[k+kOffset+5] = maxXval[k+kOffset+5]-150.;
 565     hZPClg[k]->Fit("gaus","Q","",xlow[k+kOffset+5],maxXval[k+kOffset+5]+150.);
 566     fun[k+kOffset+5] = hZPClg[k]->GetFunction("gaus");
 567     mean[k+kOffset+5]  = (Float_t) (fun[k+kOffset+5]->GetParameter(1));
 568     sigma[k+kOffset+5] = (Float_t) (fun[k+kOffset+5]->GetParameter(2));
 569     // --- ZEM1
 570     if(k+kOffset<2){
 571       det[k+kOffset+10] = 3;
 572       quad[k+kOffset+10] = k+1;
 573       maxBin[k+kOffset+10] = hZEMlg[k]->GetMaximumBin();
 574       nBin[k+kOffset+10] = (hZEMlg[k]->GetXaxis())->GetNbins();
 575       xMax[k+kOffset+10] = (hZEMlg[k]->GetXaxis())->GetXmax();
 576       if(nBin[k+kOffset+10]!=0) maxXval[k+kOffset+10] = maxBin[k+kOffset+10]*xMax[k+kOffset+10]/nBin[k+kOffset+10];
 577       if(maxXval[k+kOffset+10]-150.<0.) xlow[k+kOffset+10]=0.;
 578       else xlow[k+kOffset+10] = maxXval[k+kOffset+10]-150.;
 579       hZEMlg[k]->Fit("gaus","Q","",xlow[k+kOffset+10],maxXval[k+kOffset+10]+150.);
 580       fun[k+kOffset+10] = hZEMlg[k]->GetFunction("gaus");
 581       mean[k+kOffset+10]  = (Float_t) (fun[k+kOffset+10]->GetParameter(1));
 582       sigma[k+kOffset+10] = (Float_t) (fun[k+kOffset+10]->GetParameter(2));
 583     }
 584     // --- ZNA
 585     det[k+kOffset+12] = 4;
 586     quad[k+kOffset+12] = k;
 587     maxBin[k+kOffset+12] = hZNAlg[k]->GetMaximumBin();
 588     nBin[k+kOffset+12] = (hZNAlg[k]->GetXaxis())->GetNbins();
 589     xMax[k+kOffset+12] = (hZNAlg[k]->GetXaxis())->GetXmax();
 590     if(nBin[k+kOffset+12]!=0) maxXval[k+kOffset+12] = maxBin[k+kOffset+12]*xMax[k+kOffset+12]/nBin[k+kOffset+12];
 591     if(maxXval[k+kOffset+12]-150.<0.) xlow[k+kOffset+12]=0.;
 592     else xlow[k+kOffset+12] = maxXval[k+kOffset+12]-150.;
 593     hZNAlg[k]->Fit("gaus","Q","",xlow[k+kOffset+12],maxXval[k+kOffset+12]+150.);
 594     fun[k+kOffset+12] = hZNAlg[k]->GetFunction("gaus");
 595     mean[k+kOffset+12]  = (Float_t) (fun[k+kOffset+12]->GetParameter(1));
 596     sigma[k+kOffset+12] = (Float_t) (fun[k+kOffset+12]->GetParameter(2));
 597     // --- ZPA
 598     det[k+kOffset+17] = 5;
 599     quad[k+kOffset+17] = k;
 600     maxBin[k+kOffset+17] = hZPAlg[k]->GetMaximumBin();
 601     nBin[k+kOffset+17] = (hZPAlg[k]->GetXaxis())->GetNbins();
 602     xMax[k+kOffset+17] = (hZPAlg[k]->GetXaxis())->GetXmax();
 603     if(nBin[k+kOffset+17]!=0) maxXval[k+kOffset+17] = maxBin[k+kOffset+17]*xMax[k+kOffset+17]/nBin[k+kOffset+17];
 604     if(maxXval[k+kOffset+17]-150.<0.) xlow[k+kOffset+17]=0.;
 605     else xlow[k+kOffset+17] = maxXval[k+kOffset+17]-150.;
 606     hZPAlg[k]->Fit("gaus","Q","",xlow[k+kOffset+17],maxXval[k+kOffset+17]+150.);
 607     fun[k+kOffset+17] = hZPAlg[k]->GetFunction("gaus");
 608     mean[k+kOffset+17]  = (Float_t) (fun[k+kOffset+17]->GetParameter(1));
 609     sigma[k+kOffset+17] = (Float_t) (fun[k+kOffset+17]->GetParameter(2));
 610   }
 611   // ~~~~~~~~ PM Ref side C ~~~~~~~~
 612   det[46] = 1;
 613   quad[46] = 5;
 614   maxBin[46] = hPMRefClg->GetMaximumBin();
 615   nBin[46] = (hPMRefClg->GetXaxis())->GetNbins();
 616   xMax[46] = (hPMRefClg->GetXaxis())->GetXmax();
 617   if(nBin[46]!=0) maxXval[46] = maxBin[46]*xMax[46]/nBin[46];
 618   if(maxXval[46]-150.<0.) xlow[46]=0.;
 619   else xlow[46] = maxXval[46];
 620   hPMRefClg->Fit("gaus","Q","",xlow[46],maxXval[46]+150.);
 621   fun[46] = hPMRefClg->GetFunction("gaus");
 622   mean[46]  = (Float_t) (fun[46]->GetParameter(1));
 623   sigma[46] = (Float_t) (fun[46]->GetParameter(2));
 624   // ~~~~~~~~ PM Ref side A ~~~~~~~~
 625   det[47] = 4;
 626   quad[47] = 5;
 627   maxBin[47] = hPMRefAlg->GetMaximumBin();
 628   nBin[47] = (hPMRefAlg->GetXaxis())->GetNbins();
 629   xMax[47] = (hPMRefAlg->GetXaxis())->GetXmax();
 630   if(nBin[47]!=0) maxXval[47] = maxBin[47]*xMax[47]/nBin[47];
 631   if(maxXval[47]-100.<0.) xlow[47]=0.;
 632   else xlow[47] = maxXval[47];
 633   hPMRefAlg->Fit("gaus","Q","",xlow[47],maxXval[47]+100.);
 634   fun[47] = hPMRefAlg->GetFunction("gaus");
 635   mean[47]  = (Float_t) (fun[47]->GetParameter(1));
 636   sigma[47] = (Float_t) (fun[47]->GetParameter(2));
 637
 638   FILE *fileShuttle;
 639   fileShuttle = fopen(LASDATA_FILE,"w");
 640   for(Int_t i=0; i<2*kNChannels; i++){
 641     fprintf(fileShuttle,"\t%d\t%d\t%f\t%f\n",det[i],quad[i],mean[i], sigma[i]);
 642   }
 643   //
 644   fclose(fileShuttle);
 645   /* report progress */
 646   daqDA_progressReport(80);
 647   //
 648   TFile *histofile = new TFile(LASHISTO_FILE,"RECREATE");
 649   histofile->cd();
 650   for(int j=0; j<5; j++){
 651      hZNChg[j]->Write();
 652      hZPChg[j]->Write();
 653      hZNAhg[j]->Write();
 654      hZPAhg[j]->Write();
 655      hZNClg[j]->Write();
 656      hZPClg[j]->Write();
 657      hZNAlg[j]->Write();
 658      hZPAlg[j]->Write();
 659      if(j<2){
 660        hZEMhg[j]->Write();
 661        hZEMlg[j]->Write();
 662     }
 663   }
 664   hPMRefChg->Write();
 665   hPMRefAhg->Write();
 666   hPMRefClg->Write();
 667   hPMRefAlg->Write();
 668   //
 669   histofile->Close();
 670   //
 671   for(Int_t j=0; j<5; j++){
 672     delete hZNChg[j];
 673     delete hZPChg[j];
 674     delete hZNAhg[j];
 675     delete hZPAhg[j];
 676     delete hZNClg[j];
 677     delete hZPClg[j];
 678     delete hZNAlg[j];
 679     delete hZPAlg[j];
 680     if(j<2){
 681       delete hZEMhg[j];
 682       delete hZEMlg[j];
 683     }
 684   }
 685   delete hPMRefChg;
 686   delete hPMRefAhg;
 687   delete hPMRefClg;
 688   delete hPMRefAlg;
 689
 690   /* write report */
 691   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
 692
 693   /* close result file */
 694   fclose(fp);
 695
 696   /* report progress */
 697   daqDA_progressReport(90);
 698
 699   /* store the result file on FES */
 700   // [1] File with mapping
 701   status = daqDA_FES_storeFile(MAPDATA_FILE,MAPDATA_FILE);
 702   if(status){
 703     printf("Failed to export file : %d\n",status);
 704     return -1;
 705   }
 706   //
 707   // [2] File with laser data
 708   status = daqDA_FES_storeFile(LASDATA_FILE,LASDATA_FILE);
 709   if(status){
 710     printf("Failed to export file : %d\n",status);
 711     return -1;
 712   }
 713   // [3] File with laser histos
 714   status = daqDA_FES_storeFile(LASHISTO_FILE,LASHISTO_FILE);
 715   if(status){
 716     printf("Failed to export pedestal histos file to DAQ FES\n");
 717     return -1;
 718   }
 719
 720   /* report progress */
 721   daqDA_progressReport(100);
 722
 723   return status;
 724 }