Sana Cherif

Chiraz Trabelsi

Samy Meftali

Jean-Luc Dekeyser

HIGH LEVEL DESIGN OF ADAPTIVE DISTRIBUTED CONTROLLER FOR PARTIAL DYNAMIC RECONFIGURATION IN FPGA

Keywords: Partial Dynamic Reconfiguration, distributed control, modularity, adaptivity, high level modeling, UML MARTE

Controlling dynamic and partial reconfigurations becomes one of the most important key issues in modern embedded systems design. In fact, in such systems, the reconfiguration controller can significantly affect the system performances. Indeed, the controller has to handle efficiently three major tasks during runtime: observation (monitoring), taking reconfiguration decisions and notify decisions to the rest of the system in order to realize it. We present in this paper a novel high level approach permitting to model, using MARTE UML profile, modular and flexible distributed controllers for dynamic reconfiguration management. This approach permits components/ models reuse and allows systematic code generation. It consequently makes reconfigurable systems design less tedious and reduces time to market.

INTRODUCTION

In a reconfigurable computing system, both the hardware and software parts can be reconfigured depending upon the designer requirements. Due to their several advantages, such as flexibility and performance, reconfigurable devices, mainly FPGAs, make them well adapted for diverse application domains such as satellite based systems, telecommunications, transport, etc. As the computational power increases, more functionalities are expected to be integrated into the system. Today's FPGAs are more suitable for evolving systems with Partial Dynamic Reconfiguration (PDR) feature. It enables modification of specific regions of an FPGA on the fly while others regions remain operational and functioning. Although PDR is fastly growing and gaining popularity in the domain of reconfigurable computing, it is mandatory to predict a robust control mechanism for managing reconfiguration. A good control system should have these prerequisites : speed, accuracy and robustness.

In our work, the reconfiguration decision is taken in a distributed way, which is better in terms of performance than the centralized way. In order to decrease the design complexity of our system, we will use a model-driven design method-ology. UML (Unified Modeling Language) is considered as one of the main unified visual modeling languages in Model Driven Engineering (MDE). When UML is utilized in the framework of MDE, it is usually carried out through by its internal metamodeling mechanism termed as a profile. UML MARTE [START_REF] Omg | Modeling and analysis of real-time and embedded systems (MARTE)[END_REF](Modeling and Analysis of Real-Time and Embedded Systems) is a recent industry standard UML profile of OMG, dedicated to model-driven development of embedded systems. MARTE extends UML along with added extensions (e.g. performance and scheduling analysis). The MARTE profile extends the possibilities to model the features of software and hardware parts of a real-time embedded system and their relations. MARTE is mainly organized into two packages, the MARTE design model and the MARTE analysis model which share common concepts, grouped in the MARTE foundations package: for expressing non-functional properties (NFPs),timing notions (Time), resource modeling (GRM), components (GCM) and allocation concepts (Alloc). An additional package contains the annexes defined in the MARTE profile along with predefined libraries. The hardware concepts in MARTE are grouped in the Hardware Resource Model (HRM) package from MARTE design package [START_REF] Taha | An open framework for detailed hardware modeling[END_REF]. HRM consists of a functional and a physical views. The first one (HwLogical) classifies hardware resources depending on their functional properties, and the second view (Hw-Physical) concentrates on their physical nature. Both derive certain concepts from the HwGeneral root package in which HwResource is a core concept that denotes a generic hardware entity providing at least one HwResourceService, and may require some services from other resources. The functional view of HRM defines hardware resources as either computing, storage, communication, timing or device resources. The physical view represents physical components with details about their shape, size and power consumption among other attributes. In this paper, we propose a high level design of an adaptive distributed controller for partial dynamic reconfiguration, using UML MARTE profile. The rest of the paper is as follows. Section 2 gives an overview of some related works while section 3 explains our motivations and our proposed approach. We validate our solution with a case study in section 4 and presents several advantages in section 5. Finally, the last section presents a conclusion and future works.

RELATED WORKS

There are a few works that have carried out research using the MARTE profile for specification of reconfigurable systems. The MoPCoM project [START_REF] Koudri | Using marte in the mopcom soc/sopc co-methodology[END_REF] aims to target modeling and code generation of dynamically reconfigurable embedded systems using the MARTE UML profile for SoC Co-Design [START_REF] Vidal | A co-design approach for embedded system modeling and code generation with uml and marte[END_REF]. However, the targeted applications are extremely simplistic in nature, and do not represent complex application domains normally targeted in the SoC industry. The authors in the Over-SoC project [START_REF] Pillement | High-level model of dynamically reconfigurable architectures[END_REF] also provide a high level modeling methodology for implementing dynamic reconfigurable architectures. They integrate an operating system for providing and handling the reconfiguration mechanism.

Robust control mechanisms for managing reconfiguration must be developed in order to take into account SoC designer needs as well as QoS choices. In [START_REF] Missal | Hierarchical distributed controllers design and verification[END_REF], authors reuse the design methodology based on Automation Objects presented in [START_REF] Karras | Rapid engineering and reconfiguration of automation objects aided by formal modelling and verification[END_REF] to address formal approaches to verify the control architecture in a hierarchical way. Therefore, they propose an hierarchical control approach by means of a multi-layered architecture. Authors in [START_REF] Dekeyser | Heterogeneous Embedded Systems -Design Theory and Practice, chapter High-level modeling of dynamically reconfigurable heterogeneous systems[END_REF] present control semantics at the deployment level where each configuration can be viewed as an implementation of the modeled application: comprised of unique collection of Intellectual Properties (IPs). Each IP is related to an application elementary component. The control at deployment level can be extended for complex configurations. A comparison of several design approaches of distributed controllers in industrial automation systems is presented in [START_REF] Vyatkin | Systematic Design and Implementation of Distributed Controllers in Industrial Automation[END_REF] . Authors discussed essentially two different approaches: decomposition of a centralized controller onto several communicating distributed controllers and integration of predefined controllers of components to the control of a system. In [START_REF] Kumpati | Adaptive control using multiple models, switching, and tuning[END_REF], authors present a Multiple Models Switching and Tuning approach(MMST) for adaptive control to provide greater flexibility in design. They also consider deterministic and stochastic systems and discuss some problems of convergence and stability.

PROPOSED APPROACH

Motivations

Partial Dynamic Reconfiguration (PDR) gives the designer the ability to reconfigure a certain portion of the FPGA during run-time without influencing the other parts. This special feature allows the hardware to be adaptable to any potential situation especially when a high degree of flexibility and parallelism are required to achieve real-time constraints. Normally, the reconfiguration is carried out in complex systems by means of a control mechanism (for example a RTOS or middleware), which can depend on QoS choices: such as changes in executing functionalities due to designer requirements, or changes due to resource constraints of targeted hardware/platforms. The controller can be implemented in different manners. It can be external or internal; and can be an hardware module written in an HDL acting as a controller. It also can be integrated in different levels in a SoC Co-Design framework [START_REF] Dekeyser | Heterogeneous Embedded Systems -Design Theory and Practice, chapter High-level modeling of dynamically reconfigurable heterogeneous systems[END_REF]. Control integration at the application level has a local impact and is independent of the architecture or the allocation. The reconfiguration at this level may arise due to QoS criteria as in video processing application where switching from a high resolution mode to a lower one is required. Control integration in an architecture can be mainly local and used for QoS such as modification of the hardware parameters (voltage or frequency) for manipulating power consumption levels. The second type of control can be used to modify the system structure either globally or partially. If either the application or the architecture is changed, the allocation must be adapted accordingly. Nevertheless, dynamic reconfiguration is more architecture driven as compared to being application oriented, and is influenced by target platform low-level details. Hence, we are interested in control integration at the architecture level. Dynamic adaptation depends on the context required by designer and can be determined by different QoS criteria [START_REF] Bruneton | The fractal component model and its support in java: Experiences with auto-adaptive and reconfigurable systems[END_REF]. Self adaptation is also present in some component models. Usually this type of adaptation is executed by an adaptation manager or controller, typically a separate process or component that is responsible for managing the adaptations. The aim of our proposed approach is to model at a high abstraction level, a distributed modular controller which will carry out some of QoS factors such as changes in executing functionalities due to designer requirements, or changes due to resource constraints of targeted hardware/platforms. The changes can also take place due to other environmental criteria such as communication quality, time and area consumed for reconfiguration and energy consumption.

Controller modeling

Our proposed controller consists of a distributed and local device which manages each partial reconfigurable region (PRR) in an FPGA. While as the use of a master controller is unavoidable, a master controller (processor) is mandatory for overall system level control. The detailed block diagram of controller architecture is shown in Figure 1.

Our Controller characterized by its modularity and adaptability, will be based on the following notions:

• Ctrl2Proc Adapter: adapts signals transmitted from the controller to the processor via bus macro.

• Event Observer: has as input the processor's requests, e.g., read, write, etc. or user inputs sent from the hyper-Terminal. It treats these data and transforms them into a list of events.

• Decision module: analyzes received events coming from Observer module and decides which configuration can be run. This choice depends on defined • Notifier : this module forwards the configuration decision to the ICAP (Internal Configuration Port Access) which will load the chosen configuration.

