ZDC/ZDCPEDESTALda.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_PEDESTAL_RUN
  17 DA Type: LDC
  18 Number of events needed: no constraint (tipically ~10^3)
  19 Input Files:
  20 Output Files: ZDCPedestal.dat, ZDCChMapping.dat
  21 Trigger Types Used: Standalone Trigger
  22
  23 */
  24 #define PEDDATA_FILE  "ZDCPedestal.dat"
  25 #define MAPDATA_FILE  "ZDCChMapping.dat"
  26
  27 #include <stdio.h>
  28 #include <stdlib.h>
  29 #include <Riostream.h>
  30
  31 // DATE
  32 #include <event.h>
  33 #include <monitor.h>
  34 #include <daqDA.h>
  35
  36 //ROOT
  37 #include <TRandom.h>
  38 #include <TH1F.h>
  39 #include <TH2F.h>
  40 #include <TProfile.h>
  41 #include <TF1.h>
  42 #include <TFile.h>
  43 #include <TFitter.h>
  44
  45 //AliRoot
  46 #include <AliRawReaderDate.h>
  47 #include <AliRawEventHeaderBase.h>
  48 #include <AliZDCRawStream.h>
  49
  50
  51 /* Main routine
  52       Arguments: list of DATE raw data files
  53 */
  54 int main(int argc, char **argv) {
  55
  56   TFitter *minuitFit = new TFitter(4);
  57   TVirtualFitter::SetFitter(minuitFit);
  58
  59   int status = 0;
  60
  61   /* log start of process */
  62   printf("\n ZDC PEDESTAL program started\n");
  63
  64   /* check that we got some arguments = list of files */
  65   if (argc<2) {
  66     printf("Wrong number of arguments\n");
  67     return -1;
  68   }
  69
  70   // --- Histograms for ADC pedestals
  71   //     [22 signal channels + 2 reference PTMs]  x 2 gain chains
  72   //
  73   TH1F::AddDirectory(0);
  74   int const kNChannels = 24;
  75   TH1F *hPedhg[kNChannels], *hPedOutOfTimehg[kNChannels];
  76   TH2F *hPedCorrhg[kNChannels];
  77   TH1F *hPedlg[kNChannels], *hPedOutOfTimelg[kNChannels];
  78   TH2F *hPedCorrlg[kNChannels];
  79   //
  80   char namhist1hg[50], namhist2hg[50], namhist3hg[50];
  81   char namhist1lg[50], namhist2lg[50], namhist3lg[50];
  82   for(Int_t j=0; j<kNChannels; j++){
  83      if(j<=4){ // ZNC
  84        sprintf(namhist1hg,"PedZNChg_%d",j);
  85        sprintf(namhist2hg,"PedZNChgOutOfTime_%d",j);
  86        sprintf(namhist3hg,"PedCorrZNChg_%d",j);
  87        //
  88        sprintf(namhist1lg,"PedZNClg_%d",j);
  89        sprintf(namhist2lg,"PedZNClgOutOfTime_%d",j);
  90        sprintf(namhist3lg,"PedCorrZNClg_%d",j);
  91      }
  92      else if(j>=5 && j<=9){ // ZPC
  93        sprintf(namhist1hg,"PedZPChg_%d",j-5);
  94        sprintf(namhist2hg,"PedZPChgOutOfTime_%d",j-5);
  95        sprintf(namhist3hg,"PedCorrZPChg_%d",j-5);
  96        //
  97        sprintf(namhist1lg,"PedZPClg_%d",j-5);
  98        sprintf(namhist2lg,"PedZPClgOutOfTime_%d",j-5);
  99        sprintf(namhist3lg,"PedCorrZPClg_%d",j-5);
 100      }
 101      else if(j==10 || j==11){ // ZEM
 102        sprintf(namhist1hg,"PedZEMhg_%d",j-10);
 103        sprintf(namhist2hg,"PedZEMhgOutOfTime_%d",j-10);
 104        sprintf(namhist3hg,"PedCorrZEMhg_%d",j-10);
 105        //
 106        sprintf(namhist1lg,"PedZEMlg_%d",j-10);
 107        sprintf(namhist2lg,"PedZEMlgOutOfTime_%d",j-10);
 108        sprintf(namhist3lg,"PedCorrZEMlg_%d",j-10);
 109      }
 110      else if(j>=12 && j<=16){ // ZNA
 111        sprintf(namhist1hg,"PedZNAhg_%d",j-12);
 112        sprintf(namhist2hg,"PedZNAhgOutOfTime_%d",j-12);
 113        sprintf(namhist3hg,"PedCorrZNAhg_%d",j-12);
 114        //
 115        sprintf(namhist1lg,"PedZNAlg_%d",j-12);
 116        sprintf(namhist2lg,"PedZNAlgOutOfTime_%d",j-12);
 117        sprintf(namhist3lg,"PedCorrZNAlg_%d",j-12);
 118      }
 119      else if(j>=17 && j<=21){ // ZPA
 120        sprintf(namhist1hg,"PedZPAhg_%d",j-17);
 121        sprintf(namhist2hg,"PedZPAhgOutOfTime_%d",j-17);
 122        sprintf(namhist3hg,"PedCorrZPAhg_%d",j-17);
 123        //
 124        sprintf(namhist1lg,"PedZPAlg_%d",j-17);
 125        sprintf(namhist2lg,"PedZPAlgOutOfTime_%d",j-17);
 126        sprintf(namhist3lg,"PedCorrZPAlg_%d",j-17);
 127      }
 128      else if(j==22 || j==24){ //Reference PMs
 129        sprintf(namhist1hg,"PedRefhg_%d",j-22);
 130        sprintf(namhist2hg,"PedRefhgOutOfTime_%d",j-22);
 131        sprintf(namhist3hg,"PedCorrRefhg_%d",j-22);
 132        //
 133        sprintf(namhist1lg,"PedReflg_%d",j-22);
 134        sprintf(namhist2lg,"PedReflgOutOfTime_%d",j-22);
 135        sprintf(namhist3lg,"PedCorrReflg_%d",j-22);
 136      }
 137      // --- High gain chain histos
 138      hPedhg[j] = new TH1F(namhist1hg, namhist1hg, 200,0., 200.);
 139      hPedOutOfTimehg[j] = new TH1F(namhist2hg, namhist2hg, 200,0., 200.);
 140      hPedCorrhg[j] = new TH2F(namhist3hg,namhist3hg,100,0.,200.,100,0.,200.);
 141      // --- Low gain chain histos
 142      hPedlg[j] = new TH1F(namhist1lg, namhist1lg, 100,0., 600.);
 143      hPedOutOfTimelg[j] = new TH1F(namhist2lg, namhist2lg, 100,0., 600.);
 144      hPedCorrlg[j] = new TH2F(namhist3lg,namhist3lg,100,0.,600.,100,0.,600.);
 145   }
 146
 147
 148   /* open result file */
 149   FILE *fp=NULL;
 150   fp=fopen("./result.txt","a");
 151   if (fp==NULL) {
 152     printf("Failed to open file\n");
 153     return -1;
 154   }
 155
 156   FILE *mapFile4Shuttle;
 157
 158   /* report progress */
 159   daqDA_progressReport(10);
 160
 161
 162   /* init some counters */
 163   int nevents_physics=0;
 164   int nevents_total=0;
 165
 166   /* read the data files */
 167   int n;
 168   for(n=1;n<argc;n++){
 169
 170     status=monitorSetDataSource( argv[n] );
 171     if (status!=0) {
 172       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
 173       return -1;
 174     }
 175
 176     /* report progress */
 177     /* in this example, indexed on the number of files */
 178     daqDA_progressReport(10+80*n/argc);
 179
 180     /* read the file */
 181     for(;;) {
 182       struct eventHeaderStruct *event;
 183       eventTypeType eventT;
 184
 185       /* get next event */
 186       status=monitorGetEventDynamic((void **)&event);
 187       if(status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
 188       if(status!=0) {
 189         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 190         return -1;
 191       }
 192
 193       /* retry if got no event */
 194       if(event==NULL) {
 195         break;
 196       }
 197
 198       // Initalize raw-data reading and decoding
 199       AliRawReader *reader = new AliRawReaderDate((void*)event);
 200       reader->Select("ZDC");
 201       // --- Reading event header
 202       //UInt_t evtype = reader->GetType();
 203       //printf("\n\t ZDCPEDESTALda -> ev. type %d\n",evtype);
 204       //printf("\t ZDCPEDESTALda -> run # %d\n",reader->GetRunNumber());
 205       //
 206       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
 207
 208
 209       /* use event - here, just write event id to result file */
 210       eventT=event->eventType;
 211
 212       Int_t ich=0, adcMod[48], adcCh[48], sigCode[48], det[48], sec[48];
 213
 214       if(eventT==START_OF_DATA){
 215
 216         rawStreamZDC->SetSODReading(kTRUE);
 217
 218         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 219         else{
 220           while(rawStreamZDC->Next()){
 221             if(rawStreamZDC->IsChMapping()){
 222               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
 223               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
 224               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
 225               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
 226               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
 227               ich++;
 228             }
 229           }
 230         }
 231         // --------------------------------------------------------
 232         // --- Writing ascii data file for the Shuttle preprocessor
 233         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
 234         for(Int_t i=0; i<ich; i++){
 235            fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",i,
 236              adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
 237            //
 238            //printf("ZDCPEDESTALDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
 239            //      i,adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
 240         }
 241         fclose(mapFile4Shuttle);
 242       }
 243
 244       if(eventT==PHYSICS_EVENT){
 245         fprintf(fp,"Run #%lu, event size: %lu, BC:%u, Orbit:%u, Period:%u\n",
 246           (unsigned long)event->eventRunNb,
 247           (unsigned long)event->eventSize,
 248           EVENT_ID_GET_BUNCH_CROSSING(event->eventId),
 249           EVENT_ID_GET_ORBIT(event->eventId),
 250           EVENT_ID_GET_PERIOD(event->eventId));
 251
 252         // --- Reading data header
 253         reader->ReadHeader();
 254         const AliRawDataHeader* header = reader->GetDataHeader();
 255         if(header){
 256          UChar_t message = header->GetAttributes();
 257          if(message & 0x20){ // PEDESTAL RUN
 258             //printf("\t STANDALONE_PEDESTAL RUN raw data found\n");
 259          }
 260          else{
 261             printf("\t NO STANDALONE_PEDESTAL RUN raw data found\n");
 262             return -1;
 263          }
 264         }
 265         else{
 266            printf("\t ATTENTION! No Raw Data Header found!!!\n");
 267            return -1;
 268         }
 269
 270         rawStreamZDC->SetSODReading(kTRUE);
 271
 272         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 273         //
 274         // ----- Setting ch. mapping -----
 275         for(Int_t jk=0; jk<48; jk++){
 276           //printf("ZDCPEDESTALDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
 277           //    jk,adcMod[jk],adcCh[jk],sigCode[jk],det[jk],sec[jk]);
 278           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
 279           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
 280           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
 281           rawStreamZDC->SetMapDet(jk, det[jk]);
 282           rawStreamZDC->SetMapTow(jk, sec[jk]);
 283         }
 284         //
 285         Int_t iraw=0;
 286         Int_t RawADChg[kNChannels], RawADCoothg[kNChannels];
 287         Int_t RawADClg[kNChannels], RawADCootlg[kNChannels];
 288         for(Int_t j=0; j<kNChannels; j++){
 289            RawADChg[j]=0; RawADCoothg[j]=0;
 290            RawADClg[j]=0; RawADCootlg[j]=0;
 291         }
 292         //
 293         while(rawStreamZDC->Next()){
 294          Int_t index=-1;
 295          Int_t detector = rawStreamZDC->GetSector(0);
 296          if(rawStreamZDC->IsADCDataWord() && (detector!=-1)){
 297           if(rawStreamZDC->GetSector(1)!=5){ // Physics signals
 298             if(detector==1) index = rawStreamZDC->GetSector(1); // *** ZNC
 299             else if(detector==2) index = rawStreamZDC->GetSector(1)+5; // *** ZPC
 300             else if(detector==3) index = rawStreamZDC->GetSector(1)+9; // *** ZEM
 301             else if(detector==4) index = rawStreamZDC->GetSector(1)+12; // *** ZNA
 302             else if(detector==5) index = rawStreamZDC->GetSector(1)+17; // *** ZPA
 303           }
 304           else{ // Reference PMs
 305             index = (detector-1)/3+22;
 306           }
 307           //
 308           if(index==-1) printf("\tERROR!!! iraw %d det %d quad %d res %d index %d ADC %d\n",
 309             iraw, detector, rawStreamZDC->GetSector(1),
 310             rawStreamZDC->GetADCGain(), index, rawStreamZDC->GetADCValue());
 311
 312            //
 313            if(iraw<2*kNChannels){ // --- In-time pedestals (1st 48 raw data)
 314             if(rawStreamZDC->GetADCGain()==0){
 315               hPedhg[index]->Fill(rawStreamZDC->GetADCValue());
 316               RawADChg[index] = rawStreamZDC->GetADCValue();
 317               //
 318               //printf("\t filling histo hPedhg[%d]\n",index);
 319             }
 320             else{
 321               hPedlg[index]->Fill(rawStreamZDC->GetADCValue());
 322               RawADClg[index] = rawStreamZDC->GetADCValue();
 323               //
 324               //printf("\t filling histo hPedlg[%d]\n",index);
 325             }
 326            }
 327            else{  // --- Out-of-time pedestals
 328             if(rawStreamZDC->GetADCGain()==0){
 329               hPedOutOfTimehg[index]->Fill(rawStreamZDC->GetADCValue());
 330               RawADCoothg[index] = rawStreamZDC->GetADCValue();
 331               //
 332               //printf("\t filling histo hPedOutOfTimehg[%d]\n",index);
 333             }
 334             else{
 335               hPedOutOfTimelg[index]->Fill(rawStreamZDC->GetADCValue());
 336               RawADCootlg[index] = rawStreamZDC->GetADCValue();
 337               //
 338               //printf("\t filling histo hPedOutOfTimelg[%d]\n",index);
 339             }
 340            }
 341             iraw++;
 342           }//IsADCDataWord()
 343           //
 344           if(iraw == 4*kNChannels){ // Last ADC channel -> Filling correlation histos
 345             for(Int_t k=0; k<kNChannels; k++){
 346               hPedCorrhg[k]->Fill(RawADCoothg[k], RawADChg[k]);
 347               hPedCorrlg[k]->Fill(RawADCootlg[k], RawADClg[k]);
 348             }
 349           }
 350          }
 351          //
 352          nevents_physics++;
 353          //
 354          delete reader;
 355          delete rawStreamZDC;
 356
 357       }//(if PHYSICS_EVENT)
 358       nevents_total++;
 359
 360       /* free resources */
 361       free(event);
 362
 363     }
 364   }
 365
 366   /* Analysis of the histograms */
 367   //
 368   FILE *fileShuttle;
 369   fileShuttle = fopen(PEDDATA_FILE,"w");
 370   //
 371   Float_t MeanPed[2*kNChannels], MeanPedWidth[2*kNChannels],
 372         MeanPedOOT[2*kNChannels], MeanPedWidthOOT[2*kNChannels];
 373   // --- In-time pedestals
 374   TF1 *ADCfunchg[kNChannels];
 375   for(Int_t i=0; i<kNChannels; i++){
 376      hPedhg[i]->Fit("gaus","Q");
 377      ADCfunchg[i] = hPedhg[i]->GetFunction("gaus");
 378      MeanPed[i] = (Double_t) ADCfunchg[i]->GetParameter(1);
 379      MeanPedWidth[i] = (Double_t)  ADCfunchg[i]->GetParameter(2);
 380      fprintf(fileShuttle,"\t%f\t%f\n",MeanPed[i],MeanPedWidth[i]);
 381      printf("\t MeanPedhg[%d] = %f\n",i, MeanPed[i]);
 382   }
 383   TF1 *ADCfunclg[kNChannels];
 384   for(Int_t i=0; i<kNChannels; i++){
 385      hPedlg[i]->Fit("gaus","Q");
 386      ADCfunclg[i] = hPedlg[i]->GetFunction("gaus");
 387      MeanPed[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(1);
 388      MeanPedWidth[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(2);
 389      fprintf(fileShuttle,"\t%f\t%f\n",MeanPed[i+kNChannels],MeanPedWidth[i+kNChannels]);
 390      printf("\t MeanPedlg[%d] = %f\n",i+kNChannels, MeanPed[i+kNChannels]);
 391   }
 392   // --- Out-of-time pedestals
 393   TF1 *ADCootfunchg[kNChannels];
 394   for(Int_t i=0; i<kNChannels; i++){
 395      hPedOutOfTimehg[i]->Fit("gaus","Q");
 396      ADCootfunchg[i] = hPedOutOfTimehg[i]->GetFunction("gaus");
 397      MeanPedOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(1);
 398      MeanPedWidthOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(2);
 399      fprintf(fileShuttle,"\t%f\t%f\n",MeanPedOOT[i],MeanPedWidthOOT[i]);
 400      printf("\t MeanPedOOThg[%d] = %f\n",i, MeanPedOOT[i]);
 401   }
 402   TF1 *ADCootfunclg[kNChannels];
 403   for(Int_t i=0; i<kNChannels; i++){
 404      hPedOutOfTimelg[i]->Fit("gaus","Q");
 405      ADCootfunclg[i] = hPedOutOfTimelg[i]->GetFunction("gaus");
 406      MeanPedOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(1);
 407      MeanPedWidthOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(2);
 408      fprintf(fileShuttle,"\t%f\t%f\n",MeanPedOOT[i+kNChannels],MeanPedWidthOOT[i+kNChannels]);
 409      printf("\t MeanPedOOTlg[%d] = %f\n",i+kNChannels, MeanPedOOT[i+kNChannels]);
 410   }
 411
 412   // ***************************************************
 413   //   Unless we have a narrow correlation to fit we
 414   //    don't fit and store in-time vs. out-of-time
 415   //    histograms -> mean pedstal subtracted!!!!!!
 416   // ***************************************************
 417   // --- Correlations
 418 /*
 419   Float_t CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
 420   TProfile *hPedCorrProfhg[kNChannels], *hPedCorrProflg[kNChannels];
 421   TF1 *ffunchg[kNChannels], *ffunclg[kNChannels];
 422   char namhist4[50];
 423   for(int i=0;i<kNChannels;i++) {
 424      sprintf(namhist4,"ADCHRvsOOT%d_Prof",i);
 425      hPedCorrProfhg[i] = hPedCorrhg[i]->ProfileX(namhist4,-1,-1,"S");
 426      hPedCorrProfhg[i]->SetName(namhist4);
 427      hPedCorrProfhg[i]->Fit("pol1","Q");
 428      ffunchg[i] = hPedCorrProfhg[i]->GetFunction("pol1");
 429      CorrCoeff0[i] = (Double_t)  ffunchg[i]->GetParameter(0);
 430      CorrCoeff1[i] = (Double_t) ffunchg[i]->GetParameter(1);
 431      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,CorrCoeff0[i],CorrCoeff1[i]);
 432      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",i, CorrCoeff0[i], i, CorrCoeff1[i]);
 433   }
 434   for(int i=0;i<kNChannels;i++) {
 435      sprintf(namhist4,"ADCLRvsOOT%d_Prof",i);
 436      hPedCorrProflg[i] = hPedCorrlg[i]->ProfileX(namhist4,-1,-1,"S");
 437      hPedCorrProflg[i]->SetName(namhist4);
 438      hPedCorrProflg[i]->Fit("pol1","Q");
 439      ffunclg[i] = hPedCorrProflg[i]->GetFunction("pol1");
 440      CorrCoeff0[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(0);
 441      CorrCoeff1[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(1);
 442      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,CorrCoeff0[i+kNChannels],CorrCoeff1[i+kNChannels]);
 443      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",
 444      //         i+kNChannels, CorrCoeff0[i+kNChannels], i+kNChannels, CorrCoeff1[i+kNChannels]);
 445   }
 446 */
 447   //
 448   fclose(fileShuttle);
 449   //
 450   for(Int_t j=0; j<kNChannels; j++){
 451      delete hPedhg[j];
 452      delete hPedOutOfTimehg[j];
 453      delete hPedCorrhg[j];
 454      delete hPedlg[j];
 455      delete hPedOutOfTimelg[j];
 456      delete hPedCorrlg[j];
 457   }
 458
 459   //delete minuitFit;
 460   TVirtualFitter::SetFitter(0);
 461
 462   /* write report */
 463   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
 464
 465   /* close result file */
 466   fclose(fp);
 467
 468   /* report progress */
 469   daqDA_progressReport(90);
 470
 471   /* store the result files on FES */
 472   // [1] File with mapping
 473   status = daqDA_FES_storeFile(MAPDATA_FILE,MAPDATA_FILE);
 474   if(status){
 475     printf("Failed to export file %d to DAQ FES\n",status);
 476     return -1;
 477   }
 478   // [2] File with pedestal data
 479   status = daqDA_FES_storeFile(PEDDATA_FILE,PEDDATA_FILE);
 480   if(status){
 481     printf("Failed to export file %d to DAQ FES\n",status);
 482     return -1;
 483   }
 484
 485   /* store the result files on DB */
 486   status = daqDA_DB_storeFile(PEDDATA_FILE,PEDDATA_FILE);
 487   if(status){
 488     printf("Failed to export file %d to DAQ DB\n",status);
 489     return -1;
 490   }
 491
 492   /* report progress */
 493   daqDA_progressReport(100);
 494
 495   return status;
 496 }