[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    neventsMin = 2000;

  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

  UChar_t  kNFlagsCut   = 63;// Minimum number of TDC flags requested, replaces the trigger selection

  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 %d %hhu %d %f",
		     &neventsMin,&kStartClock,&kEndClock,&kNPreClocks,&kNPostClocks,&kNFlagsCut,&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("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",
	 neventsMin,kStartClock, kEndClock, kNPreClocks, kNPostClocks, kNFlagsCut, 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 neventsPhysicsAll=0;
  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()) {	
	     neventsPhysicsAll++;
	     UShort_t triggers = rawStream->GetTriggerInputs();
	     UChar_t nFlags = 0;
	     for(Int_t i = 0; i < 64; ++i) {
	       if (rawStream->GetBBFlag(i,10)) nFlags++;
	     }
	     if (nFlags >= kNFlagsCut) { // Check if the number of minimum number of TDC flags is reached (trigger selection)

	       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 (out of %d physics and %d total events)\n",neventsPhysics,neventsPhysicsAll,neventsTotal);
    
//___________________________________________________________________________
//  Computes regression parameters
// charge_i = p0 + charge_tot * p1

  if(neventsPhysics>neventsMin){
    /* 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
4d4eb9a3	63	Int_t neventsMin = 2000;
2ffa87b2	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
468c152f	70	UChar_t kNFlagsCut = 63;// Minimum number of TDC flags requested, replaces the trigger selection
2ffa87b2	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");
4d4eb9a3	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");
2ffa87b2	86	}
	87	fclose(fpConfig);
	88	}
	89
4d4eb9a3	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",
4d4eb9a3	91	neventsMin,kStartClock, kEndClock, kNPreClocks, kNPostClocks, kNFlagsCut, 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 */
4d4eb9a3	122	int neventsPhysicsAll=0;
2ffa87b2	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()) {
4d4eb9a3	168	neventsPhysicsAll++;
2ffa87b2	169	UShort_t triggers = rawStream->GetTriggerInputs();
468c152f	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
2ffa87b2	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
a79dac58	240	printf("%d physics events processed (out of %d physics and %d total events)\n",neventsPhysics,neventsPhysicsAll,neventsTotal);
2ffa87b2	241
	242	//___________________________________________________________________________
	243	// Computes regression parameters
	244	// charge_i = p0 + charge_tot * p1
	245
4d4eb9a3	246	if(neventsPhysics>neventsMin){
6688b0f6	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);
2ffa87b2	257
6688b0f6	258	for(Int_t i=0; i<64; i++) {
	259	fprintf(fp," %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);
	260	printf(" %d %.3f\n",GetOfflineChannel(i), beta[i]*64.);
	261	}
	262
	263	/* close local result file and FXS result file*/
	264	fclose(fp);
2ffa87b2	265	}
6688b0f6	266
2ffa87b2	267	//________________________________________________________________________
2ffa87b2	268
6688b0f6	269	/* export result file to FES */
	270	status=daqDA_FES_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
	271	if (status) {
	272	printf("Failed to export file : %d\n",status);
	273	return -1; }
	274
	275	/* store result file into Online DB */
	276	status=daqDA_DB_storeFile("./V0_EqualizationFactors.dat","V00DAEqualFactors");
	277	if (status) {
	278	printf("Failed to store file into Online DB: %d\n",status);
	279	return -1; }
2ffa87b2	280
	281	return status;
	282	}
	283
	284	Int_t GetOfflineChannel(Int_t channel) {
	285
	286	// Channel mapping Online - Offline:
	287
	288	Int_t fOfflineChannel[64] = {39, 38, 37, 36, 35, 34, 33, 32,
	289	47, 46, 45, 44, 43, 42, 41, 40,
	290	55, 54, 53, 52, 51, 50, 49, 48,
	291	63, 62, 61, 60, 59, 58, 57, 56,
	292	7, 6, 5, 4, 3, 2, 1, 0,
	293	15, 14, 13, 12, 11, 10, 9, 8,
	294	23, 22, 21, 20, 19, 18, 17, 16,
	295	31, 30, 29, 28, 27, 26, 25, 24};
	296	return fOfflineChannel[channel];
	297	}