• Ctrl2PRR Adapter :adapts signals between our modular controller and its corresponding reconfigurable region.

Our proposed controller is highly flexible thanks to its adaptivity and modularity. In fact, the Event Observer and the Decision module are adaptable and each one can contain different functionalities according to observed data and events,respectively. Figure 2 shows a high level model of our proposed modular controller using MARTE profile concepts. Each component contains the MARTE HRM HwComponent stereotype to indicate that these modules are physical in nature.

As shown in Figure 2, we use a softcore Microblaze processor as example of a master controller which will monitor all the distributed controllers. It has MARTE HRM HwProcessor stereotype and is connected to the controllerToProcessor Adapter via the Processor Local Bus (PLBbus). The Partial reconfigurable region (PRR) typed as the generic HwResource type in order to illustrate that the partially reconfigurable region can be either generic or have a specific functionality. This PRR consists of a reconfigurable Hardware Accelerator (HwAcc) defined as HwPLD, as it is reconfigurable compared to a typical hardware accelerator in a large-scale SoC design which can be seen as a HwASIC (after fabrication).

CASE STUDY

We present here a case study of a complete FPGA model to illustrate our modeling methodology.

Application level

The H.264 is currently one of the most commonly used formats for the recording, compression, and distribution of high definition video [START_REF] Bhaskaran | Image and Video Compression Standards -Algorithms and Architectures[END_REF]. This standard has already been adopted for a wide range of applications, from low bit-rate mobile video to high-definition TV. Like other previous video coding standards, the H.264 utilizes video coding algorithms based on block-based motion compensation and transform based spatial coding framework. The H.264 video encoder treats the frames as I (Integral frame) or as P (Partial frame) [START_REF] Ivanov | Real-time h.264 video encoding in software with fast mode decision and dynamic complexity control[END_REF]. The I frames are encoded without taking into consideration the previous frame, while the P frames are treated taking into account the information existing in the previous frame. Figure 3 describes the encoding principle for a PFrame. Just like encoding an I frame, the algorithm partitions the frame into 16x16 pixels MacroBlocks (MB). Motion Estimation (ME) compresses the temporal redundancies between consecutive frames [START_REF] Burleson | Dynamically parameterized architectures for power-aware video coding: Motion estimation and dct[END_REF]. The algorithm compares a Current MB (CurrMB)(16x16 pixels) from the current frame to all the candidate blocks in a search window in the previous frame and the best match is used to estimate the motion. For each CurrMB, the ME module can obtain two possible results. The first possibility is that this MB is similar to the previous one called referential MB (RefMB) (with possible presence of small differences). The second possibility is that they are completely different, i.e., the current MB contains 16x16 new pixels. In the case of similarity, the ME will provide a Motion Vector (MVect) for the CurrMB to indicate its origin from the previous decoded frame. Then, the Motion Compensation module (MC) takes as inputs this MVect and the RefMB. Then, using temporal prediction, it provides a Predicted MB (PredMB), which will store the full information needed to transform the previous frame into the next frame. Nonetheless, the similarity between CurrMB and PredMB is not absolute. This leads to use Subtractor module (Subtract) to compute differences between them in order to obtain a MB that contains only residues (extra data) called Residual MB (ResMB). Therefore, after the ME, MC and difference phases, the current MB which initially contains a 16x16 original pixels, will contain only some few data. The rest can already be found somewhere in the previous decoded frame. By using the Discrete Cosine Transform (DCT), the frame values stored into a MacroBlock will be transformed into frequency values. The DCT video coding, particularly the 2D DCT, can be decomposed as two 1D DCT computations together with a transpose memory in order to reduce the computational complexity. The DCT block implements a row-column Distributed Arithmetic (DA) [START_REF] Bhaskaran | Image and Video Compression Standards -Algorithms and Architectures[END_REF] version of the Chen fast [START_REF] Smith | A fast computational algorithm for the discrete cosine transform[END_REF] DCT algorithm. Furthermore, 1-D DCT is first performed row by row on the input data and the results are saved in a transpose memory. Then, 1-D DCT is performed column by column on the results stored in the transpose memory. The outputs obtained from the second 1-D DCT are the coefficients of the 2-D DCT, i.e., a Coefficient Matrix (Co-effMX). Using the transformed MB after the DCT function, the next step is to quantize (Quant) these values by removing a specific number of bits from all the 8x8 DCT values. This number can vary from frame to frame but is unique during the quantization of the MBs from the same frame. As a result, we obtain a quantization matrix (QuantMx). After that, an inverse transform (inverse quantization (IQuant) and inverse DCT (IDCT)) is applied to this matrix to generate a reconstructed residual MacroBlock (R.ResMB). This one in sent to an adder module (Add) and finally stored in the SRAM memory. The last step is to code quantized values using Entropy Coding (EC) techniques. Entropy coding is a term referring to a lossless data compression scheme that replace data elements with coded representations which, in combination with the previously-described predictions, transformations, and quantization, can result in significantly reduced data size [START_REF] Sullivan | The h.264/avc advanced video coding standard:overview and introduction to the fidelity range extensions[END_REF]. In H.264, two modes of EC designs are used in a context adaptive (CA) way: CA Variable Length Coding CAVLC and CA Binary Arithmetic Coding CABAC. After processing all the MBs of the frame using the P frame coding, the final result that we obtain is a bitstream. This application is modeled with UML MARTE concepts

System architecture

In Figure 5, we present the global structure of our reconfigurable architecture, implemented on the Xilinx Virtex 4 XC4VFX20. Our system can be mainly divided into two regions. The static region and the dynamically reconfigurable one. The static region mainly consists of a softcore Microblaze processor which is selected as the global reconfiguration controller and other necessary peripherals for dynamic reconfiguration. The Microblaze processor connects directly to the high speed 64-bit Processor Local Bus (PLB) and links with On-Chip Peripheral Bus (OPB) via a PLB to OPB Bridge (plb2opbBridge). Several slave peripherals are connected to the OPB bus, such as SystemACE controller for accessing the partial bitstreams placed in an external on board Com-pact Flash (CF) card, an ICAP core (Xilinx OPB HwIcap module) which carries out partial reconfiguration using the read-modify-write mechanism. An UART controller module is connected to a hyperTerminal from which the user inputs/events are sent to the processor. Finally, four processors generated from the modeled application are connected to PLB bus, i.e.,quant pr (quantization), cod pr (coding), sub pr (subtractor) and MC pr (Motion Compensation). In video coding application, ME and DCT are the most compute intensive processes. Inspired by proposed architecture in [START_REF] Huang | Scalable FPGA-based Architecture for DCT Computation Using Dynamic Partial Reconfiguration[END_REF], dynamically reconfigurable region contains eight hardware accelerators (HwAcc) connected to PLB bus via bus macros. Each one serves as a partial reconfigurable region (PRR) which can be dynamically reconfigured using partial reconfiguration. To achieve quality scalability of DCT computations, hardware Accelerators can be reconfigured in two ways. They can be removed or added to achieve zonal coding for the DCT coefficients (from 1x1 to 8x8), or we can simply change the internal logic of the HwAcc through dynamic partial reconfiguration to reduce the precision of resultant DCT coefficients. Since, ME is a highly parallel application, it should be implemented on a Hardware accelerator (HwAcc). Therefore, we can use ME computation as an example to demonstrate how to take advantage of the unused HwAccs.

Figure 6 finally shows our reconfigurable architecture (Xilinx Virtex 4 XC4VFX20 chip) using MARTE concepts in a merged functional/physical view to express all the necessary attributes related to the corresponding physical/logical stereotypes. We use 4x4 DCT with 4 HwAccs for ME as an example for modeling with MARTE profile.

ADVANTAGES OF OUR SOLUTION

Our proposed solution combines mainly three aspects : distribution, adaptability and modularity. In this section, we illustrate the advantages of our approach with some examples.

Spatial distribution:

-Performance: In modern embedded systems, the reconfiguration controller can significantly affect the system performance. When the performance of a centralized controller is not sufficient, often it is a good idea to divide the control program to several pieces and run them in different devices. Indeed, using only one global controller to manage several partial reconfigurable regions can lead to deadlocks and consequently performance degradation.

-Complexity: Partial dynamic reconfiguration becomes one of the most key issues in today's FPGAs. These systems are becoming more agile to meet the requirements of ever changing markets. Due to the presence of several partial dynamic regions (PRR) in an FPGA, it is mandatory to efficiently manage them. Therefore, implementing a distributed controller to manage each PRR enhances flexible reconfiguration and offers reusability of the controllers. Such a splitting may bring an immediate gain in the productivity and the complexity of a global reconfiguration management is reduced significantly.

This solution makes the control design less complex and increases the performance, compared to a centralized solution.

Modularity and adaptivity:

-Reconfiguration Time:

