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