VZERO/VZEROPbPbda.cxx

   1 /*********************************************************************************
   2 - Contact:    Brigitte Cheynis     b.cheynis@ipnl.in2p3.fr
   3 - Link:
   4 - Raw data test file :
   5 - Reference run number : 137366
   6 - Run Type:   PHYSICS
   7 - DA Type:    MON
   8 - Number of events needed: >=2000
   9 - Input Files:  argument list
  10 - Output Files: FXS file     V0_EqualizationFactors.dat (Channel equalization factors)
  11 - Trigger types used: PHYSICS_EVENT
  12 **********************************************************************************/
  13
  14 /**********************************************************************************
  15 *                                                                                 *
  16 * VZERO Detector Algorithm for extracting channel equalization factors for Pb-Pb  *
  17 *                                                                                 *
  18 *                                                                                 *
  19 ***********************************************************************************/
  20
  21 // DATE
  22 #include "event.h"
  23 #include "monitor.h"
  24 #include "daqDA.h"
  25
  26 //AliRoot
  27 #include <AliVZERORawStream.h>
  28 #include <AliRawReaderDate.h>
  29 #include <AliRawReader.h>
  30 #include <AliDAQ.h>
  31
  32 // standard
  33 #include <stdio.h>
  34 #include <stdlib.h>
  35
  36 //ROOT
  37 #include "TROOT.h"
  38 #include "TPluginManager.h"
  39 #include <TFile.h>
  40 #include <TH1F.h>
  41 #include <TMath.h>
  42
  43 Int_t GetOfflineChannel(Int_t channel);
  44
  45 /* Main routine --- Arguments: monitoring data source */
  46
  47 int main(int argc, char **argv) {
  48
  49 /* magic line from Cvetan */
  50   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
  51                     "*",
  52                     "TStreamerInfo",
  53                     "RIO",
  54                     "TStreamerInfo()");
  55   int status;
  56   if (argc!=2) {
  57      printf("Wrong number of arguments\n");
  58      return -1;
  59   }
  60
  61 //___________________________________________________
  62 // Get parameters from V00DAEqualFactors.config file
  63   Int_t    neventsMin = 2000;
  64
  65   Int_t    kStartClock = 9;  // First clock in the search for max adc
  66   Int_t    kEndClock = 11;   // Last clock in the search for max adc
  67   Int_t    kNPreClocks = 6;  // Number of clock before max used in the charge sum
  68   Int_t    kNPostClocks = 1; // Number of clock after max used in the charge sum
  69
  70   UChar_t  kNFlagsCut   = 63;// Minimum number of TDC flags requested, replaces the trigger selection
  71
  72   Int_t    kNBins = 10000;
  73   Float_t  kRange = 0.1;
  74
  75   status = daqDA_DB_getFile("V00DAEqualFactors.config","./V00DAEqualFactors.config");
  76   if (status) {
  77     printf("Failed to get Config file (V00DAEqualFactors.config) from DAQ DB, status=%d\n", status);
  78     printf("Take default values of parameters for pedestal calculation \n");
  79   } else {
  80     /* open the config file and retrieve cuts */
  81     FILE *fpConfig = fopen("V00DAEqualFactors.config","r");
  82     int res = fscanf(fpConfig,"%d %d %d %d %d %hhu %d %f",
  83                      &neventsMin,&kStartClock,&kEndClock,&kNPreClocks,&kNPostClocks,&kNFlagsCut,&kNBins,&kRange);
  84     if(res!=8) {
  85       printf("Failed to get values from Config file (V00DAEqualFactors.config): wrong file format - 7 integers and 1 float are expected - \n");
  86     }
  87     fclose(fpConfig);
  88   }
  89
  90   printf("Minimum number of events to process = %d; First LHC Clock = %d; Last LHC Clock = %d; N Pre Clock = %d; N Post Clock = %d; Minimum number of TDC flags = %hhu; Number of histogram bins = %d; Histogram range = %.3f\n",
  91          neventsMin,kStartClock, kEndClock, kNPreClocks, kNPostClocks, kNFlagsCut, kNBins, kRange);
  92
  93   TH1D *fMedian[64];
  94   for(Int_t j = 0; j < 64; ++j) fMedian[j] = new TH1D(Form("fMedian_%d",j),"Slopes weighted median, channel par channel",kNBins,0,kRange);
  95
  96   Bool_t fFirst = kTRUE;
  97   Float_t fPrevTotCharge = 0;
  98   Float_t fPrevadc[64];
  99   for(Int_t j = 0; j < 64; ++j) fPrevadc[j] = 0;
 100
 101 //___________________________________________________
 102
 103   /* define data source : this is argument 1 */
 104   status=monitorSetDataSource( argv[1] );
 105   if (status!=0) {
 106      printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
 107      return -1;
 108   }
 109
 110   /* declare monitoring program */
 111   status=monitorDeclareMp( __FILE__ );
 112   if (status!=0) {
 113     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
 114     return -1;
 115   }
 116
 117   /* define wait event timeout - 1s max */
 118   monitorSetNowait();
 119   monitorSetNoWaitNetworkTimeout(1000);
 120
 121   /* init counters on events */
 122   int neventsPhysicsAll=0;
 123   int neventsPhysics=0;
 124   int neventsTotal=0;
 125
 126
 127   /* loop on events (infinite) */
 128   for(;;) {
 129       struct eventHeaderStruct *event;
 130       eventTypeType eventT;
 131
 132       /* check shutdown condition */
 133       if (daqDA_checkShutdown()) {break;}
 134
 135       /* get next event (blocking call until timeout) */
 136       status=monitorGetEventDynamic((void **)&event);
 137       if (status==MON_ERR_EOF) {
 138           printf ("End of File detected\n");
 139       break; /* end of monitoring file has been reached */
 140       }
 141
 142       if (status!=0) {
 143            printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 144       break;
 145       }
 146
 147       /* retry if got no event */
 148       if (event==NULL) continue;
 149
 150       /* decode event */
 151       eventT=event->eventType;
 152
 153       switch (event->eventType){
 154
 155       case START_OF_RUN:
 156            break;
 157
 158       case END_OF_RUN:
 159            printf("End Of Run detected\n");
 160            break;
 161
 162       case PHYSICS_EVENT:
 163
 164            AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 165
 166            AliVZERORawStream* rawStream  = new AliVZERORawStream(rawReader);
 167            if (rawStream->Next()) {
 168              neventsPhysicsAll++;
 169              UShort_t triggers = rawStream->GetTriggerInputs();
 170              UChar_t nFlags = 0;
 171              for(Int_t i = 0; i < 64; ++i) {
 172                if (rawStream->GetBBFlag(i,10)) nFlags++;
 173              }
 174              if (nFlags >= kNFlagsCut) { // Check if the number of minimum number of TDC flags is reached (trigger selection)
 175
 176                neventsPhysics++;
 177
 178                Float_t adc[64];
 179                Float_t totCharge = 0.;
 180                for(Int_t i = 0; i < 64; ++i) {
 181                  adc[i] = 0.;
 182                  Float_t maxadc=0.;
 183                  Int_t imax=-1;
 184                  for(Int_t j = kStartClock; j <= kEndClock; ++j) {
 185                    Float_t charge = (Float_t)(rawStream->GetPedestal(i,j));
 186                    if(charge > maxadc){
 187                      maxadc = charge;
 188                      imax   = j;
 189                    }
 190                  }
 191
 192                  if (imax != -1) {
 193                    Int_t start = imax - kNPreClocks;
 194                    if (start < 0) start = 0;
 195                    Int_t end = imax + kNPostClocks;
 196                    if (end > 20) end = 20;
 197                    for(Int_t iClock = start; iClock <= end; iClock++) {
 198                      adc[i] += (Float_t)(rawStream->GetPedestal(i,iClock));
 199                    }
 200                  }
 201                  totCharge += adc[i];
 202                }
 203
 204                if (fFirst) {
 205                  fFirst = kFALSE;
 206                  fPrevTotCharge = totCharge;
 207                  for(int i = 0; i < 64; ++i) fPrevadc[i] = adc[i];
 208                }
 209                else {
 210                  fFirst = kTRUE;
 211                  Float_t deltaTotCharge = totCharge - fPrevTotCharge;
 212                  Float_t weight = deltaTotCharge*deltaTotCharge;
 213                  if (weight > 1) {
 214                    for(int i = 0; i < 64; ++i) {
 215                      fMedian[i]->Fill((adc[i]-fPrevadc[i])/deltaTotCharge,weight);
 216                    }
 217                  }
 218                }
 219
 220              }
 221            }   // End : if rawstream
 222            delete rawStream;
 223            rawStream = 0x0;
 224            delete rawReader;
 225            rawReader = 0x0;
 226       } // end of switch on event type
 227
 228       neventsTotal++;
 229       /* free resources */
 230       free(event);
 231
 232       /* exit when last event received, no need to wait for TERM signal */
 233       if (eventT==END_OF_RUN) {
 234         printf("End Of Run event detected\n");
 235         break;
 236       }
 237
 238   }  // loop over events
 239
 240   printf("%d physics events processed (out of %d physics and %d total events)\n",neventsPhysics,neventsPhysicsAll,neventsTotal);
 241
 242 //___________________________________________________________________________
 243 //  Computes regression parameters
 244 // charge_i = p0 + charge_tot * p1
 245
 246   if(neventsPhysics>neventsMin){
 247     /* open result file to be exported to FES */
 248     FILE *fp=NULL;
 249     fp=fopen("./V0_EqualizationFactors.dat","w");
 250     if (fp==NULL) {
 251       printf("Failed to open local result file\n");
 252       return -1;}
 253
 254     Double_t beta[64];
 255     Double_t q = 0.5;
 256     for(int i = 0; i < 64; ++i) fMedian[i]->GetQuantiles(1,&beta[i],&q);
 257
 258     for(Int_t i=0; i<64; i++) {
 259       fprintf(fp," %d %.5f\n",GetOfflineChannel(i), beta[i]*64.);
 260       printf(" %d %.5f\n",GetOfflineChannel(i), beta[i]*64.);
 261     }
 262
 263     /* close local result file and FXS result file*/
 264     fclose(fp);
 265
 266     /* export result file to FES */
 267     status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
 268     if (status)    {
 269       printf("Failed to export file : %d\n",status);
 270       return -1; }
 271
 272     /* store result file into Online DB */
 273     status=daqDA_DB_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
 274     if (status)    {
 275       printf("Failed to store file into Online DB: %d\n",status);
 276       return -1; }
 277   }
 278   else {
 279     // Take the last run's file from the DB
 280     status=daqDA_DB_getFile("V00DAEqualFactors","./V0_EqualizationFactors.dat");
 281     if (status)    {
 282       printf("Failed to get file from Online DB: %d\n",status);
 283       return -1; }
 284
 285     /* export result file to FES */
 286     status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
 287     if (status)    {
 288       printf("Failed to export file : %d\n",status);
 289       return -1; }
 290   }
 291
 292   return status;
 293 }
 294
 295  Int_t GetOfflineChannel(Int_t channel) {
 296
 297 // Channel mapping Online - Offline:
 298
 299  Int_t fOfflineChannel[64] = {39, 38, 37, 36, 35, 34, 33, 32,
 300                               47, 46, 45, 44, 43, 42, 41, 40,
 301                               55, 54, 53, 52, 51, 50, 49, 48,
 302                               63, 62, 61, 60, 59, 58, 57, 56,
 303                                7,  6,  5,  4,  3,  2,  1,  0,
 304                               15, 14, 13, 12, 11, 10,  9,  8,
 305                               23, 22, 21, 20, 19, 18, 17, 16,
 306                               31, 30, 29, 28, 27, 26, 25, 24};
 307  return fOfflineChannel[channel];
 308 }