Three main functionalities are handled by our modular controller: observation, taking reconfiguration decision and notification. The main objective for proposing a modular controller is to ease the reconfiguration of some sub-modules of our distributed controller without needing to change the others. Since an unused HwAcc can be utilized by dynamic partial reconfiguration for other functions, Motion Estimation and DCT can be performed concurrently. In this case, this accelerator controller will not totally change. Furthermore, the kind of data observed by Event Observer module is a pixel, i.e., MacroBlock. So that only the decision functionality changes and adapts to the ME's constraints. Thus, in this case, we will gain in terms of time reconfiguration.

-Energy consumption gain: The Event Observer module can detect the current activity of the processing element (hardware accelerator) and the related power consumption. Thereafter, it sends this information to the decision module. If the total consumption exceed a certain threshold, the monitor can decide that a reconfiguration is needed. So, it sends the decision to the processor responsible for the reconfiguration launching. A reconfiguration can, for example, load another version of the same algorithm but with a lower performance. For an image processing application this would be done by loading, for example, a filter that is less precise than the previous one. Moreover, in order to process a High Definition frame (Event Observer module's input is a HD frame), we are dealing then with intensive computing thus increasing parallelism degree. Therefore, we need to <<hwEndPoint>> CFlash <<hwEndPoint>> hyperteminal Fig. 6. Modeling of our reconfigurable architecture increase the number of processing elements. This leads to a higher performance. However, when there is not enough energy in the battery to carry on with the same parallelism degree, decreasing the number of the operational processing elements is required. This reconfiguration consists of switching from a high resolution mode to a lower one (Standard Definition frame). In this case, such information (battery level) is transmitted from the Event Observer module to the decision module which decides to turn off its corresponding accelerator in order to not consume clock cycles. Obviously a controller of a partial reconfigurable region (PRR) has to collaborate and communicate with other PRR controllers (as any local controller). As in many real situations, the use of a master controller is unavoidable. In our architecture, the Microblaze is considered as the master controller. While the control mechanism depends on QoS choices, the Microblaze can provide changes due to resource constraints of targeted hardware or configuration constraints to local controllers in order to turn on/off processing elements for each reconfigurable region. In the design model, these changes are modeled using Component state attribute from MARTE HRM HwComponent stereotype as shown in Figure 6. A Component state of a partial dynamic region (PRR) set to Storage instead of Operating means that this PRR is disabled or inactive and does not consume clock cycles. The changes can also depend on other environmental criteria such as communication quality, time and area consumed for reconfiguration; and energy consumption levels. In our proposed solution, reduction of accelerators number (e.g.,DCT dimension from 8x8 to 1x1) can reduce overall power consumption, total number of clocks cycles required for the running module. Thus, our proposed solution can be beneficial in terms of hardware resources and power consumption. The FPGA does not need to be stopped while changing the configuration which is important for many image/video applications. Finally, an efficient reconfiguration would find a trade-off between performance and power.

CONCLUSION AND FUTURE WORKS

Dynamic and partial reconfiguration is supported by more and more recent FPGA platforms. Controlling reconfigurations in such systems becomes one of the most difficult and important issues, because of its impact on systems performances. We presented in this paper a new high level approach for modeling and designing dynamic reconfiguration controllers, starting from UML MARTE profile. We validated our design approach with a case study which consisting of a P Frame (Partial Frame) H.264 video encoding. Finally, we illustrated the advantages of our solution with some examples. Our proposed controller is implemented as a reconfiguration distributed hardware component. Its modular and flexible structure make it reusable and permits implementation code generation using MDE approach. Thus, these properties help to decrease design effort and time to market of complex partial and dynamic reconfigurable systems. As future works, we plan to carry out MDE model transformations to enable automatic code generation of our distributed and adaptive controller in a targeted FPGA. The code can then be used as input for commercial tools for final FPGA

Fig. 1 .Fig. 2 .

 12 Fig. 1. Detailed block diagram of our proposed modular controller

Fig. 3 .

 3 Fig. 3. Simplified block diagram of a P frame H.264 video coding

Fig. 4 .

 4 Fig. 4. Modeling of the H.264 video coding

Fig. 5 .

 5 Fig. 5. Block diagram of the architecture of our reconfigurable system

synthesis.


Nicolas Brodu Ltsi

Fabien Lotte
email: fabien.lotte@gmail.com

Anatole Lécuyer
email: anatole.lecuyer@inria.fr

Comparative Study of Band-Power Extraction Techniques for Motor Imagery Classification

We review different techniques for extracting the power information contained in frequency bands in the context of electroencephalography (EEG) based Brain-Computer Interfaces (BCI). In this domain it is common to apply only one algorithm for extracting the power information. However previous work and our current study confirm that one may indeed expect varying degrees of success by choosing inadequate algorithms for the power extraction. Our results suggest that on average one algorithm seems superior for extracting the power information for Motor Imagery tasks : the application of a Morlet wavelet on the raw EEG signals, with the time-frequency resolution tradeoff selected by cross-validation.

I. INTRODUCTION

This paper proposes a comparative study of various techniques for estimating the power information contained in different frequency bands of a real signal, in the specific context of motor imagery (i.e., limb movement imagination) for Brain-Computer Interfaces (BCI) [START_REF] Pfurtscheller | Motor imagery and direct braincomputer communication[END_REF]. The specificity of the BCI field is to deal with very noisy data captured by Electro-Encephalo-Graphy (EEG). The goal is to infer the user's mental state from his/her EEG signals. In a BCI setup we consider in this study the task is to find whether the user imagines left hand or right hand movements, in order to build an interactive Human-Computer interface based only on brain activity. Practical applications are envisionned for severly disabled people who cannot move their limbs but whose brain is still functionnal [START_REF] Pfurtscheller | Rehabilitation with brain-computer interface systems[END_REF].

One major challenge of this field is thus to extract reliable information from noisy data in real time in the form of relevant features [START_REF] Bashashati | A survey of signal processing algorithms in brain-computer interfaces based on electrical brain signals[END_REF]. These can then be passed on to classification techniques for identifying the user's mental state. Physiological arguments suggest that the µ and β frequency bands (around 8-12 Hz and 16-24 Hz) are especially relevant for discriminating motor activity. A common approach in the BCI field is thus to extract the power information from the signal over these frequency bands and use that as the classification feature.

There is a wide litterature on power estimation (see [START_REF] Flandrin | Temps-Fréquence[END_REF] for pointers). Time-frequency decomposition is in itself a wide topic of research, which has been extensively studied. Different techniques allow a practitionner to select the best trade-off to make between time and frequency resolutions, the presence of artifacts (spectral bias) or increased precision, etc. Yet most BCI framework only propose or consider a few of these techniques and with a limited range of trade-offs. A thorough comparison of different power feature extraction techniques in the context of BCI was presented in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF], together with their use in classification of motor imagery. This paper aims at complementing and reproducing independently the observations made in that survey. Notable differences between this paper and [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] include:

• Eleven out of the 13 subjets we used for our benchmark are distinct from [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF]. We used a subset of the most representative power estimation techniques based on their results and as justified in the next section. • While some of our results confirmed independently the results obtained in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF], other contrasted them. More particularly, we showed that when the time-frequency resolution is properly selected the Morlet wavelet can be the most efficient band-power extraction technique, whereas it was reported the less efficient one in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF].

• In [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] the power density is extracted along each timefrequency decomposition atom (See Fig 1), and the energy is summed over frequency components at that time. In this paper we instead prioritise the frequency representation and sum the energy over the trial duration. This is justified physiologically by considering that motor imagery related brain activity is more easily discriminable in the frequency domain than in the time domain (i.e. we are looking for frequency signatures rather than time patterns). Moreover, focusing on a frequential signature rather than on a temporal signature is a popular and efficient approach that has been used and validated by several leading BCI laboratories such as the Graz BCI or the Berlin BCI [START_REF] Pfurtscheller | The Graz Brain-Computer Interface Brain-Computer Interfaces[END_REF], [START_REF] Pfurtscheller | 15 years of bci research at graz university of technology: Current projects[END_REF], [START_REF] Blankertz | The non-invasive berlin brain-computer interface: Fast acquisition of effective performance in untrained subjects[END_REF]. • We use the log-transform of the energy in each frequency band as the discriminant feature. We decided to compute the energy on 1 Hz bands from 4 to 30 Hz, for both the C3 and C4 electrodes. One can recover the physiological µ and β bands (around 8-12 Hz and 16-24 Hz) by combining the corresponding 1Hz bands to reduce the number of features, which was done in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF]. Yet we prefer to let the classifiers exploit all available information rather than fixing the µ and β bands a priori. This is indeed justified by the well known inter-subject variability in motor-imagery based BCI [START_REF] Attina | The importance of individual features for motor-imagery based bci[END_REF]. • We used only the Linear Discriminant Analysis (LDA) classifier as it is well-used in the BCI litterature [START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF].

In [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF], a genetic algorithm had to be used to fine-tune Support Vector Machines. But we want to assess the performance of the extraction methods. Using a nonparametric classifier like LDA offers a stable basis of comparison and avoids any interference from the SVM training procedure. In the next section we detail how we extracted the features and present the corresponding techniques we used to extract the band-power information. Section III presents the data set we used and the results of our experiments. The paper ends with a conclusion and some recommendations based on our results.

