ZDC/ZDCEMDda.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 CALIBRATION_EMD runs
  13
  14 Contact: Chiara.Oppedisano@to.infn.it
  15 Link:
  16 Run Type: STANDALONE_EMD_RUN
  17 DA Type: LDC
  18 Number of events needed: at least ~5*10^3
  19 Input Files: ZDCPedestal.dat
  20 Output Files: ZDCEMDCalib.dat, ZDCChMapping.dat
  21 Trigger Types Used: Standalone Trigger
  22
  23 */
  24 #define PEDDATA_FILE  "ZDCPedestal.dat"
  25 #define MAPDATA_FILE  "ZDCChMapping.dat"
  26 #define ENCALIBDATA_FILE   "ZDCEnergyCalib.dat"
  27 #define TOWCALIBDATA_FILE  "ZDCTowerCalib.dat"
  28
  29 #include <stdio.h>
  30 #include <Riostream.h>
  31 #include <Riostream.h>
  32
  33 // DATE
  34 #include <daqDA.h>
  35 #include <event.h>
  36 #include <monitor.h>
  37
  38 //ROOT
  39 #include <TROOT.h>
  40 #include <TPluginManager.h>
  41 #include <TH1F.h>
  42 #include <TH2F.h>
  43 #include <TProfile.h>
  44 #include <TF1.h>
  45 #include <TFile.h>
  46 #include <TFitter.h>
  47 #include "TMinuitMinimizer.h"
  48
  49 //AliRoot
  50 #include <AliRawReaderDate.h>
  51 #include <AliRawEventHeaderBase.h>
  52 #include <AliZDCRawStream.h>
  53
  54
  55 /* Main routine
  56       Arguments:
  57       1- monitoring data source
  58 */
  59 int main(int argc, char **argv) {
  60
  61
  62   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
  63                                         "*",
  64                                         "TStreamerInfo",
  65                                         "RIO",
  66                                         "TStreamerInfo()");
  67
  68   TMinuitMinimizer m;
  69   gROOT->GetPluginManager()->AddHandler("ROOT::Math::Minimizer", "Minuit","TMinuitMinimizer",
  70       "Minuit", "TMinuitMinimizer(const char *)");
  71   TVirtualFitter::SetDefaultFitter("Minuit");
  72
  73   int status = 0;
  74   // No. of ZDC cabled ch.
  75   int const kNChannels = 24;
  76
  77   /* log start of process */
  78   printf("ZDC EMD program started\n");
  79
  80   /* check that we got some arguments = list of files */
  81   if (argc<2) {
  82     printf("Wrong number of arguments\n");
  83     return -1;
  84   }
  85
  86   //
  87   // --- Preparing histos for EM dissociation spectra
  88   //
  89   TH1F* histoEMDRaw[4];
  90   TH1F* histoEMDCorr[4];
  91
  92   char namhistr[50], namhistc[50];
  93   for(Int_t i=0; i<4; i++) {
  94      if(i==0){
  95        sprintf(namhistr,"ZN%d-EMDRaw",i+1);
  96        sprintf(namhistc,"ZN%d-EMDCorr",i+1);
  97      }
  98      else if(i==1){
  99        sprintf(namhistr,"ZP%d-EMDRaw",i);
 100        sprintf(namhistc,"ZP%d-EMDCorr",i);
 101      }
 102      else if(i==2){
 103        sprintf(namhistr,"ZN%d-EMDRaw",i);
 104        sprintf(namhistc,"ZN%d-EMDCorr",i);
 105      }
 106      else if(i==3){
 107        sprintf(namhistr,"ZP%d-EMDRaw",i-1);
 108        sprintf(namhistc,"ZP%d-EMDCorr",i-1);
 109      }
 110      histoEMDRaw[i] = new TH1F(namhistr,namhistr,100,0.,4000.);
 111      histoEMDCorr[i] = new TH1F(namhistc,namhistc,100,0.,4000.);
 112   }
 113
 114   /* open result file */
 115   FILE *fp=NULL;
 116   fp=fopen("./result.txt","a");
 117   if (fp==NULL) {
 118     printf("Failed to open file\n");
 119     return -1;
 120   }
 121
 122   FILE *mapFile4Shuttle;
 123
 124   // *** To analyze EMD events you MUST have a pedestal data file!!!
 125   // *** -> check if a pedestal run has been analyzed
 126   int read = 0;
 127   read = daqDA_DB_getFile(PEDDATA_FILE,PEDDATA_FILE);
 128   if(read){
 129     printf("\t ERROR!!! ZDCPedestal.dat file NOT FOUND in DAQ db!!!\n");
 130     return -1;
 131   }
 132   else printf("\t ZDCPedestal.dat file retrieved from DAQ db\n");
 133
 134   FILE *filePed = fopen(PEDDATA_FILE,"r");
 135   if (filePed==NULL) {
 136     printf("\t ERROR!!! Can't open ZDCPedestal.dat file!!!\n");
 137     return -1;
 138   }
 139
 140   // 144 = 48 in-time + 48 out-of-time + 48 correlations
 141   Float_t readValues[2][3*2*kNChannels];
 142   Float_t MeanPed[2*kNChannels];
 143   Float_t CorrCoeff0[2*kNChannels], CorrCoeff1[2*kNChannels];
 144   // ***************************************************
 145   //   Unless we have a narrow correlation to fit we
 146   //    don't fit and store in-time vs. out-of-time
 147   //    histograms -> mean pedstal subtracted!!!!!!
 148   // ***************************************************
 149   //
 150   for(int jj=0; jj<6*kNChannels; jj++){
 151     for(int ii=0; ii<2; ii++){
 152        fscanf(filePed,"%f",&readValues[ii][jj]);
 153     }
 154     if(jj<kNChannels && jj<2*kNChannels){
 155       MeanPed[jj] = readValues[0][jj];
 156       //printf("\t MeanPedhg[%d] = %1.1f\n",jj, MeanPedhg[jj]);
 157     }
 158     else if(jj>2*kNChannels && jj>4*kNChannels){
 159       CorrCoeff0[jj-4*kNChannels] = readValues[0][jj];
 160       CorrCoeff1[jj-4*kNChannels] = readValues[1][jj];;
 161     }
 162   }
 163
 164   /* report progress */
 165   daqDA_progressReport(10);
 166
 167
 168   /* init some counters */
 169   int nevents_physics=0;
 170   int nevents_total=0;
 171
 172   /* read the data files */
 173   int n;
 174   for(n=1;n<argc;n++){
 175
 176     status=monitorSetDataSource( argv[n] );
 177     if (status!=0) {
 178       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
 179       return -1;
 180     }
 181
 182     /* report progress */
 183     /* in this example, indexed on the number of files */
 184     daqDA_progressReport(10+80*n/argc);
 185
 186     /* read the file */
 187     for(;;){
 188       struct eventHeaderStruct *event;
 189       eventTypeType eventT;
 190
 191       /* get next event */
 192       status=monitorGetEventDynamic((void **)&event);
 193       if(status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
 194       if(status!=0) {
 195         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 196         return -1;
 197       }
 198
 199       /* retry if got no event */
 200       if(event==NULL) {
 201         break;
 202       }
 203
 204       // Initalize raw-data reading and decoding
 205       AliRawReader *reader = new AliRawReaderDate((void*)event);
 206       reader->Select("ZDC");
 207       // --- Reading event header
 208       //UInt_t evtype = reader->GetType();
 209       //printf("\n\t ZDCEMDda -> ev. type %d\n",evtype);
 210       //printf("\t ZDCEMDda -> run # %d\n",reader->GetRunNumber());
 211       //
 212       AliZDCRawStream *rawStreamZDC = new AliZDCRawStream(reader);
 213
 214
 215       /* use event - here, just write event id to result file */
 216       eventT=event->eventType;
 217
 218       Int_t ich=0;
 219       Int_t adcMod[2*kNChannels], adcCh[2*kNChannels];
 220       Int_t sigCode[2*kNChannels], det[2*kNChannels], sec[2*kNChannels];
 221
 222       if(eventT==START_OF_DATA){
 223
 224         rawStreamZDC->SetSODReading(kTRUE);
 225
 226         // --------------------------------------------------------
 227         // --- Writing ascii data file for the Shuttle preprocessor
 228         mapFile4Shuttle = fopen(MAPDATA_FILE,"w");
 229         if(!rawStreamZDC->Next()) printf(" \t No raw data found!! \n");
 230         else{
 231           while(rawStreamZDC->Next()){
 232             if(rawStreamZDC->IsChMapping()){
 233               adcMod[ich] = rawStreamZDC->GetADCModFromMap(ich);
 234               adcCh[ich] = rawStreamZDC->GetADCChFromMap(ich);
 235               sigCode[ich] = rawStreamZDC->GetADCSignFromMap(ich);
 236               det[ich] = rawStreamZDC->GetDetectorFromMap(ich);
 237               sec[ich] = rawStreamZDC->GetTowerFromMap(ich);
 238               //
 239               fprintf(mapFile4Shuttle,"\t%d\t%d\t%d\t%d\t%d\t%d\n",
 240                 ich,adcMod[ich],adcCh[ich],sigCode[ich],det[ich],sec[ich]);
 241               //
 242               //printf("ZDCEMDDA.cxx ->  ch.%d mod %d, ch %d, code %d det %d, sec %d\n",
 243               //           ich,adcMod[ich],adcCh[ich],sigCode[ich],det[ich],sec[ich]);
 244               ich++;
 245             }
 246           }
 247         }
 248         fclose(mapFile4Shuttle);
 249       }//SOD event
 250
 251
 252     if(eventT==PHYSICS_EVENT){
 253       // --- Reading data header
 254       reader->ReadHeader();
 255       const AliRawDataHeader* header = reader->GetDataHeader();
 256       if(header){
 257          UChar_t message = header->GetAttributes();
 258          if(message & 0x70){ // DEDICATED EMD RUN
 259             //printf("\t STANDALONE_EMD_RUN raw data found\n");
 260             continue;
 261          }
 262          else{
 263             printf("\t NO STANDALONE_EMD_RUN raw data found\n");
 264             return -1;
 265          }
 266       }
 267       else{
 268          printf("\t ATTENTION! No Raw Data Header found!!!\n");
 269          return -1;
 270       }
 271
 272       rawStreamZDC->SetSODReading(kTRUE);
 273
 274       if (!rawStreamZDC->Next()) printf(" \t No raw data found!! ");
 275       //
 276       // ----- Setting ch. mapping -----
 277       for(Int_t jk=0; jk<2*kNChannels; 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       Float_t ZDCRawADC[4], ZDCCorrADC[4], ZDCCorrADCSum[4];
 286       for(Int_t g=0; g<4; g++){
 287            ZDCCorrADCSum[g] = 0.;
 288            ZDCRawADC[g] = 0.;
 289       }
 290       //
 291       while(rawStreamZDC->Next()){
 292         Int_t det = rawStreamZDC->GetSector(0);
 293         Int_t quad = rawStreamZDC->GetSector(1);
 294
 295         if(rawStreamZDC->IsADCDataWord() && !(rawStreamZDC->IsUnderflow())
 296              && !(rawStreamZDC->IsOverflow()) && det!=-1
 297              && (rawStreamZDC->GetADCGain() == 1)){ // Selecting LOW RES ch.s
 298
 299           //printf("  IsADCWord %d, IsUnderflow %d, IsOverflow %d\n",
 300           //  rawStreamZDC->IsADCDataWord(),rawStreamZDC->IsUnderflow(),rawStreamZDC->IsOverflow());
 301
 302           // Taking LOW RES channels -> channel+kNChannels !!!!
 303           Int_t DetIndex=999, PedIndex=999;
 304           if(det != 3 && quad != 5){ // Not ZEM nor PMRef
 305             if(det == 1){
 306               DetIndex = det-1;
 307               PedIndex = quad+kNChannels;
 308             }
 309             else if(det==2){
 310               DetIndex = det-1;
 311               PedIndex = quad+5+kNChannels;
 312             }
 313             else if(det == 4){
 314               DetIndex = det-2;
 315               PedIndex = quad+12+kNChannels;
 316             }
 317             else if(det == 5){
 318               DetIndex = det-2;
 319               PedIndex = quad+17+kNChannels;
 320             }
 321             //EMD -> LR ADCs
 322             if(rawStreamZDC->GetADCGain() == 1 && (DetIndex!=999 || PedIndex!=999)){
 323               //
 324               ZDCRawADC[DetIndex] += (Float_t) rawStreamZDC->GetADCValue();
 325               //
 326               // Mean pedestal subtraction
 327               Float_t Pedestal = MeanPed[PedIndex];
 328               // Pedestal subtraction from correlation with out-of-time signals
 329               //Float_t Pedestal = CorrCoeff0[PedIndex]+CorrCoeff1[PedIndex]*MeanPedOOT[PedIndex];
 330               //
 331               ZDCCorrADC[DetIndex] = (rawStreamZDC->GetADCValue()) - Pedestal;
 332               ZDCCorrADCSum[DetIndex] += ZDCCorrADC[DetIndex];
 333               //
 334               /*printf("\t det %d quad %d res %d pedInd %d detInd %d:"
 335                  "ADCCorr = %d, ZDCCorrADCSum = %d\n",
 336                  det,quad,rawStreamZDC->GetADCGain(),PedIndex,DetIndex,
 337                  (Int_t) ZDCCorrADC[DetIndex],(Int_t) ZDCCorrADCSum[DetIndex]);
 338               */
 339             }
 340             if(DetIndex==999 || PedIndex==999)
 341                 printf(" WARNING! Detector a/o pedestal index are WRONG!!!\n");
 342
 343           }
 344         }//IsADCDataWord()
 345          //
 346        }
 347        //
 348        nevents_physics++;
 349        //
 350        for(Int_t j=0; j<4; j++){
 351           histoEMDRaw[j]->Fill(ZDCRawADC[j]);
 352           histoEMDCorr[j]->Fill(ZDCCorrADCSum[j]);
 353        }
 354     }
 355
 356     nevents_total++;
 357
 358     /* free resources */
 359     free(event);
 360
 361     /* exit when last event received, no need to wait for TERM signal */
 362     if (eventT==END_OF_RUN) {
 363       printf("EOR event detected\n");
 364       break;
 365     }
 366
 367    }
 368   }
 369
 370   /* Analysis of the histograms */
 371   //
 372   FILE *fileShuttle1 = fopen(ENCALIBDATA_FILE,"w");
 373   //
 374   Int_t BinMax[4];
 375   Float_t YMax[4];
 376   Int_t NBinsx[4];
 377   Float_t MeanFitVal[4];
 378   TF1 *fitfun[4];
 379   for(Int_t k=0; k<4; k++){
 380      if(histoEMDCorr[k]->GetEntries() == 0){
 381        printf("\n WARNING! Empty histos -> ending DA WITHOUT writing output\n\n");
 382        return -1;
 383      }
 384      BinMax[k] = histoEMDCorr[k]->GetMaximumBin();
 385      YMax[k] = (histoEMDCorr[k]->GetXaxis())->GetXmax();
 386      NBinsx[k] = (histoEMDCorr[k]->GetXaxis())->GetNbins();
 387 //     printf("\n\t Det%d -> BinMax = %d, ChXMax = %f\n", k+1, BinMax[k], BinMax[k]*YMax[k]/NBinsx[k]);
 388      histoEMDCorr[k]->Fit("gaus","Q","",BinMax[k]*YMax[k]/NBinsx[k]*0.7,BinMax[k]*YMax[k]/NBinsx[k]*1.25);
 389      fitfun[k] = histoEMDCorr[k]->GetFunction("gaus");
 390      MeanFitVal[k] = (Float_t) (fitfun[k]->GetParameter(1));
 391      //printf("\n\t Mean Value from gaussian fit = %f\n", MeanFitVal[k]);
 392   }
 393   //
 394   Float_t CalibCoeff[6];
 395   Float_t icoeff[5];
 396   //
 397   for(Int_t j=0; j<6; j++){
 398      if(j<4){
 399        CalibCoeff[j] = MeanFitVal[j];
 400        fprintf(fileShuttle1,"\t%f\n",CalibCoeff[j]);
 401      }
 402      // ZEM energy calib. coeff. = 1
 403      else if(j==4 || j==5){
 404        CalibCoeff[j] = 1.;
 405        fprintf(fileShuttle1,"\t%f\n",CalibCoeff[j]);
 406      }
 407   }
 408   fclose(fileShuttle1);
 409   //
 410   FILE *fileShuttle2 = fopen(TOWCALIBDATA_FILE,"w");
 411   for(Int_t j=0; j<4; j++){
 412      // Note -> For the moment the inter-calibration coeff. are set to 1
 413      for(Int_t k=0; k<5; k++){
 414        icoeff[k] = 1.;
 415        fprintf(fileShuttle2,"\t%f",icoeff[k]);
 416        if(k==4) fprintf(fileShuttle2,"\n");
 417      }
 418   }
 419   fclose(fileShuttle2);
 420
 421   for(Int_t ij=0; ij<4; ij++){
 422     delete histoEMDRaw[ij];
 423     delete histoEMDCorr[ij];
 424   }
 425
 426   //delete minuitFit;
 427   TVirtualFitter::SetFitter(0);
 428
 429   /* write report */
 430   fprintf(fp,"Run #%s, received %d physics events out of %d\n",getenv("DATE_RUN_NUMBER"),nevents_physics,nevents_total);
 431
 432   /* close result file */
 433   fclose(fp);
 434
 435   /* report progress */
 436   daqDA_progressReport(90);
 437
 438   /* store the result file on FES */
 439   status = daqDA_FES_storeFile(MAPDATA_FILE, MAPDATA_FILE);
 440   if(status){
 441     printf("Failed to export file : %d\n",status);
 442     return -1;
 443   }
 444   //
 445   status = daqDA_FES_storeFile(ENCALIBDATA_FILE, ENCALIBDATA_FILE);
 446   if(status){
 447     printf("Failed to export file : %d\n",status);
 448     return -1;
 449   }
 450   //
 451   status = daqDA_FES_storeFile(TOWCALIBDATA_FILE, TOWCALIBDATA_FILE);
 452   if(status){
 453     printf("Failed to export file : %d\n",status);
 454     return -1;
 455   }
 456
 457   /* report progress */
 458   daqDA_progressReport(100);
 459
 460   return status;
 461 }