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