II. POWER IN FREQUENCY BANDS

The next subsections describe each band-power estimation technique we used and how they were extracted in the present study. We considered using the periodogram, the spectrogram and the Morlet wavelet scalogram to cover a diversity of timefrequency samplings (see Fig 1); the Wigner-Ville distribution following [START_REF] Flandrin | Temps-Fréquence[END_REF] because it represents an alternative compromise to the spectrogram (see Section II-B); the power spectral density estimation using an auto-regressive model with the Yule-Walker algorithm for its good results in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] and because it brings another way of estimating the spectrum compared to the aformentioned approaches; and the Butterworth filter since it is a commonly used feature in the BCI litterature [START_REF] Pfurtscheller | Motor imagery and direct braincomputer communication[END_REF], [START_REF] Renard | Openvibe: An open-source software platform to design, test and use brain-computer interfaces in real and virtual environments[END_REF].

A. The spectrogram method

A usual way to extract the power in frequency bands is to use the spectrogram method [START_REF] Flandrin | Temps-Fréquence[END_REF] to get a time-frequency signal decomposition:

• The signal is convolved with a windowing function that is localised in time around t . A Fourier transform is then performed on the convolved signal, non null only on a sliding chunk centred on t, leading to a time-frequency representation X(f, t) of the original signal x(t) (see Fig

1, second case).

• The energy marginals in the frequency domain [START_REF] Auger | The timefrequency toolbox tutorial[END_REF] allow us to estimate a spectrum:

E f = t |X(f, t)| 2 •
The energy in each band B is obtained by summing the marginals in that band: E B = f ∈B E f . As explained in the introduction, a feature vector is then formed from the logarithms of the estimated energy in each 1Hz band.

B. The Wigner-Ville distribution

As shown in [START_REF] Flandrin | Temps-Fréquence[END_REF] there is a trade-off to make between increased precision for the energy estimation in each decomposition region of the time-frequency plane (see Fig 1), and the occurrence of artifacts. These artifacts denote alterations due to the frequency decomposition and not physiological artifacts such as EOG or EMG. The spectrogram introduces little artifacts, and the Wigner-Ville distribution achieves the opposite trade-off: improved accuracy over each decomposition atom, at the cost of more artifacts.

An hypothesis is that the improved localisation in frequency could compensate the presence of the artifacts for the classification tasks.

C. The Morlet wavelet scalogram

The scalogram is the third splicing of the time-frequency plane represented on Fig 1 . A mother wavelet ψ(x) function is scaled and shifted and forms a basis for the decomposition. The frequency interpretation comes from the fact that if the mother wavelet is localised at a given frequency f 0 , then a scaled version ψ(ax) is localised at frequency f0 a by the scaling property of the Fourier transform. The wavelet decomposition formula can thus be rewritten as we did in Fig. 1 to highlight f f0 as a scaling factor. In [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] the mother wavelet is carefully chosen so the relevant µ and β bands are well covered in the time-frequency plane, since they combine the energy contributions in these bands into a single feature. In the present study we use the facility provided by [START_REF] Auger | The timefrequency toolbox tutorial[END_REF] to adapt the wavelet so that f = f0 a ranges between 4 and 30 Hz: For a Morlet wavelet f 0 can be chosen arbitrarily, and the scale parameter a is adapted to maintain the selected frequency range. The influence of this free parameter however is to shape the effective resolutions on both axis in the time-frequency plane as shown in 1. In the present study we use a 10-fold cross-validation on the training set in order to select the best time-frequency resolution compromise.

D. Periodogram on the whole signal

The spectrogram method considered sliding chunks because it is primarily used for time-frequency decomposition: the energy is localised with a fixed resolution in both the frequency and time domain. However here we only perform a classification step at the end of trials. Thus we are only interested in the frequency components over a fixed signal duration, not at each instant. Applying an FFT transform to the whole signal and proceeding directly to the power estimation better matches the structure of the time-frequency decomposition we need, as depicted in Fig. 1. Technically the periodogram is thus not a time-frequency decomposition method, only a frequency decomposition on the whole signal.

Unfortunately as the signal is usually finite and non-periodic this technique introduces spurious artifacts in otherwise unrelated frequency bands (known as spectral bias). A workaround is to apply a window function on the whole signal, but we loose some information as the window would smoothly damp it to 0 on each end. The spectrogram copes with this by sliding each chunk in time, thus getting a correct estimate over each chunk centre.

A second workaround is to remove the 0 and low-frequency components so as to detrend the signal: This operation centres the signal around 0, and the periodisation implicitly operated by the FFT results in less artifacts. We applied this second workaround by pre-filtering the signals with a linear phaseresponse high-pass FIR filter: at least 50dB attenuation at 0Hz and less than 1dB change at 4Hz in order to preserve the physiologically relevant bands.

The periodogram is then simply computed:

• The signal x(t) is taken as a whole, without applying a window function, but filtered as aforementioned. • A Discrete Fourier Transform is applied, giving a frequency representation X(f ) of the whole signal, including artifacts due to the non-periodic nature of the signal (attenuated by the FIR pre-filtering). • Power is extracted in narrow bands B (of 1 Hz in this example):

E B (X) = f ∈B |X(f )| 2 .

E. Band-pass filtering with a Butterworth IIR filter

It is also possible to directly apply a band-pass filter to the signal x(t) in order to damp out frequencies outside each selected band B: the resulting filtered signal y(t) presumably contains only the selected frequencies. Assuming a perfect block filter, it is then possible to estimate the energy in these bands by squaring the magnitude of the time series samples directly:

E band = t |y(t)| 2 .
This is a frequent feature found in the BCI field [START_REF] Pfurtscheller | Motor imagery and direct braincomputer communication[END_REF]. This is the method that is performed by default in the OpenViBE software [START_REF] Renard | Openvibe: An open-source software platform to design, test and use brain-computer interfaces in real and virtual environments[END_REF], with a Butterworth IIR filter of order 4, and that was also used for [START_REF] Brodu | Multifractal feature vectors for brain-computer interfaces[END_REF].

In the present study we select the order of the Butterworth filter on the training data set in a 10-fold cross-validation step.

F. Power spectral density estimation using an auto-regressive model

The Yule-Walker algorithm can be used to train an autoregressive model. The power spectral density is then derived from the model's response. As in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] we estimate the order of the auto-regressive model with a 10-fold cross-validation on the training data set (using a least-squarred error hyperplane fit as the classifier). However the feature we extract is the log-transform of the energy in each of the bands. We did not consider using directly the auto-regressive model coefficients as the feature vector as this was covered already in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF].

III. EVALUATION

A. Data sets

We have applied the power estimation techniques on 13 subjets. Twelve subjets are taken from the BCI competitions II, III and IV1 [START_REF] Blankertz | The BCI competition 2003: Progress and perspectives in detection and discrimination of EEG single trials[END_REF] [START_REF] Blankertz | The BCI competition III: Validating alternative approaches to actual BCI problems[END_REF]. We also captured data with the software platform OpenViBE [START_REF] Renard | Openvibe: An open-source software platform to design, test and use brain-computer interfaces in real and virtual environments[END_REF] 2 .

1) BCI competition data sets: These data set were captured and provided by the Department of Medical Informatics, Institute for Biomedical Engineering, University of Technology Graz [START_REF]Estimating the mutual information of an eegbased brain-computer-interface[END_REF]. All the data we used follow the same experimental protocol, with minor differences explained below. In a first step a visual clue is presented to the subjet indicating left or right. The subject shall then imagine moving the corresponding hand for a fixed duration. The results presented below estimate the power over all this imaginated period. In some experiments a visual feedback is provided to the user during this period to give him a clue as to how well the signal can be detected with an internal classifier. A few seconds separate each trial, and the whole experiment is repeated with another random left or right direction.

We used data from three installments of the competitions:

• BCI competition II, data set III: This set contains data captured from only one subjet. 280 trials were equally spread in the left and right hand class. The C3, C4 and Cz electrodes were sampled at 128 Hz. Following [START_REF] Lemm | BCI competition 2003-data set III: probabilistic modeling of sensorimotor mu rhythms for classification of imaginary hand movements[END_REF] we used only the C3 and C4 electrodes, during the trial period 3s to 9s. • BCI competition III, data set IIIb. Three subjets are provided, but we used only the subjects labelled S4 and X11 due to errors in the data of the remaining subject 3 . Again we used only the C3 and C4 electrodes, which were sampled in this data set at 125 Hz, during the presentation period 4s to 7s. • BCI competition IV, data set IIb. Nine subjets are provided. Although the experimental protocol remains the same the competition organisers warned against the presence of occular artefacts that are not present in the other two previous competitions. One goal of this data set was precisely to assert how well techniques can cope with the presence of the artifacts. In the current study we simply ignore them: we want to compare various ways to extract the power in frequency bands, irrespectively of whether the data contain artifacts or not. Should some power extraction techniques be less sensitive than the others to this kind of noise, then we are happy to have our benchmark include this more realistic case. The total of the visual cue and imaginary durations in the provided signals is 4s or 4.5s depending on the trials. All these data sets were preprocessed by the Graz team with a band-pass filter between 0.5 and 30 Hz. We unfortunately do not have access to the raw data for comparing directly the various power extraction techniques. We thus assume that the filter they used does not interfere too much in the bands we used for comparison: between 4 and 30 Hz.

