[u/mrichter/AliRoot.git] / VZERO / VZEROPbPbda.cxx

/*********************************************************************************
- Contact:    Brigitte Cheynis     b.cheynis@ipnl.in2p3.fr
- Link:
- Raw data test file :          
- Reference run number : 137366	      
- Run Type:   PHYSICS
- DA Type:    MON
- Number of events needed: >=2000
- Input Files:  argument list
- Output Files: FXS file     V0_EqualizationFactors.dat (Channel equalization factors)
- Trigger types used: PHYSICS_EVENT
**********************************************************************************/

/**********************************************************************************
*                                                                                 *
* VZERO Detector Algorithm for extracting channel equalization factors for Pb-Pb  *
*                                                                                 *
*                                                                                 *
***********************************************************************************/

// DATE
#include "event.h"
#include "monitor.h"
#include "daqDA.h"

//AliRoot
#include <AliVZERORawStream.h>
#include <AliRawReaderDate.h>
#include <AliRawReader.h>
#include <AliDAQ.h>

// standard
#include <stdio.h>
#include <stdlib.h>

//ROOT
#include "TROOT.h"
#include "TPluginManager.h"
#include <TFile.h>
#include <TH1F.h>
#include <TMath.h>

Int_t GetOfflineChannel(Int_t channel);

/* Main routine --- Arguments: monitoring data source */
      
int main(int argc, char **argv) {

/* magic line from Cvetan */
  gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
                    "*",
                    "TStreamerInfo",
                    "RIO",
                    "TStreamerInfo()");
  int status;
  if (argc!=2) {
     printf("Wrong number of arguments\n");
     return -1;
  }
  
//___________________________________________________
// Get parameters from V00DAEqualFactors.config file

  Int_t    kStartClock = 9;  // First clock in the search for max adc
  Int_t    kEndClock = 11;   // Last clock in the search for max adc
  Int_t    kNPreClocks = 6;  // Number of clock before max used in the charge sum
  Int_t    kNPostClocks = 1; // Number of clock after max used in the charge sum

  UShort_t    kTriggerAcc = 64;    // Trigger mask for accepted events (64 = CTA1 & CTC1)
  UShort_t    kTriggerRej = 256;   // Trigger mask for rejected events (256 = CTA2 & CTC2)

  Int_t    kNBins = 10000;
  Float_t  kRange = 0.1;

  status = daqDA_DB_getFile("V00DAEqualFactors.config","./V00DAEqualFactors.config");
  if (status) {
    printf("Failed to get Config file (V00DAEqualFactors.config) from DAQ DB, status=%d\n", status);
    printf("Take default values of parameters for pedestal calculation \n");
  } else {
    /* open the config file and retrieve cuts */
    FILE *fpConfig = fopen("V00DAEqualFactors.config","r");
    int res = fscanf(fpConfig,"%d %d %d %d %u %u %d %f",
		     &kStartClock,&kEndClock,&kNPreClocks,&kNPostClocks,&kTriggerAcc,&kTriggerRej,&kNBins,&kRange);
    if(res!=8) {
      printf("Failed to get values from Config file (V00DAEqualFactors.config): wrong file format - 7 integers and 1 float are expected - \n");
    }
    fclose(fpConfig);
  }
  
  printf("First LHC Clock = %d; Last LHC Clock = %d; N Pre Clock = %d ; N Post Clock = %d; Trigger mask for accepted events = %u; Trigger mask for rejected events = %u; Number of histogram bins = %d; Histogram range = %.3f\n",
	 kClockStart, kClockStop, kNPreClock, kNPostClock, kTriggerAcc, kTriggerRej, kNBins, kRange);

  TH1D *fMedian[64];
  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);

  Bool_t fFirst = kTRUE;
  Float_t fPrevTotCharge = 0;
  Float_t fPrevadc[64];
  for(Int_t j = 0; j < 64; ++j) fPrevadc[j] = 0;

