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