2) OpenViBE data set: We captured data from one subject using the same protocol as for the BCI competition experiments.

We recorded 560 trials over two weeks, equally spread on the left and right hand movement classes, each trial lasting 5s. We then randomly selected half the trials for training with the constraint of maintaining the same number of left and right instances. The remaining trials were used as the test set. We recorded 10 electrodes at 512 Hz : C3, C4, FC3, FC4, C5, C1, C2, C6, CP3, CP4 (nasion reference). We then applied a discrete Laplacian spatial filter [START_REF] Hjorth | An on-line transformation of eeg scalp potentials into orthogonal source derivations[END_REF], [START_REF] Mcfarland | Spatial filter selection for EEG-based communication[END_REF] over the C3 and C4 electrodes.

B. Comparison of the power extraction techniques

Table I shows the result of classifying the signals from the test sets using each of the techniques presented in Section II with a Linear Discriminant Analysis classifier, as this is a classifier of choice with good success records in the BCI litterature [START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF]. The parameters described in each respective method subsection (ex: Butterworth filter length, auto-regressive model order, etc) were determined by searching the parameter space and optimising the 10-fold cross-validation accuracy on the training set.

We compare each method using the classification accuracy of the LDA classifier. In addition we also note for each subject 3 See the notes at http://www.bbci.de/competition/iii/desc_IIIb_ps.html which is the maximal performance obtained with any of the power-band extraction method. The difference to that maximal performance, on average, allows us to assert the performance of the power-band method itself: an ideal extraction technique would always be the best one. In a real case there are fluctuations between subjects and experimental frameworks, but a method that is better than the other on average (and reliably so, i.e. with low deviation from that average) shall still be prefered. These results are reported on the last two lines of Table I.

We also compared the influence of performing a pre-filtering on the data. Some methods indeed benefit from such prefiltering. As explained in section II-D, detrending the data by removing low-frequency components increase the periodogram accuracy by removing spectral bias. Similarly, training an auto-regressive model only on the frequencies of interest leads to better results than training it on irrelevant frequency components. On the other hand some methods like the Morlet Wavelet are intrinsically localised in time and frequency, so the pre-filtering does not necessary help in this case. Morevover the filters we used incur a 250ms delay on the data, which can be penalizing for real-time applications. Table I mentions "hp" for a high-pass filter where the 0Hz component is nullified, and a "bp" filter with at least -50 dB attenuation at 0Hz and above 50Hz, in each case with less than 1dB change in the 6-30 Hz band.

Only 2 subjects out of the 13 subjects considered here are common with [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF]. Our results both contrast with and confirm some points of [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF]. One conclusion there is that Yule-Walker estimation of an auto-regressive model, which is then used to estimate the power density spectrum, is adequate for BCI tasks. In our experiments this technique performed well only when pre-filtering the data with a low or band-pass filter. It globally turns as the second-best technique (3.8% below the maximum accuracy on average) with a relative stability (3.4% standard deviation) of that performance, ex-aequo with the standard Spectrogram estimation. The Spectrogram leads here to a a better stabililty (3% standard deviation), but given the low number of subjects to perform the statistical analysis both methods seem equally suited.

However our results also contrast with [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF]. They noticed that "Morlet CWT produced clearly the least separable feature" amongst the category of features they considered. On the contrary, we here find that the Morlet wavelet produces very good results, better than the auto-regressive model technique on the test sets, and directly on the raw data (2.6% below the maximal accuracy on average vs 3.8%, and 2.4% vs 3% for the stability). This difference is probably explained by the way we use the wavelet. In [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] the approach is to pre-tune the wavelet time-frequency resolution in order to concentrate on specific frequency µ and β bands a priori, and only keep the extracted values at these frequencies. In the present study we cross-validate the wavelet time-frequency resolution to adapt it to the data, and we keep all frequencies available for the classifier. Thus we claim that with adequate pre-processing, the Morlet wavelet shall be considered as a technique of choice Bold figures indicate the maximum accuracy for each subject, gray background results within 1.5% to the maximum. "hp" stands for highpass prefiltering of the data (see the main text), "bp" indicates a band-pass filtering, and "raw" represents the unfiltered data. δmax is the difference to the maximum accuracy obtained for each subject. for extracting power information for BCI tasks.

IV. CONCLUSION

This study has reproduced and extended the comparisons performed in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] to identify the best method for extracting the power in given frequency bands, for BCI purposes. We have evaluated the features obtained with the various band power estimation methods on motor imagery data from several EEG data sets, including the effect of pre-filtering the data. Our results both confirm and contrast with [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF]. We confirmed that using the Yule-Walker algorithm for estimating an autoregressive model is a good technique for estimating the power in selected frequency bands for BCI tasks. But we have found, unlike [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF], that the Morlet wavelet seems better suited. This difference between the two comparative studies is probably due to our cross-validation of the Morlet wavelet parameters and keeping more frequency bands.

In any case the Butterworth filter, commonly used in the litterature [START_REF] Pfurtscheller | Motor imagery and direct braincomputer communication[END_REF], is here amongst the worse techniques. We thus stress as in [START_REF] Herman | Comparative analysis of spectral approaches to feature extraction for eeg-based motor imagery classification[END_REF] that better power estimation algorithms shall be considered for BCI tasks. Software for computing all the features we used is widely available on the Internet [START_REF] Auger | The timefrequency toolbox tutorial[END_REF]. The code for all the experiments performed in this study is also provided on the main author's web page 4 .

Figure 1 .

 1 Figure 1. Splicing the time-frequency plane

Table I CLASSIFICATION

 I ACCURACY (%) OF POWER EXTRACTION TECHNIQUES ON THE TEST DATA SETS

http://www.bbci.de/competition/

http://openvibe.inria.fr/?q=datasets

http://nicolas.brodu.numerimoire.net/

ACKNOWLEDGMENTS

This work has been partially supported by the French National Research Agency within the "OpenViBE2" project.


A Espargilière
email: ambroise.espargiliere@lapp.in2p3.fr

Generator level study of a Z ′ boson decaying into non-SUSY dark matter

A Randal-Sundrum inspired model predicts a dark matter candidate together with an additional Z ′ boson. The detection of this boson and the measurement of its mass and production cross section at a future linear collider is investigated at the generator level. The expected accuracy of the cross section production and of the mass of the Z ′ boson measurements are reported.

Introduction

Heavy Dirac neutrinos with Standard Model (SM) interactions are not viable dark matter candidates. Their large coupling to the Z 0 boson prevents them from providing the right Dark Matter (DM) relic density because of a much too high annihilation rate. This coupling would as well have made them easily detected in direct detection experiments. In the model described in [?], another type of interaction for Dirac neutrinos is proposed involving a new neutral gauge boson coupling to the right handed states, a Z ′ that could come from an additional U (1) ′ or SU (2)R. This Z ′ is suposed to have a small mixing with the SM Z 0 boson. In this model, called Right Handed Neutrino Model (RHNM), the right handed heavy Dirac neutrino, considered as a Weakly Interactive Massive Particle (WIMP) and denoted ν ′ in the following, becomes an interesting DM candidate because there is no coupling with the Z 0 (except through the small Z 0 /Z ′ mixing). Therefore, with a suitable parametrisation, the model leads to the correct DM relic density and complies with the experimental direct detection constraints. The Z ′ may be observed at the Large Hadron Collider (LHC) [?] or at future linear colliders as studied hereafter. The Z ′ boson is supposed to couple mainly to the heaviest states of the SM, thus to the third SM fermion family. As top factories, the LHC and future linear colliders should be the best laboratories to discover and study those processes.

In the following, the work focuses on the Z ′ production in e + e -collisions. In the first section, a cross section study leads to the choice of the channel that will be studied at the generator level. The detailed analysis is done at a centre of mass energy of 3 TeV.

Cross section study

Three channels have been considered for the study of the Z ′ production. They are discussed below and their feynman diagrams are displayed in figure 1. Cross sections have been computed using the software calcHEP (version 2.5.5) described in [?] and references therein. The model is named: "RH neutrino/RS-inspired". Centre of mass energy ranging from 0.5 -3 TeV and Z ′ mass (MZP) from 200 GeV/c 2 to 3 TeV/c 2 have been considered. The out-going Z ′ boson decays in heavy leptons or in WIMPs. At this step, all model parameters have been taken to arbitrary default values (couplings at 1, Z ′ /Z 0 mixing at 1%).

the Z ′ production cross section, for the complete MZP and centre of mass energy ranges. Some values are tabulated in table 1. The dependency of the production cross section on the Z 0 /Z ′ mixing is linear, the values are comupted with a 1% mixing. The Z 0 /Z ′ mixing and the Z ′ mass are the only non SM parameters involved in this process. The fact that the cross section keeps rising with the centre of mass energy tends to indicate a UV divergence. This hypothesis has been verified by computing the cross section at energies up to 30 TeV. The result is plotted in figure 3 and clearly confirms the UV divergence. Therefore, this effective model is not appropriate to describe this process. 2.2 Channel e + e -→ HZ ′

