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
  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 <TH2F.h>
  38 #include <TProfile.h>
  39 #include <TF1.h>
  40 #include <TFile.h>
  41 #include <TFitter.h>
  42
  43 //AliRoot
  44 #include <AliRawReaderDate.h>
  45 #include <AliRawEventHeaderBase.h>
  46 #include <AliZDCRawStream.h>
  47
  48
  49 /* Main routine
  50       Arguments: list of DATE raw data files
  51 */
  52 int main(int argc, char **argv) {
  53
  54   TFitter *minuitFit = new TFitter(4);
  55   TVirtualFitter::SetFitter(minuitFit);
  56
  57   int status = 0;
  58
  59   /* log start of process */
  60   printf("ZDC PEDESTAL program started\n");
  61
  62   /* check that we got some arguments = list of files */
  63   if (argc<2) {
  64     printf("Wrong number of arguments\n");
  65     return -1;
  66   }
  67
  68   // --- Histograms for ADC pedestals
  69   //     [22 signal channels + 2 reference PTMs]  x 2 gain chains
  70   //
  71   TH1F::AddDirectory(0);
  72   int const kNChannels = 24;
  73   TH1F *hPedhg[kNChannels], *hPedOutOfTimehg[kNChannels];
  74   TH2F *hPedCorrhg[kNChannels];
  75   TH1F *hPedlg[kNChannels], *hPedOutOfTimelg[kNChannels];
  76   TH2F *hPedCorrlg[kNChannels];
  77   //
  78   char namhist1hg[50], namhist2hg[50], namhist3hg[50];
  79   char namhist1lg[50], namhist2lg[50], namhist3lg[50];
  80   for(Int_t j=0; j<kNChannels; j++){
  81      if(j<5){ // ZN1
  82        sprintf(namhist1hg,"PedZN1hg_%d",j);
  83        sprintf(namhist2hg,"PedZN1hgOutOfTime_%d",j);
  84        sprintf(namhist3hg,"PedCorrZN1hg_%d",j);
  85        //
  86        sprintf(namhist1lg,"PedZN1lg_%d",j);
  87        sprintf(namhist2lg,"PedZN1lgOutOfTime_%d",j);
  88        sprintf(namhist3lg,"PedCorrZN1lg_%d",j);
  89      }
  90      else if(j>=5 && j<10){ // ZP1
  91        sprintf(namhist1hg,"PedZP1hg_%d",j-5);
  92        sprintf(namhist2hg,"PedZP1hgOutOfTime_%d",j-5);
  93        sprintf(namhist3hg,"PedCorrZP1hg_%d",j-5);
  94        //
  95        sprintf(namhist1lg,"PedZP1lg_%d",j-5);
  96        sprintf(namhist2lg,"PedZP1lgOutOfTime_%d",j-5);
  97        sprintf(namhist3lg,"PedCorrZP1lg_%d",j-5);
  98      }
  99      else if(j>=10 && j<12){ // ZEM
 100        sprintf(namhist1hg,"PedZEMhg_%d",j-10);
 101        sprintf(namhist2hg,"PedZEMhgOutOfTime_%d",j-10);
 102        sprintf(namhist3hg,"PedCorrZEMhg_%d",j-10);
 103        //
 104        sprintf(namhist1lg,"PedZEMlg_%d",j-10);
 105        sprintf(namhist2lg,"PedZEMlgOutOfTime_%d",j-10);
 106        sprintf(namhist3lg,"PedCorrZEMlg_%d",j-10);
 107      }
 108      else if(j>=12 && j<17){ // ZN2
 109        sprintf(namhist1hg,"PedZN2hg_%d",j-12);
 110        sprintf(namhist2hg,"PedZN2hgOutOfTime_%d",j-12);
 111        sprintf(namhist3hg,"PedCorrZN2hg_%d",j-12);
 112        //
 113        sprintf(namhist1lg,"PedZN2lg_%d",j-12);
 114        sprintf(namhist2lg,"PedZN2lgOutOfTime_%d",j-12);
 115        sprintf(namhist3lg,"PedCorrZN2lg_%d",j-12);
 116      }
 117      else if(j>=17 && j<22){ // ZP2
 118        sprintf(namhist1hg,"PedZP2hg_%d",j-17);
 119        sprintf(namhist2hg,"PedZP2hgOutOfTime_%d",j-17);
 120        sprintf(namhist3hg,"PedCorrZP2hg_%d",j-17);
 121        //
 122        sprintf(namhist1lg,"PedZP2lg_%d",j-17);
 123        sprintf(namhist2lg,"PedZP2lgOutOfTime_%d",j-17);
 124        sprintf(namhist3lg,"PedCorrZP2lg_%d",j-17);
 125      }
 126      else if(j>=22 && j<24){ //Reference PMs
 127        sprintf(namhist1hg,"PedRefhg_%d",j-22);
 128        sprintf(namhist2hg,"PedRefhgOutOfTime_%d",j-22);
 129        sprintf(namhist3hg,"PedCorrRefhg_%d",j-22);
 130        //
 131        sprintf(namhist1lg,"PedReflg_%d",j-22);
 132        sprintf(namhist2lg,"PedReflgOutOfTime_%d",j-22);
 133        sprintf(namhist3lg,"PedCorrReflg_%d",j-22);
 134      }
 135      // --- High gain chain histos
 136      hPedhg[j] = new TH1F(namhist1hg, namhist1hg, 200,0., 200.);
 137      hPedOutOfTimehg[j] = new TH1F(namhist2hg, namhist2hg, 200,0., 200.);
 138      hPedCorrhg[j] = new TH2F(namhist3hg,namhist3hg,100,0.,200.,100,0.,200.);
 139      // --- Low gain chain histos
 140      hPedlg[j] = new TH1F(namhist1lg, namhist1lg, 100,0., 600.);
 141      hPedOutOfTimelg[j] = new TH1F(namhist2lg, namhist2lg, 100,0., 600.);
 142      hPedCorrlg[j] = new TH2F(namhist3lg,namhist3lg,100,0.,600.,100,0.,600.);
 143   }
 144
 145
 146   /* open result file */
 147   FILE *fp=NULL;
 148   fp=fopen("./result.txt","a");
 149   if (fp==NULL) {
 150     printf("Failed to open file\n");
 151     return -1;
 152   }
 153
 154   FILE *mapFile4Shuttle;
 155   const char *mapfName = "ZDCChMapping.dat";
 156
 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       if(eventT==START_OF_DATA){
 214
 215         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 216         else{
 217           while(rawStreamZDC->Next()){
 218             if(rawStreamZDC->IsChMapping()){
 219               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
 220               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
 221               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
 222               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
 223               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
 224               ich++;
 225             }
 226           }
 227         }
 228         // --------------------------------------------------------
 229         // --- Writing ascii data file for the Shuttle preprocessor
 230         mapFile4Shuttle = fopen(mapfName,"w");
 231         for(Int_t i=0; i<ich; i++){
 232            fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",i,
 233              adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
 234            //
 235            //printf("ZDCPEDESTALDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
 236            //      i,adcMod[i],adcCh[i],sigCode[i],det[i],sec[i]);
 237         }
 238         fclose(mapFile4Shuttle);
 239       }
 240
 241       if(eventT==PHYSICS_EVENT){
 242         fprintf(fp,"Run #%lu, event size: %lu, BC:%u, Orbit:%u, Period:%u\n",
 243           (unsigned long)event->eventRunNb,
 244           (unsigned long)event->eventSize,
 245           EVENT_ID_GET_BUNCH_CROSSING(event->eventId),
 246           EVENT_ID_GET_ORBIT(event->eventId),
 247           EVENT_ID_GET_PERIOD(event->eventId));
 248
 249         // --- Reading data header
 250         reader->ReadHeader();
 251         const AliRawDataHeader* header = reader->GetDataHeader();
 252         if(header){
 253          UChar_t message = header->GetAttributes();
 254          if(message & 0x20){ // PEDESTAL RUN
 255             //printf("\t STANDALONE_PEDESTAL RUN raw data found\n");
 256          }
 257          else{
 258             printf("\t NO STANDALONE_PEDESTAL RUN raw data found\n");
 259             return -1;
 260          }
 261         }
 262         else{
 263            printf("\t ATTENTION! No Raw Data Header found!!!\n");
 264            return -1;
 265         }
 266
 267         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 268         //
 269         // ----- Setting ch. mapping -----
 270         for(Int_t jk=0; jk<48; jk++){
 271           rawStreamZDC->SetMapADCMod(jk, adcMod[jk]);
 272           rawStreamZDC->SetMapADCCh(jk, adcCh[jk]);
 273           rawStreamZDC->SetMapADCSig(jk, sigCode[jk]);
 274           rawStreamZDC->SetMapDet(jk, det[jk]);
 275           rawStreamZDC->SetMapTow(jk, sec[jk]);
 276         }
 277         //
 278         Int_t iraw=0;
 279         Int_t RawADChg[kNChannels], RawADCoothg[kNChannels];
 280         Int_t RawADClg[kNChannels], RawADCootlg[kNChannels];
 281         for(Int_t j=0; j<kNChannels; j++){
 282            RawADChg[j]=0; RawADCoothg[j]=0;
 283            RawADClg[j]=0; RawADCootlg[j]=0;
 284         }
 285         //
 286         while(rawStreamZDC->Next()){
 287          Int_t index=-1;
 288          if(rawStreamZDC->IsADCDataWord()){
 289           if(rawStreamZDC->GetSector(1)!=5){ // Physics signals
 290             if(rawStreamZDC->GetSector(0)==1) index = rawStreamZDC->GetSector(1); // *** ZNC
 291             else if(rawStreamZDC->GetSector(0)==2) index = rawStreamZDC->GetSector(1)+5; // *** ZPC
 292             else if(rawStreamZDC->GetSector(0)==3) index = rawStreamZDC->GetSector(1)+9; // *** ZEM
 293             else if(rawStreamZDC->GetSector(0)==4) index = rawStreamZDC->GetSector(1)+12; // *** ZNA
 294             else if(rawStreamZDC->GetSector(0)==5) index = rawStreamZDC->GetSector(1)+17; // *** ZPA
 295           }
 296           else{ // Reference PMs
 297             index = (rawStreamZDC->GetSector(0)-1)/3+22;
 298           }
 299           //
 300           /*printf("\t iraw %d index %d det %d quad %d res %d ADC %d\n", iraw, index,
 301                 rawStreamZDC->GetSector(0), rawStreamZDC->GetSector(1),
 302                 rawStreamZDC->GetADCGain(), rawStreamZDC->GetADCValue());
 303           */
 304            //
 305            if(iraw<2*kNChannels){ // --- In-time pedestals (1st 48 raw data)
 306             if(rawStreamZDC->GetADCGain()==0){
 307               hPedhg[index]->Fill(rawStreamZDC->GetADCValue());
 308               RawADChg[index] = rawStreamZDC->GetADCValue();
 309               //
 310               //printf("\t filling histo hPedhg[%d]\n",index);
 311             }
 312             else{
 313               hPedlg[index]->Fill(rawStreamZDC->GetADCValue());
 314               RawADClg[index] = rawStreamZDC->GetADCValue();
 315               //
 316               //printf("\t filling histo hPedlg[%d]\n",index);
 317             }
 318            }
 319            else{  // --- Out-of-time pedestals
 320             if(rawStreamZDC->GetADCGain()==0){
 321               hPedOutOfTimehg[index]->Fill(rawStreamZDC->GetADCValue());
 322               RawADCoothg[index] = rawStreamZDC->GetADCValue();
 323               //
 324               //printf("\t filling histo hPedOutOfTimehg[%d]\n",index);
 325             }
 326             else{
 327               hPedOutOfTimelg[index]->Fill(rawStreamZDC->GetADCValue());
 328               RawADCootlg[index] = rawStreamZDC->GetADCValue();
 329               //
 330               //printf("\t filling histo hPedOutOfTimelg[%d]\n",index);
 331             }
 332            }
 333             iraw++;
 334           }//IsADCDataWord()
 335           //
 336           if(iraw == 4*kNChannels){ // Last ADC channel -> Filling correlation histos
 337             for(Int_t k=0; k<kNChannels; k++){
 338               hPedCorrhg[k]->Fill(RawADCoothg[k], RawADChg[k]);
 339               hPedCorrlg[k]->Fill(RawADCootlg[k], RawADClg[k]);
 340             }
 341           }
 342          }
 343          //
 344          nevents_physics++;
 345          //
 346          delete reader;
 347          delete rawStreamZDC;
 348
 349       }//(if PHYSICS_EVENT)
 350       nevents_total++;
 351
 352       /* free resources */
 353       free(event);
 354
 355     }
 356   }
 357
 358   /* Analysis of the histograms */
 359   //
 360   FILE *fileShuttle;
 361   const char *fName = "ZDCPedestal.dat";
 362   fileShuttle = fopen(fName,"w");
 363   //
 364   Float_t MeanPed[2*kNChannels], MeanPedWidth[2*kNChannels],
 365         MeanPedOOT[2*kNChannels], MeanPedWidthOOT[2*kNChannels],
 366         CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
 367   // --- In-time pedestals
 368   TF1 *ADCfunchg[kNChannels];
 369   for(Int_t i=0; i<kNChannels; i++){
 370      hPedhg[i]->Fit("gaus","Q");
 371      ADCfunchg[i] = hPedhg[i]->GetFunction("gaus");
 372      MeanPed[i] = (Double_t) ADCfunchg[i]->GetParameter(1);
 373      MeanPedWidth[i] = (Double_t)  ADCfunchg[i]->GetParameter(2);
 374      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,MeanPed[i],MeanPedWidth[i]);
 375      //printf("\t MeanPed[%d] = %f\n",i, MeanPed[i]);
 376   }
 377   TF1 *ADCfunclg[kNChannels];
 378   for(Int_t i=0; i<kNChannels; i++){
 379      hPedlg[i]->Fit("gaus","Q");
 380      ADCfunclg[i] = hPedlg[i]->GetFunction("gaus");
 381      MeanPed[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(1);
 382      MeanPedWidth[i+kNChannels] = (Double_t)  ADCfunclg[i]->GetParameter(2);
 383      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,MeanPed[i+kNChannels],MeanPedWidth[i+kNChannels]);
 384      //printf("\t MeanPed[%d] = %f\n",i+kNChannels, MeanPed[i+kNChannels]);
 385   }
 386   // --- Out-of-time pedestals
 387   TF1 *ADCootfunchg[kNChannels];
 388   for(Int_t i=0; i<kNChannels; i++){
 389      hPedOutOfTimehg[i]->Fit("gaus","Q");
 390      ADCootfunchg[i] = hPedOutOfTimehg[i]->GetFunction("gaus");
 391      MeanPedOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(1);
 392      MeanPedWidthOOT[i] = (Double_t)  ADCootfunchg[i]->GetParameter(2);
 393      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,MeanPedOOT[i],MeanPedWidthOOT[i]);
 394      //printf("\t MeanPedOOT[%d] = %f\n",i, MeanPedOOT[i]);
 395   }
 396   TF1 *ADCootfunclg[kNChannels];
 397   for(Int_t i=0; i<kNChannels; i++){
 398      hPedOutOfTimelg[i]->Fit("gaus","Q");
 399      ADCootfunclg[i] = hPedOutOfTimelg[i]->GetFunction("gaus");
 400      MeanPedOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(1);
 401      MeanPedWidthOOT[i+kNChannels] = (Double_t)  ADCootfunclg[i]->GetParameter(2);
 402      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,MeanPedOOT[i+kNChannels],MeanPedWidthOOT[i+kNChannels]);
 403      //printf("\t MeanPedOOT[%d] = %f\n",i+kNChannels, MeanPedOOT[i+kNChannels]);
 404   }
 405   //
 406   // --- Correlations
 407   TProfile *hPedCorrProfhg[kNChannels], *hPedCorrProflg[kNChannels];
 408   TF1 *ffunchg[kNChannels], *ffunclg[kNChannels];
 409   char namhist4[50];
 410   for(int i=0;i<kNChannels;i++) {
 411      sprintf(namhist4,"ADCHRvsOOT%d_Prof",i);
 412      hPedCorrProfhg[i] = hPedCorrhg[i]->ProfileX(namhist4,-1,-1,"S");
 413      hPedCorrProfhg[i]->SetName(namhist4);
 414      hPedCorrProfhg[i]->Fit("pol1","Q");
 415      ffunchg[i] = hPedCorrProfhg[i]->GetFunction("pol1");
 416      CorrCoeff0[i] = (Double_t)  ffunchg[i]->GetParameter(0);
 417      CorrCoeff1[i] = (Double_t) ffunchg[i]->GetParameter(1);
 418      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i,CorrCoeff0[i],CorrCoeff1[i]);
 419      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",i, CorrCoeff0[i], i, CorrCoeff1[i]);
 420   }
 421   for(int i=0;i<kNChannels;i++) {
 422      sprintf(namhist4,"ADCLRvsOOT%d_Prof",i);
 423      hPedCorrProflg[i] = hPedCorrlg[i]->ProfileX(namhist4,-1,-1,"S");
 424      hPedCorrProflg[i]->SetName(namhist4);
 425      hPedCorrProflg[i]->Fit("pol1","Q");
 426      ffunclg[i] = hPedCorrProflg[i]->GetFunction("pol1");
 427      CorrCoeff0[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(0);
 428      CorrCoeff1[i+kNChannels] =  (Double_t) ffunclg[i]->GetParameter(1);
 429      fprintf(fileShuttle,"\t%d\t%f\t%f\n",i+kNChannels,CorrCoeff0[i+kNChannels],CorrCoeff1[i+kNChannels]);
 430      //printf("\t CorrCoeff0[%d] = %f, CorrCoeff1[%d] = %f\n",
 431      //         i+kNChannels, CorrCoeff0[i+kNChannels], i+kNChannels, CorrCoeff1[i+kNChannels]);
 432   }
 433   //
 434   fclose(fileShuttle);
 435   //
 436   for(Int_t j=0; j<kNChannels; j++){
 437      delete hPedhg[j];
 438      delete hPedOutOfTimehg[j];
 439      delete hPedCorrhg[j];
 440      delete hPedlg[j];
 441      delete hPedOutOfTimelg[j];
 442      delete hPedCorrlg[j];
 443   }
 444
 445   //delete minuitFit;
 446   TVirtualFitter::SetFitter(0);
 447
 448   /* write report */
 449   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
 450
 451   /* close result file */
 452   fclose(fp);
 453
 454   /* report progress */
 455   daqDA_progressReport(90);
 456
 457   /* store the result files on FES */
 458   status = daqDA_FES_storeFile(mapfName,"ZDCCHMAPPING_data");
 459   if(status){
 460     printf("Failed to export file : %d\n",status);
 461     return -1;
 462   }
 463   //
 464   status = daqDA_FES_storeFile(fName,"ZDCPEDESTAL_data");
 465   if(status){
 466     printf("Failed to export file : %d\n",status);
 467     return -1;
 468   }
 469
 470   /* report progress */
 471   daqDA_progressReport(100);
 472
 473   return status;
 474 }