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:
  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
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27 #include <Riostream.h>
  28
  29 // DATE
  30 #include <event.h>
  31 #include <monitor.h>
  32 #include <daqDA.h>
  33
  34 //ROOT
  35 #include <TRandom.h>
  36 #include <TH1F.h>
  37 #include <TF1.h>
  38 #include <TFile.h>
  39
  40 //AliRoot
  41 #include <AliRawReaderDate.h>
  42 #include <AliZDCRawStream.h>
  43
  44
  45 /* Main routine
  46       Arguments: list of DATE raw data files
  47 */
  48 int main(int argc, char **argv) {
  49
  50   int status = 0;
  51
  52   /* log start of process */
  53   printf("ZDC LASER program started\n");
  54
  55   /* check that we got some arguments = list of files */
  56   if (argc<2) {
  57     printf("Wrong number of arguments\n");
  58     return -1;
  59   }
  60
  61   // --- Histograms for LASER runs
  62   //     20 signal channels + 2 reference PTMs
  63   //
  64   TH1F::AddDirectory(0);
  65   TH1F *hZDCsideC[10], *hZDCsideA[10];
  66   char nhistZDCC[50], nhistZDCA[50];
  67   for(Int_t j=0; j<10; j++){
  68      if(j<5){ // ZNs
  69        sprintf(nhistZDCC,"ZNCtow%d",j);
  70        sprintf(nhistZDCA,"ZNAtow%d",j);
  71      }
  72      else if(j>=5 && j<10){ // ZPs
  73        sprintf(nhistZDCC,"ZPCtow%d",j);
  74        sprintf(nhistZDCA,"ZPAtow%d",j);
  75      }
  76      hZDCsideC[j] = new TH1F(nhistZDCC, nhistZDCC, 100, 0., 1000.);
  77      hZDCsideA[j] = new TH1F(nhistZDCA, nhistZDCA, 100, 0., 1000.);
  78   }
  79   TH1F *hPMRefsideC = new TH1F("hPMRefsideC","hPMRefsideC", 100,0.,1000.);
  80   TH1F *hPMRefsideA = new TH1F("hPMRefsideA","hPMRefsideA", 100,0.,1000.);
  81
  82
  83   /* open result file */
  84   FILE *fp=NULL;
  85   fp=fopen("./result.txt","a");
  86   if (fp==NULL) {
  87     printf("Failed to open file\n");
  88     return -1;
  89   }
  90
  91
  92   /* report progress */
  93   daqDA_progressReport(10);
  94
  95
  96   /* init some counters */
  97   int nevents_physics=0;
  98   int nevents_total=0;
  99
 100   /* read the data files */
 101   int n;
 102   for (n=1;n<argc;n++) {
 103
 104     status=monitorSetDataSource( argv[n] );
 105     if (status!=0) {
 106       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
 107       return -1;
 108     }
 109
 110     /* report progress */
 111     /* in this example, indexed on the number of files */
 112     daqDA_progressReport(10+80*n/argc);
 113
 114     /* read the file */
 115     for(;;) {
 116       struct eventHeaderStruct *event;
 117       eventTypeType eventT;
 118
 119       /* get next event */
 120       status=monitorGetEventDynamic((void **)&event);
 121       if (status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
 122       if (status!=0) {
 123         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 124         return -1;
 125       }
 126
 127       /* retry if got no event */
 128       if (event==NULL) {
 129         break;
 130       }
 131
 132
 133       /* use event - here, just write event id to result file */
 134       eventT=event->eventType;
 135
 136       if(eventT==PHYSICS_EVENT){
 137         //
 138         // *** To analyze EMD events you MUST have a pedestal data file!!!
 139         // *** -> check if a pedestal run has been analyzied
 140         FILE *filePed=NULL;
 141         filePed=fopen("./ZDCPedestal.dat","r");
 142         if (filePed==NULL) {
 143           printf("\t ERROR!!! You MUST have a ZDCPedestal.dat file!!!\n");
 144           return -1;
 145         }
 146         // 144 = 48 in-time + 48 out-of-time + 48 correlations
 147         Float_t readValues[2][144], MeanPed[44], MeanPedWidth[44],
 148               MeanPedOOT[44], MeanPedWidthOOT[44],
 149               CorrCoeff0[44], CorrCoeff1[44];
 150         //
 151         for(int jj=0; jj<144; jj++){
 152           for(int ii=0; ii<2; ii++){
 153              fscanf(filePed,"%f",&readValues[ii][jj]);
 154           }
 155           if(jj<48){
 156             MeanPed[jj] = readValues[0][jj];
 157             MeanPedWidth[jj] = readValues[1][jj];
 158           }
 159           else if(jj>48 && jj<96){
 160             MeanPedOOT[jj-48] = readValues[0][jj];
 161             MeanPedWidthOOT[jj-48] = readValues[1][jj];
 162           }
 163           else if(jj>144){
 164             CorrCoeff0[jj-96] = readValues[0][jj];
 165             CorrCoeff1[jj-96] = readValues[1][jj];;
 166           }
 167         }
 168         //
 169         // Initalize raw-data reading and decoding
 170         AliRawReader *reader = new AliRawReaderDate((void*)event);
 171         const AliRawDataHeader* header = reader->GetDataHeader();
 172         if(header) {
 173          UChar_t message = header->GetAttributes();
 174          if(message & 0x20){ // DEDICATED LASER RUN
 175             printf("\t STANDALONE_LASER_RUN raw data found\n");
 176             continue;
 177          }
 178          else{
 179             printf("\t NO STANDALONE_LASER_RUN raw data found\n");
 180             return -1;
 181          }
 182         }
 183         //Commented until we won't have true Raw Data Header...
 184         //else{
 185         //   printf("\t ATTENTION! No Raw Data Header found!!!\n");
 186           // return -1;
 187         //}
 188         //
 189         AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
 190         //
 191         if (!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 192         //
 193         while(rawStreamZDC->Next()){
 194           Int_t index=-1;
 195           // Implemented only for HIGH gain chain
 196           if((rawStreamZDC->IsADCDataWord()) && (rawStreamZDC->GetADCGain()==0)){
 197             index = rawStreamZDC->GetADCChannel();
 198             Float_t Pedestal = MeanPed[index];
 199             Float_t CorrADC = rawStreamZDC->GetADCValue() - Pedestal;
 200             if(rawStreamZDC->GetSector(0)==1){
 201               if(rawStreamZDC->GetSector(1)==5) hPMRefsideC->Fill(CorrADC);
 202               else hZDCsideC[rawStreamZDC->GetSector(1)]->Fill(CorrADC);
 203             }
 204             else if(rawStreamZDC->GetSector(0)==2){
 205               hZDCsideC[rawStreamZDC->GetSector(1)+5]->Fill(CorrADC);
 206             }
 207             else if(rawStreamZDC->GetSector(0)==4){
 208               if(rawStreamZDC->GetSector(1)==5) hPMRefsideA->Fill(CorrADC);
 209               else hZDCsideA[rawStreamZDC->GetSector(1)]->Fill(CorrADC);
 210             }
 211             else if(rawStreamZDC->GetSector(0)==5){
 212               hZDCsideA[rawStreamZDC->GetSector(1)+5]->Fill(CorrADC);
 213             }
 214           }//IsADCDataWord()
 215           //
 216          }
 217          //
 218          nevents_physics++;
 219          //
 220          delete reader;
 221          delete rawStreamZDC;
 222
 223       }//(if PHYSICS_EVENT)
 224       nevents_total++;
 225
 226       /* free resources */
 227       free(event);
 228
 229     }
 230   }
 231
 232   /* Analysis of the histograms */
 233   //
 234   Int_t maxBinC[10], maxBinA[10], maxBinRef[2];
 235   Int_t nBinC[10], nBinA[10], nBinRef[2];
 236   Float_t xMaxC[10], xMaxA[10], xMaxRef[2];
 237   Float_t maxXvalC[10], maxXvalA[10], maxXvalRef[2];
 238   Float_t xlowC[10], xlowA[10], xlowRef[10];
 239   TF1 *funA[10], *funC[10], *funRef[2];
 240   //
 241   Float_t meanC[10], meanA[10], meanRef[2];
 242   Float_t sigmaA[10], sigmaC[10], sigmaRef[10];
 243   //
 244   for(Int_t k=0; k<10; k++){
 245      maxBinC[k] = hZDCsideC[k]->GetMaximumBin();
 246      nBinC[k] = (hZDCsideC[k]->GetXaxis())->GetNbins();
 247      xMaxC[k] = (hZDCsideC[k]->GetXaxis())->GetXmax();
 248      maxXvalC[k] = maxBinC[k]*xMaxC[k]/nBinC[k];
 249      //
 250      if(maxXvalC[k]-100.<0.) {xlowC[k]=0.;}
 251      else xlowC[k] = maxXvalC[k];
 252      hZDCsideC[k]->Fit("gaus","Q","",xlowC[k],maxXvalC[k]+100.);
 253      funC[k] = hZDCsideC[k]->GetFunction("gaus");
 254      meanC[k] = (Float_t) (funC[k]->GetParameter(1));
 255      sigmaC[k] = (Float_t) (funC[k]->GetParameter(2));
 256      //
 257      maxBinA[k] = hZDCsideA[k]->GetMaximumBin();
 258      nBinA[k] = (hZDCsideA[k]->GetXaxis())->GetNbins();
 259      xMaxA[k] = (hZDCsideA[k]->GetXaxis())->GetXmax();
 260      maxXvalA[k] = maxBinA[k]*xMaxA[k]/nBinA[k];
 261      //
 262      if(maxXvalA[k]-100.<0.) {xlowA[k]=0.;}
 263      else xlowA[k] = maxXvalA[k];
 264      hZDCsideA[k]->Fit("gaus","Q","",xlowA[k],maxXvalA[k]+100.);
 265      funA[k] = hZDCsideC[k]->GetFunction("gaus");
 266      meanA[k] = (Float_t) (funA[k]->GetParameter(1));
 267      sigmaA[k] = (Float_t) (funA[k]->GetParameter(2));
 268   }
 269   //
 270   maxBinRef[0] = hPMRefsideC->GetMaximumBin();
 271   nBinRef[0] = (hPMRefsideC->GetXaxis())->GetNbins();
 272   xMaxRef[0] = (hPMRefsideC->GetXaxis())->GetXmax();
 273   maxXvalRef[0] = maxBinRef[0]*xMaxRef[0]/nBinRef[0];
 274   //
 275   if(maxXvalRef[0]-100.<0.) {xlowRef[0]=0.;}
 276   else xlowRef[0] = maxXvalRef[0];
 277   hPMRefsideC->Fit("gaus","Q","",xlowRef[0],maxXvalRef[0]+100.);
 278   funRef[0] = hPMRefsideC->GetFunction("gaus");
 279   meanRef[0] = (Float_t) (funRef[0]->GetParameter(1));
 280   sigmaRef[0] = (Float_t) (funRef[0]->GetParameter(2));
 281   //
 282   maxBinRef[1] = hPMRefsideA->GetMaximumBin();
 283   nBinRef[1] = (hPMRefsideA->GetXaxis())->GetNbins();
 284   xMaxRef[1] = (hPMRefsideA->GetXaxis())->GetXmax();
 285   maxXvalRef[1] = maxBinRef[1]*xMaxRef[1]/nBinRef[1];
 286   //
 287   if(maxXvalRef[1]-100.<0.) {xlowRef[1]=0.;}
 288   else xlowRef[1] = maxXvalRef[1];
 289   hPMRefsideA->Fit("gaus","Q","",xlowRef[1],maxXvalRef[1]+100.);
 290   funRef[1] = hPMRefsideA->GetFunction("gaus");
 291   meanRef[1] = (Float_t) (funRef[1]->GetParameter(1));
 292   sigmaRef[1] = (Float_t) (funRef[1]->GetParameter(2));
 293   //
 294   FILE *fileShuttle;
 295   const char *fName = "ZDCLaser.dat";
 296   fileShuttle = fopen(fName,"w");
 297   for(Int_t i=0; i<10; i++) fprintf(fileShuttle,"\t%f\t%f\n",meanC[i], sigmaC[i]);
 298   for(Int_t i=0; i<10; i++) fprintf(fileShuttle,"\t%f\t%f\n",meanA[i], sigmaA[i]);
 299   for(Int_t i=0; i<2; i++)  fprintf(fileShuttle,"\t%f\t%f\n",meanRef[i], sigmaRef[i]);
 300   //
 301   fclose(fileShuttle);
 302   //
 303   for(Int_t j=0; j<10; j++){
 304      delete hZDCsideC[j];
 305      delete hZDCsideA[j];
 306      delete hPMRefsideC;
 307      delete hPMRefsideA;
 308   }
 309
 310   /* write report */
 311   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
 312
 313   /* close result file */
 314   fclose(fp);
 315
 316   /* report progress */
 317   daqDA_progressReport(90);
 318
 319   /* store the result file on FES */
 320   status = daqDA_FES_storeFile(fName,"ZDCLASER_data");
 321   if(status){
 322     printf("Failed to export file : %d\n",status);
 323     return -1;
 324   }
 325
 326   /* report progress */
 327   daqDA_progressReport(100);
 328
 329
 330   return status;
 331 }