The Higgs boson mass has been set to 120 GeV/c 2 . The cross section of the Z ′ production in this channels is mapped in figure 4 and tabulated in table 2. This process depends on the Higgs mass and Z ′ coupling to the Higgs boson as displayed in table 3. The dependency on the Higgs mass is not critical but constitutes an extra parameter in the study. Provided the Higgs coupling to the Z ′ is not too low, this channel would have a large cross sec-tion and therefore be of great interest when results about Higgs are made available by the LHC.

Choice of the channel

The three channels considered to study the Z ′ production in high energy e + e -collisions have been discussed in the previous sections. Its production through W + W -fusion has been discarded because of its UV divergence.

Choice of the decay mode

The Z ′ has two decay modes:

1. Z ′ → invisible, where "invisible" stands for the DM candidate, since the model parameter settings have suppressed the other possibilities (if one neglects the effect of the small Z 0 /Z ′ mixing).

2. Z ′ → t t, leading to four-top final states.

The corresponding cross sections and branching ratios times cross section are shown in table 6. The SM cross section for four-top events in e + e - collision at 3 TeV is ≈ 20 ab, thus, provided the four top events can be correctly tagged and reconstructed, this channel will have almost no background. In addition, as the DM candidate has been taken with a mass below MZP/2, the decay of the Z ′ into ν ′ ν′ will be always allowed. Therefore for MZP below 2 × mtop, invisible decay channel will be the only available. In the following the study focuses on the invisible decay of the Z ′ .

Event selection

The main characteristics of signal and backgrounds are summarised below. A method to discriminate the signal events is proposed.

Signal

The signal consists of multi-jet events with a large number of particles and a certain amount of missing energy. It is expected to be very close to SM s-channel t t events which constitute the main background. The cross sections for different mass of the Z ′ are shown in table 6.

Background

The signal is extracted from a sample of events tagged as t t events. In a first step, the tagging of the top events is assumed to have a full efficiency and purity. Therefore, the backgrounds considered for this study are only the SM inclusive t t events and the following background channels are leftover:

e + e -→ W + W -(576 fb) e + e -→ Z 0 Z 0 (32 fb) e + e -→ W + W -Z 0 (33 fb) e + e -→ Z 0 Z 0 Z 0 (< 1 fb)
The present results shall then be scaled according to the real top tagging performance.

Both s and t channels t t production are considered, their respective Feynman diagrams are displayed below.

Event generation

The program calcHEP has been used to simulate the e + e -collisions and generate the interaction out-going particles. The information has been stored in the "Les Houches Event" (LHE) format and fed to a Pythia8 instance [?, ?] which propagates and decays the particles.

The statistics generated correspond to the integrated luminosity of 1 ab -1 , i.e.: about three years of operation at a luminosity of ≈ 10 34 cm -2 s -1 . In this study, the initial state energy spread takes into account the Initial State Radiation (ISR) only.

Discriminating variables

The final state variables listed in annex E are evaluated for each channel and their values are stored in a ROOT tree [?]. To give an intuitive idea of which variables are the most discriminating, the superimposition of the final state variables for the signal and the two backgrounds together with their correlations are displayed in annex D in the form of a correlation matrix with the variable spectra on the diagonal cells. Such a view lead to the conclusion that no direct cut on the variables would lead to a satisfactory separation between the signal and the two backgrounds. Some attempts not reported here have been made to apply slanted cuts in the correlation planes of the variables. This increase in complexity called for the use of Multi-Variate Analysis (MVA) tools. The Toolkit for Multi-Variate Analysis (TMVA) framework [?, ?] offers many functionalities which will be used all along this study.

Background as well as signal events are basically t t events. It is rather obvious that variables like jet masses, or the number of visible particles (or any particle counter from annex C) will not discriminate between signal and background. The differ-ence rather lies in the event dynamics, shape and variables related to missing energy. The final variable set retained for performing the event selection is shown in the table of annex E and superimpositions of their spectrum for the signal and he two background channels, are displayed in annex H.

Classifier training and testing

A Boosted Decision Tree (BDT) has been trained separately on the two backgrounds. The input variables spectra are displayed in annex H.1 for background 1, and in H.2 for background 2. The training and testing have been done on 2×4500 events of each type. The BDT settings are shown in table 8 and the performance of the training are given in table 9. The Receiver Operating Characteristic (ROC) integral is a good indicator of the classifier's performance but, in the present case, the most important aspect is the discriminating power of the classifier when cutting at the maximum significance. Table 9 shows that in our case the classifier gives good performance compared to a simple cut on the most discriminating variable. The significance is given by:

S = ǫ N signal ǫ N signal + B N bg.
,

where ǫ is the signal efficiency at the given cut, N signal is the total number of event signal, B is the background efficiency and N bg. the total number of background events. The BDT output for both discrimination steps are shown in figures 6(a) and 6(b) for signal vs. background 1 and signal vs. background 2 respectively. The outputs computed with the training event sample are superimposed to those computed with the testing one. A small difference between the two spectra is noticeable and is reflected by the low Kolmogorov-Smirnov probability. This indicates a slight over-training of the algorithm. A higher prune strength can avoid this (a value of 95 instead of 90 is sufficient here) but at the price of lower separation even on the test sample. Therefore, it has been considered reasonable to bear this slight overtraining. The global efficiency of the event selection is 0.884 ± 0.004. In the next section, the trained MVA methods are used to apply the two successive cuts given in table 9 on the different event samples to select the event sample from which the Z ′ production cross section and mass is to be measured.

Measurements

Method

The mass of the Z ′ is measured using the event invariant mass spectrum obtained after the event selection described above. Given that the Z ′ decays to invisible and tend to be emitted at rest or very low energy, the event total invariant mass spectrum should show a sharp upper edge corresponding to the value of the nominal centre of mass energy of 3 TeV/c 2 minus the mass of the Z ′ . Nevertheless, the ISR induces a smearing of the event total invariant mass spectrum and then a smoothing of the spectrum upper edge. A function, called "smooth gate" in the following, defined in eq. 1, is fitted to the event total invariant mass spectrum without presence of the signal. This way the background contribution is modelized and is subtracted from the total spectrum. The area of the remaining spectrum is then proportional to the cross section of the e + e -→ t tZ ′ process. To measure the Z ′ mass, the "smooth gate" is fitted to this spectrum and the s2 parameter, corresponding to the upper slope inflexion point, is related to MZP through the relation 4 established in section 4.2.2. The "smooth gate" is defined by:

f (x) = A (e s 1 -x p 1 + 1)(e x-s 2 p 2 + 1) , (1)
where si, i ∈ {1, 2} are the inflection point abscissae and pi, i ∈ {1, 2} are the slope steepnesses, see graphical representation with arbitrary values of the parameters in figure 7. After event selection, the cross section is calculated subtracting the background and correcting for event selection efficiency and then dividing by the 1 ab -1 luminosity.

Results

The analysis process is detailed in the next section for MZP = 200 GeV/c 2 . Then, in section 4.2.2, this analysis is repeated for several values of MZP and its measurement is described. Finally, in section 4.2.3, the resolution on the Z ′ mass is estimated.

Details for MZ

P = 200 GeV/c 2
All events from the two backgrounds (20000 and 5000 for background 1 and 2 respectively) and the signal (15000) are tested with the two BDT cuts. The cut for background 2 rejection is applied first. The event invariant mass spectrum is plotted with the remaining events and displayed in figure 8 The background contribution alone is plotted in figure 9 and the "smooth gate" function is fitted to it. The fit results are displayed in figure 9 and summarized in table 10. The total histogram (corresponding to the addition of the three contributions of the stack histogram from figure 8) is corrected by removing from each bin the value returned by the fitted "smooth gate" function evaluated at the bin centre abscissa. The corrected histogram is displayed in figure 10 It may be not obvious at first glance how the Z ′ mass can be read from those fit parameters, but remember that si , i ∈ {1, 2} are the inflexion point abscissae. With no smearing at all the upper limit of the spectrum would show a steep end-point at 3 TeV-MZP. This corresponds to the limit when p2 goes to 0. Then p2 can be considered as an indicator of the smearing strength. In the same view, s2 would match ECM -MZP in the low smearing limit. Therefore, there should be a set of parameters relating s2 to MZP as:

s2 = k -αMZP , (2)
where k and α are parameters to determine.

In the next section, the same procedure has been repeated for various values of MZP to determine the values of the parameters of eq. 2.

Sensitivity to MZ P

The procedure described in the previous section is applied for MZP from 200 to 700 GeV/c 2 with a constant statistics. The use of a constant statistics is meant to probe the s2 sensitivity to MZP only. The resulting event invariant mass spectra are displayed in figure 11

Resolution on MZ P

Thanks to eq. 4, it is now possible to predict the Z ′ mass from a t t event invariant mass spectrum to which a "smooth gate" function is fitted (see definition in eq. 1). As MZP = (k -s2)/α, the error margin on this prediction reads:

