[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",
	 kStartClock, kEndClock, kNPreClocks, kNPostClocks, 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;

//___________________________________________________ 
   
  /* 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

  if(neventsPhysics>2000){
    /* 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;}

    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);
  }

//________________________________________________________________________
   
  /* 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",
8e2ea3c5	91	kStartClock, kEndClock, kNPreClocks, kNPostClocks, kTriggerAcc, kTriggerRej, kNBins, kRange);
2ffa87b2	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	//___________________________________________________
2ffa87b2	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 neventsPhysics=0;
	123	int neventsTotal=0;
	124
	125
	126	/* loop on events (infinite) */
	127	for(;;) {
	128	struct eventHeaderStruct *event;
	129	eventTypeType eventT;
	130
	131	/* check shutdown condition */
	132	if (daqDA_checkShutdown()) {break;}
	133
	134	/* get next event (blocking call until timeout) */
	135	status=monitorGetEventDynamic((void **)&event);
	136	if (status==MON_ERR_EOF) {
	137	printf ("End of File detected\n");
	138	break; /* end of monitoring file has been reached */
	139	}
	140
	141	if (status!=0) {
	142	printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
	143	break;
	144	}
	145
	146	/* retry if got no event */
	147	if (event==NULL) continue;
	148
	149	/* decode event */
	150	eventT=event->eventType;
	151
	152	switch (event->eventType){
	153
	154	case START_OF_RUN:
	155	break;
	156
	157	case END_OF_RUN:
	158	printf("End Of Run detected\n");
	159	break;
	160
	161	case PHYSICS_EVENT:
	162
	163	AliRawReader rawReader = new AliRawReaderDate((void)event);
	164
	165	AliVZERORawStream* rawStream = new AliVZERORawStream(rawReader);
166	if (rawStream->Next()) {
167	UShort_t triggers = rawStream->GetTriggerInputs();
168	if (((triggers & kTriggerAcc) == kTriggerAcc) && // Check if the requested trigger(s) is fired
169	((triggers & kTriggerRej) == 0)) { // Check if the requested trigger(s) is NOT fired
170	neventsPhysics++;
171
172	Float_t adc[64];
173	Float_t totCharge = 0.;
174	for(Int_t i = 0; i < 64; ++i) {
175	adc[i] = 0.;
176	Float_t maxadc=0.;
177	Int_t imax=-1;
178	for(Int_t j = kStartClock; j <= kEndClock; ++j) {
179	Float_t charge = (Float_t)(rawStream->GetPedestal(i,j));
180	if(charge > maxadc){
181	maxadc = charge;
182	imax = j;
183	}
184	}
185
186	if (imax != -1) {
187	Int_t start = imax - kNPreClocks;
188	if (start < 0) start = 0;
189	Int_t end = imax + kNPostClocks;
190	if (end > 20) end = 20;
191	for(Int_t iClock = start; iClock <= end; iClock++) {
192	adc[i] += (Float_t)(rawStream->GetPedestal(i,iClock));
193	}
194	}
195	totCharge += adc[i];
196	}
197
198	if (fFirst) {
199	fFirst = kFALSE;
200	fPrevTotCharge = totCharge;
201	for(int i = 0; i < 64; ++i) fPrevadc[i] = adc[i];
202	}
203	else {
204	fFirst = kTRUE;
205	Float_t deltaTotCharge = totCharge - fPrevTotCharge;
206	Float_t weight = deltaTotCharge*deltaTotCharge;
207	if (weight > 1) {
208	for(int i = 0; i < 64; ++i) {
209	fMedian[i]->Fill((adc[i]-fPrevadc[i])/deltaTotCharge,weight);
210	}
211	}
212	}
213
214	}
215	} // End : if rawstream
216	delete rawStream;
217	rawStream = 0x0;
218	delete rawReader;
219	rawReader = 0x0;
220	} // end of switch on event type
221
222	neventsTotal++;
223	/* free resources */
224	free(event);
225
226	/* exit when last event received, no need to wait for TERM signal */
227	if (eventT==END_OF_RUN) {
228	printf("End Of Run event detected\n");
229	break;
230	}
231
232	} // loop over events
233
234	printf("%d physics events processed\n",neventsPhysics);
235
236	//___________________________________________________________________________
237	// Computes regression parameters
238	// charge_i = p0 + charge_tot * p1
239
6688b0f6	240	if(neventsPhysics>2000){
	241	/* open result file to be exported to FES */
	242	FILE *fp=NULL;
	243	fp=fopen("./V0_EqualizationFactors.dat","w");
	244	if (fp==NULL) {
	245	printf("Failed to open local result file\n");
	246	return -1;}
	247
	248	Double_t beta[64];
	249	Double_t q = 0.5;
	250	for(int i = 0; i < 64; ++i) fMedian[i]->GetQuantiles(1,&beta[i],&q);
2ffa87b2	251
6688b0f6	252	for(Int_t i=0; i<64; i++) {
	253	fprintf(fp," %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);
	254	printf(" %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);
	255	}
	256
	257	/* close local result file and FXS result file*/
	258	fclose(fp);
2ffa87b2	259	}
6688b0f6	260
2ffa87b2	261	//________________________________________________________________________
2ffa87b2	262
6688b0f6	263	/* export result file to FES */
	264	status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
	265	if (status) {
	266	printf("Failed to export file : %d\n",status);
	267	return -1; }
	268
	269	/* store result file into Online DB */
	270	status=daqDA_DB_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
	271	if (status) {
	272	printf("Failed to store file into Online DB: %d\n",status);
	273	return -1; }
2ffa87b2	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	}