//___________________________________________________ 

  /* open result file to be exported to FES */
  FILE *fp=NULL;
  fp=fopen("./V0_EqualizationFactors.dat","w");
  if (fp==NULL) {
      printf("Failed to open local result file\n");
      return -1;}
   
  /* define data source : this is argument 1 */  
  status=monitorSetDataSource( argv[1] );
  if (status!=0) {
     printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
     return -1;
  }

  /* declare monitoring program */
  status=monitorDeclareMp( __FILE__ );
  if (status!=0) {
    printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
    return -1;
  }

  /* define wait event timeout - 1s max */
  monitorSetNowait();
  monitorSetNoWaitNetworkTimeout(1000);
  
  /* init counters on events */
  int neventsPhysics=0;
  int neventsTotal=0;


  /* loop on events (infinite) */
  for(;;) {
      struct eventHeaderStruct *event;
      eventTypeType eventT;

      /* check shutdown condition */
      if (daqDA_checkShutdown()) {break;}   

      /* get next event (blocking call until timeout) */
      status=monitorGetEventDynamic((void **)&event);
      if (status==MON_ERR_EOF) {
          printf ("End of File detected\n");
      break; /* end of monitoring file has been reached */
      }
    
      if (status!=0) {
           printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
      break;
      }

      /* retry if got no event */
      if (event==NULL) continue;
               
      /* decode event */
      eventT=event->eventType;
	
      switch (event->eventType){
      
      case START_OF_RUN:
           break;
      
      case END_OF_RUN:
           printf("End Of Run detected\n");
           break;
      
      case PHYSICS_EVENT:
 	     		 
	   AliRawReader *rawReader = new AliRawReaderDate((void*)event);
  
	   AliVZERORawStream* rawStream  = new AliVZERORawStream(rawReader); 
	   if (rawStream->Next()) {	
	     UShort_t triggers = rawStream->GetTriggerInputs();
	     if (((triggers & kTriggerAcc) == kTriggerAcc) &&  // Check if the requested trigger(s) is fired
		 ((triggers & kTriggerRej) == 0)) { // Check if the requested trigger(s) is NOT fired
	       neventsPhysics++;

	       Float_t adc[64];
	       Float_t totCharge = 0.;
	       for(Int_t i = 0; i < 64; ++i) {
		 adc[i] = 0.;
		 Float_t maxadc=0.;
		 Int_t imax=-1;
		 for(Int_t j = kStartClock; j <= kEndClock; ++j) {
		   Float_t charge = (Float_t)(rawStream->GetPedestal(i,j));
		   if(charge > maxadc){
		     maxadc = charge;
		     imax   = j;
		   }
		 }

		 if (imax != -1) {
		   Int_t start = imax - kNPreClocks;
		   if (start < 0) start = 0;
		   Int_t end = imax + kNPostClocks;
		   if (end > 20) end = 20;
		   for(Int_t iClock = start; iClock <= end; iClock++) {
		     adc[i] += (Float_t)(rawStream->GetPedestal(i,iClock));
		   }
		 }
		 totCharge += adc[i];
	       }

	       if (fFirst) {
		 fFirst = kFALSE;
		 fPrevTotCharge = totCharge;
		 for(int i = 0; i < 64; ++i) fPrevadc[i] = adc[i];
	       }
	       else {
		 fFirst = kTRUE;
		 Float_t deltaTotCharge = totCharge - fPrevTotCharge;
		 Float_t weight = deltaTotCharge*deltaTotCharge;
		 if (weight > 1) {
		   for(int i = 0; i < 64; ++i) {
		     fMedian[i]->Fill((adc[i]-fPrevadc[i])/deltaTotCharge,weight);
		   }
		 }
	       }

	     }
	   }   // End : if rawstream
           delete rawStream;
           rawStream = 0x0;      
           delete rawReader;
           rawReader = 0x0;	     	 						         
      } // end of switch on event type 
	
      neventsTotal++;
      /* free resources */
      free(event);
	
      /* exit when last event received, no need to wait for TERM signal */
      if (eventT==END_OF_RUN) {
	printf("End Of Run event detected\n");
	break;
      }

  }  // loop over events
  
  printf("%d physics events processed\n",neventsPhysics);
    