∆MZP = ∆k α constant term ⊕ MZP∆α α linear term ⊕ ∆s2 α fit term . (5) = 7.8 GeV/c 2 ⊕ 0.95∆s2 ⊕ 1.9%MZP
See annex G for the detailed calculation of this error.

The errors ∆k and ∆α are inherent to the measurement method itself. Hence, the constant and the linear term of eq. 5 describe systematic uncertainties. They may be reduced by refining the analysis on the s2 dependency on MZP described in the previous section (4.2.2). For instance, using more statistics for each value of MZP and/or using more values of MZP are possible refinements. Whereas the error ∆s2 comes from the event invariant mass spectra fitted with the "smooth gate" function and can only be improved by the accumulation of more statistics. Thus, the so-called "fit term" reflects the statistical uncertainty on MZP.

To evaluate correctly the expectable resolution on MZP, the statistics predicted by the model must be used. Table 13 sums up the main features of the analysis using the model predicted statistics: for each mass, the table shows the number of events passing the MVA cuts and surviving the background subtraction, the value of the s2 parameter with its error and the significance of the new physics signal. It appears here that above MZP =500 GeV/c 2 the event yield is not sufficient to perform a reliable fit and the s2 returned values do not make sense any more, thus the Z ′ mass would not be measurable. Nevertheless, given that the SM predicts about 4700 events passing the cuts, the Z ′ discovery may still be allowed in virtue of the obtained significance, the measured cross sections are shown in table 14. The present method tend to underestimate the process cross section for low MZP and slightly over-estimate it at higher values. For MZP between 200 and 500 GeV/c 2 a resolution is predicted using eq. 5 and values from table 13. The result is plotted in figure 13.

Conclusion

A generator level analysis of the predictions of a Randall-Sundrum-inspired model furnishing right handed neutrino as a dark matter candidate coupling to an additional gauge boson Z ′ (either SU (2)R or U (1) ′ ) has been carried out. The study focused on the e + e -→ t tZ ′ channel, with the Z ′ decaying into WIMPs. A method to measure the Z ′ mass has been proposed. If the Z ′ has a mass below 500 GeV/c 2 , it may be measured with a precision better than 10% with 1 ab -1 of data, and therefore may guide a threshold energy scan to reach better resolution. Shall the Z ′ mass be around 600 GeV/c 2 or higher, a discovery would still be possible but no mass estimation could be cast. The results given here are subject to review when a full simulation or at least a realistic smearing of the MVA input variables is available.

like especially to thank Jean-Jacques Blaising for advices, useful discussions and precious help.

A Dark matter relic density and detection rate

The DM relic density and the direct detection event rate are calculated using micrOMEGAS2.2, with MZP = 200 GeV/c 2 and m ν ′ = 85 GeV/c 2 .

Ω ν ′ = 0.11 WIMP-nucleon spin independent scattering cross section:

σ ν ′ -p = 5.8 • 10 -9 σ ν ′ -n = 5.3 • 10 -7
Direct detection event rate for 73 Ge and 131 Xe targets:

73 Ge: 0.058 /day/kg 131 Xe: 0.12 /day/kg B Summary of the cross sections and statistics for the signal.

