HMPID/HMPIDda.cxx

   1 /*
   2 *********************************************************
   3 Author:                                                 *
   4 this file provides the detector algorithm for HMPID.    *
   5 *********************************************************
   6 */
   7 extern "C" {
   8 #include <daqDA.h>
   9 }
  10
  11 #include "event.h"
  12 #include "monitor.h"
  13
  14 #include <Riostream.h>
  15 #include <stdio.h>
  16 #include <stdlib.h>
  17
  18 //AliRoot
  19 #include "AliHMPIDRawStream.h"
  20 #include "AliRawReaderDate.h"
  21 #include "AliBitPacking.h"
  22 #include "TMath.h"
  23
  24 //ROOT
  25 #include "TFile.h"
  26 #include "TKey.h"
  27 #include "TH2S.h"
  28 #include "TObject.h"
  29 #include "TBenchmark.h"
  30 #include "TMath.h"
  31 #include "TRandom.h"
  32
  33 int main(int argc, char **argv){
  34
  35   int status;
  36
  37   /* log start of process */
  38   printf("HMPID DA program started\n");
  39
  40   /* check that we got some arguments = list of files */
  41   if (argc<2) {
  42     printf("Wrong number of arguments\n");
  43     return -1;
  44   }
  45
  46   /* copy locally a file from daq detector config db
  47   status=daqDA_DB_getFile("myconfig","./myconfig.txt");
  48   if (status) {
  49     printf("Failed to get config file : %d\n",status);
  50     return -1;
  51   }
  52   and possibly use it */
  53
  54   /* report progress */
  55   daqDA_progressReport(10);
  56
  57
  58   /* init some counters */
  59   Double_t SummQ[14][48][11][25], SummQ2[14][48][11][25];
  60   Int_t   nEntries[14][48][11][25];
  61   Bool_t  isDDLOn[14];
  62   for(Int_t ddl=0;ddl<=13;ddl++) {
  63     isDDLOn[ddl] = kFALSE;
  64     for(Int_t row=1;row<=24;row++)
  65       for(Int_t dil=1;dil<=10;dil++)
  66         for(Int_t adr=0;adr<=47;adr++)
  67           {
  68             SummQ[ddl][adr][dil][row]=0;
  69             SummQ2[ddl][adr][dil][row]=0;
  70             nEntries[ddl][adr][dil][row]=0;
  71           }
  72   }
  73
  74   /* init event counter */
  75   Int_t iEvtNcal=0;
  76
  77   int n;
  78   for (n=1;n<argc;n++) {
  79
  80     status=monitorSetDataSource( argv[n] );
  81     if (status!=0) {
  82       printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
  83       return -1;
  84     }
  85
  86     /* report progress */
  87     /* in this example, indexed on the number of files */
  88     daqDA_progressReport(10+80*n/argc);
  89
  90     for(;;) { // infinite loop
  91       struct eventHeaderStruct *event;
  92       eventTypeType eventT;
  93
  94       /* get next event */
  95       status=monitorGetEventDynamic((void **)&event);
  96       if (status==MON_ERR_EOF) break; /* end of monitoring file has been reached */
  97       if (status!=0) {
  98         printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
  99         return -1;
 100       }
 101
 102       /* retry if got no event */
 103       if (event==NULL) {
 104         break;
 105       }
 106
 107       /* use event - here, just write event id to result file */
 108       eventT=event->eventType;
 109
 110       if (eventT==PHYSICS_EVENT) {
 111
 112         iEvtNcal++;
 113
 114         AliRawReader *reader = new AliRawReaderDate((void*)event);
 115
 116         // Temporary there. Shall be removed as soon as the equipment ID is set correctly
 117         // For the moment ddl.map file contains one line which maps
 118         // the observed eqID=225 to the first HMPID DDL with ID=1536
 119         //      reader->LoadEquipmentIdsMap("ddl.map");
 120
 121         AliHMPIDRawStream stream(reader);
 122
 123         while(stream.Next()) {
 124
 125           Int_t ddl = stream.GetDDLNumber();
 126
 127           for(Int_t row=1; row < 25; row++)
 128             for(Int_t dil=1; dil < 11; dil++)
 129               for(Int_t adr=0; adr < 48; adr++) {
 130                 Short_t q = stream.GetCharge(row,dil,adr);
 131                 if (q<0) continue;
 132
 133                 SummQ[ddl][adr][dil][row] += q;
 134                 SummQ2[ddl][adr][dil][row] += (q*q);
 135                 nEntries[ddl][adr][dil][row]++;
 136                 isDDLOn[ddl] = kTRUE;
 137               }
 138         } //raw data loop
 139
 140         delete reader;
 141
 142       }// if CALIBRATION_EVENT
 143
 144       /* exit when last event received, no need to wait for TERM signal */
 145       if (eventT==END_OF_RUN) {
 146         printf("EOR event detected\n");
 147         break;
 148
 149       } // events loop
 150
 151       free(event);
 152     }
 153
 154   }
 155
 156   /* write report */
 157   printf("Run #%s, received %d calibration events\n",getenv("DATE_RUN_NUMBER"),iEvtNcal);
 158
 159   if (!iEvtNcal) {
 160     printf("No calibration events have been read. Exiting\n");
 161     return -1;
 162   }
 163
 164   /* report progress */
 165   daqDA_progressReport(90);
 166
 167   for(Int_t ddl=0; ddl < 14; ddl++) {
 168     if (!isDDLOn[ddl]) continue;
 169
 170     ofstream out;
 171     out.open(Form("./HmpidPedDdl%02i.txt",ddl));
 172
 173     for(Int_t row=1; row < 25; row++)
 174       for(Int_t dil=1; dil < 11; dil++)
 175         for(Int_t adr=0; adr < 48; adr++) {
 176
 177           Int_t n = nEntries[ddl][adr][dil][row];
 178           if (!n) {
 179             printf("No data for channel: %d %d %d %d\n",ddl,adr,dil,row);
 180             continue;
 181           }
 182           Double_t mean = SummQ[ddl][adr][dil][row]/n;
 183           if ((SummQ2[ddl][adr][dil][row]/n - (SummQ[ddl][adr][dil][row]/n)*(SummQ[ddl][adr][dil][row]/n)) < 0) {
 184             printf("Invalid sums: %d %d %d   %d  %f %f\n",row,dil,adr,n,SummQ[ddl][adr][dil][row],SummQ2[ddl][adr][dil][row]);
 185             return -1;
 186           }
 187           Float_t sigma = TMath::Sqrt(SummQ2[ddl][adr][dil][row]/n - (mean*mean));
 188           Int_t inhard=((Int_t(mean))<<9)+Int_t(mean+3*sigma);
 189
 190           out << Form("%2i %2i %2i %5.2f %5.2f %x\n",row,dil,adr,mean,sigma,inhard);
 191         }
 192
 193     /* store the result file on FES */
 194     status=daqDA_FES_storeFile(Form("./HmpidPedDdl%02i.txt",ddl),Form("HMPID_DA_Pedestals_ddl=%02i",ddl));
 195     if (status) {
 196       printf("Failed to export file : %d\n",status);
 197       return -1;
 198     }
 199   }
 200
 201
 202   /* report progress */
 203   daqDA_progressReport(100);
 204
 205   return status;
 206 }