//___________________________________________________________________________
//  Computes regression parameters
// charge_i = p0 + charge_tot * p1

  Double_t beta[64];
  Double_t q = 0.5;
  for(int i = 0; i < 64; ++i) fMedian[i]->GetQuantiles(1,&beta[i],&q);

  for(Int_t i=0; i<64; i++) {
    fprintf(fp," %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);				       
    printf(" %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);				       
  }
 
//________________________________________________________________________
   
  /* close local result file and FXS result file*/
  fclose(fp);

  if(neventsPhysics>2000){
    /* export result file to FES */
    status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
    if (status)    {
      printf("Failed to export file : %d\n",status);
      return -1; }

    /* store result file into Online DB */
    status=daqDA_DB_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
    if (status)    {
      printf("Failed to store file into Online DB: %d\n",status);
      return -1; }
  }

  return status;
}

 Int_t GetOfflineChannel(Int_t channel) {

// Channel mapping Online - Offline:
 
 Int_t fOfflineChannel[64] = {39, 38, 37, 36, 35, 34, 33, 32, 
                              47, 46, 45, 44, 43, 42, 41, 40, 
			      55, 54, 53, 52, 51, 50, 49, 48, 
			      63, 62, 61, 60, 59, 58, 57, 56,
			       7,  6,  5,  4,  3,  2,  1,  0, 
			      15, 14, 13, 12, 11, 10,  9,  8,
			      23, 22, 21, 20, 19, 18, 17, 16, 
			      31, 30, 29, 28, 27, 26, 25, 24};
 return	fOfflineChannel[channel];		      
}
Commit	Line	Data
2ffa87b2	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
	64	Int_t kStartClock = 9; // First clock in the search for max adc
65	Int_t kEndClock = 11; // Last clock in the search for max adc
66	Int_t kNPreClocks = 6; // Number of clock before max used in the charge sum
67	Int_t kNPostClocks = 1; // Number of clock after max used in the charge sum
68
69	UShort_t kTriggerAcc = 64; // Trigger mask for accepted events (64 = CTA1 & CTC1)
70	UShort_t kTriggerRej = 256; // Trigger mask for rejected events (256 = CTA2 & CTC2)
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 %u %u %d %f",
83176645	83	&kStartClock,&kEndClock,&kNPreClocks,&kNPostClocks,&kTriggerAcc,&kTriggerRej,&kNBins,&kRange);
2ffa87b2	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("First LHC Clock = %d; Last LHC Clock = %d; N Pre Clock = %d ; N Post Clock = %d; Trigger mask for accepted events = %u; Trigger mask for rejected events = %u; Number of histogram bins = %d; Histogram range = %.3f\n",
	91	kClockStart, kClockStop, kNPreClock, kNPostClock, kTriggerAcc, kTriggerRej, 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
83176645	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
2ffa87b2	101	//___________________________________________________
	102
	103	/* open result file to be exported to FES */
	104	FILE *fp=NULL;
	105	fp=fopen("./V0_EqualizationFactors.dat","w");
	106	if (fp==NULL) {
	107	printf("Failed to open local result file\n");
	108	return -1;}
	109
	110	/* define data source : this is argument 1 */
	111	status=monitorSetDataSource( argv[1] );
	112	if (status!=0) {
	113	printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
	114	return -1;
	115	}
	116
	117	/* declare monitoring program */
	118	status=monitorDeclareMp( __FILE__ );
	119	if (status!=0) {
	120	printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
	121	return -1;
	122	}
	123
	124	/* define wait event timeout - 1s max */
	125	monitorSetNowait();
	126	monitorSetNoWaitNetworkTimeout(1000);
	127
	128	/* init counters on events */
	129	int neventsPhysics=0;
	130	int neventsTotal=0;
	131
	132
	133	/* loop on events (infinite) */
	134	for(;;) {
	135	struct eventHeaderStruct *event;
	136	eventTypeType eventT;
	137
	138	/* check shutdown condition */
	139	if (daqDA_checkShutdown()) {break;}
	140
	141	/* get next event (blocking call until timeout) */
	142	status=monitorGetEventDynamic((void **)&event);
	143	if (status==MON_ERR_EOF) {
	144	printf ("End of File detected\n");
	145	break; /* end of monitoring file has been reached */
	146	}
	147
	148	if (status!=0) {
	149	printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
	150	break;
	151	}
	152
	153	/* retry if got no event */
	154	if (event==NULL) continue;
	155
	156	/* decode event */
	157	eventT=event->eventType;
	158
	159	switch (event->eventType){
	160
	161	case START_OF_RUN:
	162	break;
	163
	164	case END_OF_RUN:
165	printf("End Of Run detected\n");
166	break;
167
168	case PHYSICS_EVENT:
169
170	AliRawReader rawReader = new AliRawReaderDate((void)event);
171
172	AliVZERORawStream* rawStream = new AliVZERORawStream(rawReader);
173	if (rawStream->Next()) {
174	UShort_t triggers = rawStream->GetTriggerInputs();
175	if (((triggers & kTriggerAcc) == kTriggerAcc) && // Check if the requested trigger(s) is fired
176	((triggers & kTriggerRej) == 0)) { // Check if the requested trigger(s) is NOT fired
177	neventsPhysics++;
178
179	Float_t adc[64];
180	Float_t totCharge = 0.;
181	for(Int_t i = 0; i < 64; ++i) {
182	adc[i] = 0.;
183	Float_t maxadc=0.;
184	Int_t imax=-1;
185	for(Int_t j = kStartClock; j <= kEndClock; ++j) {
186	Float_t charge = (Float_t)(rawStream->GetPedestal(i,j));
187	if(charge > maxadc){
188	maxadc = charge;
189	imax = j;
190	}
191	}
192
193	if (imax != -1) {
194	Int_t start = imax - kNPreClocks;
195	if (start < 0) start = 0;
196	Int_t end = imax + kNPostClocks;
197	if (end > 20) end = 20;
198	for(Int_t iClock = start; iClock <= end; iClock++) {
199	adc[i] += (Float_t)(rawStream->GetPedestal(i,iClock));
200	}
201	}
202	totCharge += adc[i];
203	}
204
205	if (fFirst) {
206	fFirst = kFALSE;
207	fPrevTotCharge = totCharge;
208	for(int i = 0; i < 64; ++i) fPrevadc[i] = adc[i];
209	}
210	else {
211	fFirst = kTRUE;
212	Float_t deltaTotCharge = totCharge - fPrevTotCharge;
213	Float_t weight = deltaTotCharge*deltaTotCharge;
214	if (weight > 1) {
215	for(int i = 0; i < 64; ++i) {
216	fMedian[i]->Fill((adc[i]-fPrevadc[i])/deltaTotCharge,weight);
217	}
218	}
219	}
220
221	}
222	} // End : if rawstream
223	delete rawStream;
224	rawStream = 0x0;
225	delete rawReader;
226	rawReader = 0x0;
227	} // end of switch on event type
228
229	neventsTotal++;
230	/* free resources */
231	free(event);
232
233	/* exit when last event received, no need to wait for TERM signal */
234	if (eventT==END_OF_RUN) {
235	printf("End Of Run event detected\n");
236	break;
237	}
238
239	} // loop over events
240
241	printf("%d physics events processed\n",neventsPhysics);
242
243	//___________________________________________________________________________
244	// Computes regression parameters
245	// charge_i = p0 + charge_tot * p1
246
247	Double_t beta[64];
248	Double_t q = 0.5;
249	for(int i = 0; i < 64; ++i) fMedian[i]->GetQuantiles(1,&beta[i],&q);
250
251	for(Int_t i=0; i<64; i++) {
252	fprintf(fp," %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);
253	printf(" %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);
254	}
255
256	//________________________________________________________________________
257
258	/* close local result file and FXS result file*/
259	fclose(fp);
260
261	if(neventsPhysics>2000){
262	/* export result file to FES */
263	status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
264	if (status) {
265	printf("Failed to export file : %d\n",status);
266	return -1; }
267
268	/* store result file into Online DB */
269	status=daqDA_DB_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
270	if (status) {
271	printf("Failed to store file into Online DB: %d\n",status);
272	return -1; }
273	}
274
275	return status;
276	}
277
278	Int_t GetOfflineChannel(Int_t channel) {
279
280	// Channel mapping Online - Offline:
281
282	Int_t fOfflineChannel[64] = {39, 38, 37, 36, 35, 34, 33, 32,
283	47, 46, 45, 44, 43, 42, 41, 40,
284	55, 54, 53, 52, 51, 50, 49, 48,
285	63, 62, 61, 60, 59, 58, 57, 56,
286	7, 6, 5, 4, 3, 2, 1, 0,
287	15, 14, 13, 12, 11, 10, 9, 8,
288	23, 22, 21, 20, 19, 18, 17, 16,
289	31, 30, 29, 28, 27, 26, 25, 24};
290	return fOfflineChannel[channel];
291	}