MZP ( GeV/c G MZ P measurement error calculation.

As described in section 4.2.2, the Z ′ mass relates to the upper inflection point abscissa of the "smooth gate" fitted on the event invariant mass spectrum, denoted s2. The error is derived from this relation:

MZP(s2) = k -s2 α ⇒ ∆ 2 MZP (MZP) 2 = ∆ 2 (k -s2) (k -s2) 2 + ∆ 2 α α 2 with ∆ 2 (k -s2) = ∆ 2 k + ∆ 2 s2 ⇒ ∆ 2 MZP (MZP) 2 = ∆ 2 k (k -s2) 2 + ∆ 2 s2 (k -s2) 2 + ∆ 2 α α 2 but (k -s2) 2 = α 2 (MZP) 2 ⇒ ∆ 2 MZP = ∆ 2 k α 2 + ∆ 2 s2 α 2 + MZP∆ 2 α α 2 ⇒ ∆MZP = ∆ 2 k α 2 + ∆ 2 s2 α 2 + MZP∆ 2 α α 2
which also reads :

∆MZP = ∆k α ⊕ ∆s2 α ⊕ MZP∆α α
H MVA input variable spectrum

Figure 1 : Z ′ production channels 2 . 1

 121 Figure 1: Z ′ production channels

Figure 2 :

 2 Figure 2: Map of the Z ′ cross section production in W + W -fusion process.

Figure 3 :

 3 Figure 3: Divergence of the Z ′ cross section production in W + W -fusion process in the RHNM.

Figure 4 :

 4 Figure 4: Map of the Z ′ production cross section via "pseudo-higgstrahlung"

Figure 5 Figure 5 :

 55 Figure5shows the maps of the Z ′ production via emission by an out-going top quark (fig.5(a)) or an out-going τ lepton (fig.5(b)). Cross section values are also tabulated in table4. The behaviour of the Z ′ production cross section is very similar in both cases. The difference resides in the mass of the involved fermion. The cross section maps clearly indicate that the detection of a rather light Z ′ (below 1 TeV) is favoured.

Figure 6 :

 6 Figure 6: Response of BDT trained for bg1 rejection (6(a)) and for bg2 rejection (6(b)) .

Figure 7 :

 7 Figure 7: Smooth gate function with arbitrary parameter values.

  . The contributions from the different sources are represented with different colours. The second background has almost vanished and the signal represents a clear excess with respect to the SM background.

Figure 8 :

 8 Figure 8: Event total invariant mass spectrum for signal and both backgrounds (MZ P =200 GeV/c 2 ).

Figure 9 :

 9 Figure9: Event total invariant mass spectrum after event selection, for background only, fitted with the "smooth gate" function.

Figure 11 :

 11 Figure 11: Event invariant mass spectra for various Z ′ masses.

Figure 12 :

 12 Figure 12: s 2 parameter versus MZ P .

Figure 13 :

 13 Figure 13: MZ P error vs. MZ P .

H. 1 Signal versus background 1 Missing

 11

Table 1 :

 1 Cross section of the Z ′ production through W + W -fusion.

Table 2 :

 2 Cross section of the Z ′ production through "pseudo-higgstrahlung".

	MZP √ s = 3 TeV √ s = 1 TeV √ s = 0.5 TeV 157.5 fb 200 2.2 fb 25.6 fb	500 0.4 fb 0.2 fb 1000 8.1 fb ---

Table 3 :

 3 Cross section dependency on Higgs mass and Higgs coupling to the Z ′ .

	H 0 mass	120	150	300
	gZH = 0.5	6.4 fb	6.3 fb	5.2 fb
	gZH = 1	25.6 fb 25.1 fb 20.8 fb
	gZH = 2	102 fb	100 fb	83 fb

2.3 Channel e

+ e -→ l + l -Z ′

Table 4 :

 4 Cross section of the Z ′ production via emission by an out-going fermion.

	case of the top quark
	MZP √ s = 3 TeV √ s = 1 TeV √ s = 0.5 TeV	200 3.4 fb 3.0 fb 0.04 fb 500 1.0 fb --	1000 0.2 fb --
	case of the τ lepton
	MZP √ s = 3 TeV √ s = 1 TeV √ s = 0.5 TeV 2.3 fb 200 2.5 fb 4.4 fb	500 0.7 fb 0.2 fb -	1000 0.1 fb --

Table 5 :

 5 Summary of the model parameters set to non default values

	parameter name	description	value
	mixzzp	Z 0 /Z ′ mixing	0.01
	gZp	Z ′ -ν ′ coupling	3
	gtr	Z ′ -tR coupling	3
	gUp	Z ′ -τ coupling	0
	gtl	Z-(tL, bL) coupling	0
	gll	Z ′ -τ coupling	0

Table 6 :

 6 Total process cross sections and Z ′ decay branching ratios times cross section.

	MZP	Cross sect. σ Z ′ →inv. σ Z ′ →t	t
	200 GeV/c 2	14.9 fb	14.9 fb	0
	300 GeV/c 2	8.7 fb	8.7 fb	0
	400 GeV/c 2	5.6 fb	2.4 fb	3.2 fb
	500 GeV/c 2	3.7 fb	1.2 fb	2.5 fb
	600 GeV/c 2	2.6 fb	0.8 fb	1.8 fb
	700 GeV/c 2	1.8 fb	0.5 fb	1.3 fb

  1. background 1: e + e -→ t t

	e -		t
	γ, Z
	e +		t
	2. background 2: e + e -→ t t + ν	ν
	e -	νe
	W -	γ, Z	t
	W +			t
	e +	νe
	Their respective cross sections are summed up
	in table 7.

Table 7 :

 7 background channel cross sections.

	Channel	Cross sect.

e + e -→ t t 19.87 fb e + e -→ t t + ν ν 5 fb

Table 8 :

 8 BDT parameters

	Boost type Adaptative
	Number of trees 1000
	Maximum depth 20
	Prune method Cost complexity
	Prune strength 90

Table 9 :

 9 BDT performance. ROC stands for the ROC curve integral, S max denotes the maximal significance, ǫ is the signal efficiency and B is the background efficiency.

		BDT-BG1 cut BDT-BG2 cut
		Sig. vs Bg1	Sig. vs Bg2
	ROC	0.955	0.971
	Smax	99.6	113.8
	cut @Smax	-0.0378	-0.3391
	ǫ@Smax	90.8%	97.4%
	B@Smax	18.0%	14.6%

  and the "smooth gate" is fitted to it. The corresponding fit results are displayed on the picture and summarized in table 11.

	) 2 Number of events / (30 GeV/c	Entries Mean RMS / ndf 2 χ A 1 s 1 p 2 s 2 p 160 Entries 40 60 80 100 120 Mean 140 RMS / ndf 2 χ A 1 s 1 p 2 s 2 p		4714 1967 407.5 118.5 / 73 3.6 ± 129.8 19.8 ± 1491 8.1 ± 197.8 8.4 ± 2556 3.87 ± 62.55 4714 1967 407.5 118.5 / 73 3.6 ± 129.8 19.8 ± 1491 8.1 ± 197.8 8.4 ± 2556 3.87 ± 62.55
		20
		0 0	500	1000	] 2 3000 Event invariant mass [GeV/c 1500 2000 2500

Table 10 :

 10 Fit results for background contribution

			χ 2 /ndf 118.5/73
				A s1 p1 s2 p2	129.8 ± 3.6 1491 ± 19.8 197.8 ± 8.1 2556 ± 8.4 62.6 ± 3.9
	) 2 Number of events / (30 GeV/c	100 150 200 250 300	Entries Mean RMS / ndf 2 χ A 1 s 1 p 2 s 2 p Entries Mean RMS / ndf 2 χ A 1 s 1 p 2 s 2 p	12094 1824 473.9 433.8 / 95 4.5 ± 286.4 14.2 ± 1263 6.2 ± 245.9 5.1 ± 2529 2.59 ± 54.65 12094 1824 473.9 433.8 / 95 4.5 ± 286.4 14.2 ± 1263 6.2 ± 245.9 5.1 ± 2529 2.59 ± 54.65
		50
		0 0	500	1000	] 2 3000 Event invariant mass [GeV/c 1500 2000 2500
	Figure 10: Event total invariant mass spectrum
	after event selection and background subtrac-
	tion fitted with the "smooth gate" function for
	MZ P = 200 GeV/c 2 .

Table 11 :

 11 Fit results for the corrected event invariant mass spectrum.

	χ 2 /ndf 118.5/73
	A s1 p1 s2 p2	286.4 ± 4.5 1263 ± 15 245.9 ± 6.2 2529 ± 5.1 54.7 ± 2.6

Table 12 :

 12 Summary of the "smooth gate" upper inflexion point abscissae (s 2 )

	MZP	s2
	200 300 400 500 600 700	2530 ± 5.4 2458 ± 5.6 2340 ± 6.6 2235 ± 6.7 2119 ± 7.5 2005 ± 9.4

Table 13 :

 13 Z ′ events using model predicted statistics

	MZP events yield	s2	significance
	200 300 400 500 600 700	12094 7411 2150 1192 815 642	2529 ± 5.1 2448 ± 7.7 2342 ± 19.7 2282 ± 39.5 913 ± 82 962 ± 78	92.7 66.6 24.6 13.5 8.4 5.6

Table 14 :

 14 Z ′ cross section measurement result

	MZP True cross sect. Measured cross sect.
	200	14.9 fb	13.7 fb
	300	8.7 fb	8.4 fb
	400	2.4 fb	2.4 fb
	500	1.2 fb	1.3 fb
	600	0.8 fb	0.9 fb
	700	0.5 fb	0.7 fb

  Channel e + e -→ ν νZ ′ . . . . . . . . 2 2.2 Channel e + e -→ HZ ′ . . . . . . . . 2 2.3 Channel e + e -→ l + l -Z ′ . . . . . . . Method . . . . . . . . . . . . . . . . 5 4.2 Results . . . . . . . . . . . . . . . . . 5 4.2.1 Details for MZP = 200 GeV/c 2 5 4.2.2 Sensitivity to MZP . . . . . . 6 4.2.3 Resolution on MZP . . . . . . 7 Signal versus background 1 . . . . . 12 H.2 Signal versus background 2 . . . . . 13

		F BDT variable importance ranking	11
		G MZP measurement error calculation. 11
	Contents	H MVA input variable spectrum H.1	12
	1 Introduction	1
	2 Cross section study	1
	2.1 2
	2.4 Choice of the channel . . . . . . . .	3
	2.5 Choice of the decay mode . . . . . .	3
	3 Event selection	3
	3.1 Signal . . . . . . . . . . . . . . . . .	3
	3.2 Background . . . . . . . . . . . . . .	3
	3.3 Event generation . . . . . . . . . . .	4
	3.4 Discriminating variables . . . . . . .	4
	3.5 Classifier training and testing . . . .	4
	4 Measurements	5
	4.1 5 Conclusion	7
	Acknowledgment	8
	Bibliography	8
	A Dark matter relic density and detec-
	tion rate	9
	B Summary of the cross sections and
	statistics for the signal.	9
	C Index of the final state variables	9
	D Final state correlation matrix	10
	E Summary of the analysis variables	11

  2 ) Cross sect. Width(Z ′ → inv.) Width(Z ′ → t t) Br(Z ′ → inv.) Z ′ → inv.

	D Final state correlation matrix E Summary of the analysis variables
	Variable name		Description	Discriminating power
	missPT		Missing transverse momentum		3.462e-01
	missE		Missing energy		3.114e-01
	jet1E		First jet energy		3.012e-01
	jet2E		Second jet energy		2.496e-01
	jet2P		Second jet momentum		2.365e-01
	missP		Missing momentum		2.250e-01
	jet1P		First jet momentum		2.138e-01
	jetPTdiff:=jet1P-jet2P	jet PT diffrence		2.122e-01
	missA		Angle of missing momentum		1.134e-01
	jet2PT		Second jet transverse momentum		1.078e-01
	sphericity		Sphericity		1.024e-01
	aplanarity		Aplanarity		9.443e-02
	jetEratio:=jet2E/jet1E	Jets energy ratio		9.302e-02
	Jet1PT		First jet transverse momentum		7.763e-02
	jetEdiff:=jet1E-jet2E	Jets energy difference		6.089e-02
	200 300 F BDT variable importance ranking 14.9 fb 19.4 GeV 0 GeV 8.7 fb 32.8 GeV 0 GeV		100% 100%	15000 8700
	400 500 600 700	5.6 fb 3.7 fb 2.6 fb 1.8 fb		45.5 GeV 57.9 GeV Rank Variable 70.1 GeV 1 missPT 82.3 GeV 2 missE 3 jet1P	60.3 GeV 115.0 GeV Importance 162.0 GeV 1.929e-01 205.5 GeV 1.366e-01 1.069e-01		43.0% 33.0% 30.2% 28.6%	2408 1240 785 515
	4 C Index of the final state variables jet1E 8.365e-02 5 jet1PT 6.978e-02
	Variable name ttM evtM missE missP missA missPT missM sphericity aplanarity	Description missA t t system invariant mass 5.944e-02 missP 5.346e-02 Event invariant mass jet2E 4.576e-02 Missing energy jetEratio 4.288e-02 Missing momentum jetPTdiff 4.001e-02 Cosine of polar angle of missing momentum 6 7 8 9 10 11 jet2PT 3.919e-02 Missing transverse momentum 12 jet2P 3.636e-02 Missing mass 13 sphericity 3.534e-02 Event sphericity 14 jetEdiff 2.973e-02 Event aplanarity 15 aplanarity 2.799e-02	Range 0-3000 GeV/c 2 0-3000 GeV/c 2 0-3000 GeV 0-3000 GeV/c -1-1 0-3000 GeV/c 0-3000 GeV/c 2 0-1 0-0.5
	jet1P			First jet momentum		0-3000 GeV/c
	jet1PT			First jet transverse momentum		0-3000 GeV/c
	jet1E			First jet energy		0-3000 GeV
	jet1M			First jet mass		0-3000 GeV/c 2
	jet2P			Second jet momentum		0-3000 GeV/c
	jet2PT			Second jet transverse momentum		0-3000 GeV/c
	jet2E			Second jet energy		0-3000 GeV
	jet2M			Second jet mass		0-3000 GeV/c 2
	nPart nTrack nHadr nHadrC nGamma nElect nMu			Number of particles Number of charged particles Number of hadrons Number of charged hadrons Number of photons Number of electrons Number of muons		∈ N ∈ N ∈ N ∈ N ∈ N ∈ N ∈ N

Acknowledgements

The author would like to thank Geneviève Bélanger, Marco Battaglia and Géraldine Servant for their help in guiding this study at the beginning (the choice of the channel, the model parameter values to use...). Many thanks as well for Sasha Puckov and Geneviève Bélanger again, for answering all my questions about the RHNM and calcHEP and to Peter Skands for helping me fix the interface between calcHEP and Pythia8. The author would