propagate cluster error parametrization
[u/mrichter/AliRoot.git] / TRD / AliTRDtestBeam.cxx
CommitLineData
0eedb47d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
00d52803 18////////////////////////////////////////////////////////////////////////////
19// //
20// Class to handle the test beam data of 2007 //
21// //
22// Authors: //
23// Sylwester Radomski (radomski@physi.uni-heidelberg.de) //
24// Anton Andronic (A.Andronic@gsi.de) //
25// //
26////////////////////////////////////////////////////////////////////////////
27
28//#include <iostream>
29//#include <fstream>
30//#include <sys/types.h>
31//#include <sys/stat.h>
32//#include <fcntl.h>
33//include <unistd.h>
0eedb47d 34
987ba9a3 35#include "AliTRDrawStream.h"
0eedb47d 36#include "AliRawReaderMemory.h"
00d52803 37#include "AliTRDtestBeam.h"
0eedb47d 38
39ClassImp(AliTRDtestBeam)
40
6511d650 41const Long_t AliTRDtestBeam::fgkFileHeadSize = 544; // ?
42const Long_t AliTRDtestBeam::fgkEventHeadSize = 68; //?
43const Long_t AliTRDtestBeam::fgkLdcHeadSize = 68; //?
44const Long_t AliTRDtestBeam::fgkEquipHeadSize = 28; //
45const Int_t AliTRDtestBeam::fgkVmeIn =1; //VME event in
46const Int_t AliTRDtestBeam::fgkSimIn =1; //Si-strips in
0eedb47d 47
48//typedef char byte;
49
50//offsets in bytes
6511d650 51const Int_t AliTRDtestBeam::fgkPosRun = 20; //run nr. (in file and event header)
52const Int_t AliTRDtestBeam::fgkPosLength = 0; //event/equip. length
53const Int_t AliTRDtestBeam::fgkEqId = 8; //equipment id.
54const Int_t AliTRDtestBeam::fgkPosSiOff = 12; //Si data size offset (3 extra words!!!)
0eedb47d 55
56using namespace std;
57
58//____________________________________________________________________________
59AliTRDtestBeam::AliTRDtestBeam() :
60 fDataStream(0),
61 fHeaderIsRead(0),
62 fEventCount(0),
63 fLimit(4),
64 fCurrent(0),
65 fDdlOff(0),
66 fSiOff(0),
67 fQdcOff(0),
68 fDdlSize(0),
69 fFileHeader(0),
70 fEventHeader(0),
71 fEventData(0),
72 fNSi1(0),
73 fNSi2(0),
74 fCher(0),
75 fPb(0)
76{
6511d650 77 //
78 // Standard construction
79 //
0eedb47d 80
81}
01abcaa3 82//____________________________________________________________________________
056c75e1 83AliTRDtestBeam::AliTRDtestBeam(const char *filename) :
01abcaa3 84 fDataStream(0),
85 fHeaderIsRead(0),
86 fEventCount(0),
87 fLimit(4),
88 fCurrent(0),
89 fDdlOff(0),
90 fSiOff(0),
91 fQdcOff(0),
92 fDdlSize(0),
93 fFileHeader(0),
94 fEventHeader(0),
95 fEventData(0),
96 fNSi1(0),
97 fNSi2(0),
98 fCher(0),
99 fPb(0)
100{
101 //
056c75e1 102 // AliTRDtestBeam constructor
01abcaa3 103 //
104
056c75e1 105 fDataStream = new ifstream(filename, ifstream::in | ifstream::binary );
106 cout << fDataStream->is_open() << endl;
107 //fHeaderIsRead = kTRUE;
108 fHeaderIsRead = kTRUE;
109
110 fFileHeader = new Char_t[fgkFileHeadSize];
111 fEventHeader = new Char_t[fgkEventHeadSize];
112 fEventData = new Char_t[fLimit];
113
01abcaa3 114}
115
056c75e1 116//____________________________________________________________________________
117AliTRDtestBeam::AliTRDtestBeam(const AliTRDtestBeam &tb)
118 :TObject(tb),
0eedb47d 119 fDataStream(0),
120 fHeaderIsRead(0),
121 fEventCount(0),
056c75e1 122 fLimit(4),
0eedb47d 123 fCurrent(0),
124 fDdlOff(0),
125 fSiOff(0),
126 fQdcOff(0),
127 fDdlSize(0),
128 fFileHeader(0),
129 fEventHeader(0),
130 fEventData(0),
131 fNSi1(0),
132 fNSi2(0),
133 fCher(0),
134 fPb(0)
135{
6511d650 136 //
056c75e1 137 // Copy constructor
6511d650 138 //
0eedb47d 139
0eedb47d 140}
6511d650 141
0eedb47d 142//____________________________________________________________________________
6511d650 143AliTRDtestBeam::~AliTRDtestBeam()
144{
145 //
146 // Destructor
147 //
0eedb47d 148
6511d650 149 if (fDataStream) delete fDataStream;
150 if (fEventHeader) delete fEventHeader;
151 if (fFileHeader) delete fFileHeader;
152 if (fEventData) delete fEventData;
153
154}
155
156//____________________________________________________________________________
157Int_t AliTRDtestBeam::NextEvent()
158{
159 //
160 // Read the next event
161 //
0eedb47d 162
6511d650 163 Long_t dataSize=0,ldcOff; //,ldc_id,ldc2_id;
164 Long_t ldcSize,eqId; //,ev_l2;
165 Long_t eventNr,evL1;
0eedb47d 166 Long_t word;
167
168 if ( !fHeaderIsRead ) {
6511d650 169 fDataStream->read(fFileHeader, fgkFileHeadSize);
0eedb47d 170 if(fDataStream->fail()) {
171 cerr << "Error reading file header! " << endl;
172 return false;
173 }
6511d650 174 cout << " Run nr. " << Int(fgkPosRun, fFileHeader) << endl;
0eedb47d 175 fHeaderIsRead=kTRUE;
176 }
177
6511d650 178 fDataStream->read(fEventHeader, fgkEventHeadSize);
0eedb47d 179 if(fDataStream->fail()) {
180 cerr << "End of file, Event " << fEventCount << endl;
181 return false;
182 }
183
6511d650 184 dataSize = Int(fgkPosLength, fEventHeader)-fgkEventHeadSize; //?
185 eventNr = Int((4+fgkPosRun), fEventHeader); //ev.nr.
186 //cout << " Event " << eventNr <<" size "<< dataSize <<endl;
0eedb47d 187
6511d650 188 if (eventNr <= fEventCount-1) { //watch-out ...event counter starts at 1?
0eedb47d 189 cout << fEventCount << " End of file?, Event " << fEventCount << endl;
190 return false;
191 }
6511d650 192 //cout << "Run " << Int(fgkPosRun, header)<< " , Event " <<eventNr <<endl;
0eedb47d 193
194 // enough space for data?
6511d650 195 if (fLimit < dataSize) {
0eedb47d 196 delete[] fEventData;
6511d650 197 fEventData = new Char_t[dataSize];
198 fLimit = dataSize;
0eedb47d 199 }
200
6511d650 201 fDataStream->read(fEventData, dataSize);
0eedb47d 202
203 if(fDataStream->fail()) {
204 cerr << "End of file, Event " << fEventCount; // << endl;
205 return false;
206 }
207
208 //cout << " ...IDs (size) : ";
209
6511d650 210 ldcOff=0; // size of data from one DDL link
0eedb47d 211
212 for ( size_t k = 0; k < 2; k++ ) { // 2 LDCs (DDL & VME)
213
6511d650 214 ldcSize = Int(ldcOff+fgkPosLength, fEventData); //
215 //ldcSize1=(ldcSize-fgkLdcHeadSize);
216 eqId = Int(ldcOff+fgkLdcHeadSize+fgkEqId, fEventData);
217 //cout << eqId <<" ("<<ldcSize<<") ";
0eedb47d 218
6511d650 219 evL1 = Int((4+ldcOff+fgkPosRun), fEventData); //ev.nr.
220 if ( evL1 != eventNr ){
221 //cerr << "eqId " <<eqId<<" event nr. mismatch? " << eventNr <<" / "<< evL1 <<" ...LDC data size (header:68) " <<ldcSize<<endl;
0eedb47d 222 }
223
6511d650 224 if (eqId == 1024) { //DDL data
225 fDdlOff = ldcOff; //+fgkLdcHeadSize+fgkEquipHeadSize + 32;
226 fDdlSize = ldcSize;
0eedb47d 227 }
228
6511d650 229 if (eqId == 550) { //Si-strip data (+QDC)
0eedb47d 230 //cout << "550" << endl;
6511d650 231 fSiOff=ldcOff+fgkLdcHeadSize+fgkEquipHeadSize+fgkPosSiOff;
0eedb47d 232 word = Int(fSiOff, fEventData);
6511d650 233 Short_t lenSi1 = (word >> 16) & 0xffff;
234 Short_t lenSi2 = word & 0xffff;
235 fQdcOff=fSiOff+4*(lenSi1+lenSi2+1)+fgkEquipHeadSize+4;
0eedb47d 236 }
6511d650 237 else if (eqId == 1182) { //QDC first...
0eedb47d 238 //cout << "1182" << endl;
6511d650 239 fQdcOff=ldcOff+fgkLdcHeadSize+fgkEquipHeadSize+fgkPosSiOff;
240 fSiOff=fQdcOff+fgkEquipHeadSize+4;
0eedb47d 241 }
242
6511d650 243 ldcOff=ldcSize;
0eedb47d 244
245 }
246 //cout << endl;
247
248 //cout << "DDL = " << fDdlOff << endl;
249 // cout << "Si = " << fSiOff << endl;
250 //cout << "QDC = " << fQdcOff << endl;
251
252 DecodeSi();
253
254 fEventCount++; //event counter
255 return true;
256}
6511d650 257
0eedb47d 258//____________________________________________________________________________
6511d650 259Int_t AliTRDtestBeam::DecodeSi()
260{
261 //
262 // Decode the silicon detector
263 //
0eedb47d 264
265 if (fSiOff < 0) return 0;
266
267 // cout << "decoding Si data" << endl;
268
269 Long_t word;
270
271 word=Int(fSiOff, fEventData);
272 fNSi1 = (word >> 16) & 0xffff;
273 fNSi2 = word & 0xffff;
274
275 Int_t cSi=fSiOff; //
276 for (int i = 0; i < fNSi1; i++) {
277 fSi1Address[i] = ( Int(cSi, fEventData) >> 12 ) & 0x7ff;
278 fSi1Charge[i] = Int(cSi, fEventData) & 0xfff;
279 cSi+=4;
280 }
281
282 for (int i = 0; i < fNSi2; i++) { //1,for Date!
283 fSi2Address[i] = ( Int(cSi, fEventData) >> 12 ) & 0x7ff;
284 fSi2Charge[i] = Int(cSi, fEventData) & 0xfff;
285 cSi+=4;
286 }
287
288 // reconstruction
289
6511d650 290 int aLenSiX = 640;
0eedb47d 291
056c75e1 292 int amaxX=0;
293 int amaxY=0;
294
295 Int_t q, a;
00d52803 296 Int_t nst1=0,nst2=0;
297 Int_t qclX=0,qclY=0, nclX=0,nclY=0, nstX=0,nstY=0;
298 const Int_t kThr = 20;
299
300 nst1=0;
301 nstX=0;
302 nstY=0;
303 nclX=0;
304 nclY=0;
305 qclX=0;
306 qclY=0;
0eedb47d 307
0eedb47d 308 for( int i = 0; i < GetNSi1(); i++ ) {
309
310 if (fSi1Address[i] == 0) continue; // noise
056c75e1 311
312 q = fSi1Charge[i];
313 a = fSi1Address[i];
314
00d52803 315 if ( q > kThr )
056c75e1 316 {
317 if ( i > 0 && i < (GetNSi1()-1) ) {
318
319 if ( (a-fSi1Address[i+1]) == -1 &&
320 (a-fSi1Address[i-1]) == 1)
321 {
00d52803 322 nst1++;
056c75e1 323 if (a < aLenSiX) {
00d52803 324 qclX = q+fSi1Charge[i+1]+fSi1Charge[i-1];
325 nclX++;
326 nstX+=3;
056c75e1 327 amaxX = a;
328 }
329 else {
00d52803 330 qclY = q+fSi1Charge[i+1]+fSi1Charge[i-1];
331 nclY++;
332 nstY+=3;
056c75e1 333 amaxY = a;
334 }
335 i+=1;
336 }
337 else if ( (a-fSi1Address[i-1]) == 1)
338 {
00d52803 339 nst1++;
056c75e1 340 if (a < aLenSiX) {
00d52803 341 qclX = q+fSi1Charge[i-1];
342 nclX++;
343 nstX+=2;
056c75e1 344 amaxX = a;
345 }
346 else {
00d52803 347 qclY = q+fSi1Charge[i-1];
348 nclY++;
349 nstY+=2;
056c75e1 350 amaxY = a;
351 }
352 }
353 else if ( (a-fSi1Address[i+1]) == -1)
354 {
00d52803 355 nst1++;
056c75e1 356 if (a < aLenSiX) {
00d52803 357 qclX = q+fSi1Charge[i+1];
358 nclX++;
359 nstX+=2;
056c75e1 360 amaxX = a;
361 }
362 else {
00d52803 363 qclY = q+fSi1Charge[i+1];
364 nclY++;
365 nstY+=2;
056c75e1 366 amaxY = a;
367 }
368 i+=1;
369 }
370 }
0eedb47d 371 }
372 }
00d52803 373 if (nst1==2 && nstX<4 && nstY<4 ) {
056c75e1 374 fX[0] = (float)(amaxX*0.05); // [mm]
375 fY[0] = (float)((amaxY-aLenSiX)*0.05);
00d52803 376 fQx[0] = (float)qclX;
377 fQy[0] = (float)qclY;
056c75e1 378 }
379 else {
380 fX[0] = -1.;
381 fY[0] = -1.;
382 fQx[0] = 0.;
383 fQy[0] = 0.;
384 }
0eedb47d 385
056c75e1 386 // ...and Si2
387
00d52803 388 nst2=0;
389 nstX=0;
390 nstY=0;
391 nclX=0;
392 nclY=0;
393 qclX=0;
394 qclY=0;
0eedb47d 395
396 for( int i = 0; i < GetNSi2(); i++ ) {
397
398 if (fSi2Address[i] == 1279) continue; // noise
399 if (fSi2Address[i] == 0) continue; // noise
400
056c75e1 401 q = fSi2Charge[i];
402 a = fSi2Address[i];
403
00d52803 404 if ( q > kThr/2 ) //...as Si2 has 1/2 gain!
056c75e1 405 {
406 if ( i > 0 && i < (GetNSi2()-1) ) {
407
408 if ( (a-fSi2Address[i+1]) == -1 &&
409 (a-fSi2Address[i-1]) == 1)
410 {
00d52803 411 nst2++;
056c75e1 412 if (a < aLenSiX) {
00d52803 413 qclX = q+fSi2Charge[i+1]+fSi2Charge[i-1];
414 nclX++;
415 nstX+=3;
056c75e1 416 amaxX = a;
417 }
418 else {
00d52803 419 qclY = q+fSi2Charge[i+1]+fSi2Charge[i-1];
420 nclY++;
421 nstY+=3;
056c75e1 422 amaxY = a;
423 }
424 i+=1;
425 }
426 else if ( (a-fSi2Address[i-1]) == 1)
427 {
00d52803 428 nst2++;
056c75e1 429 if (a < aLenSiX) {
00d52803 430 qclX = q+fSi2Charge[i-1];
431 nclX++;
432 nstX+=2;
056c75e1 433 amaxX = a;
434 }
435 else {
00d52803 436 qclY = q+fSi2Charge[i-1];
437 nclY++;
438 nstY+=2;
056c75e1 439 amaxY = a;
440 }
441 }
442 else if ( (a-fSi2Address[i+1]) == -1)
443 {
00d52803 444 nst2++;
056c75e1 445 if (a < aLenSiX) {
00d52803 446 qclX = q+fSi2Charge[i+1];
447 nclX++;
448 nstX+=2;
056c75e1 449 amaxX = a;
450 }
451 else {
00d52803 452 qclY = q+fSi2Charge[i+1];
453 nclY++;
454 nstY+=2;
056c75e1 455 amaxY = a;
456 }
457 i+=1;
458 }
459 }
0eedb47d 460 }
461 }
462
00d52803 463 if (nst2==2 && nstX<4 && nstY<4 ) {
056c75e1 464 fX[1] = (float)(amaxX*0.05); // [mm]
465 fY[1] = (float)((amaxY-aLenSiX)*0.05);
00d52803 466 fQx[1] = (float)qclX;
467 fQy[1] = (float)qclY;
056c75e1 468 }
469 else {
470 fX[1] = -1.;
471 fY[1] = -1.;
472 fQx[1] = 0.;
473 fQy[1] = 0.;
474 }
0eedb47d 475
476 if (fQdcOff < 0) return 0;
477
478 word=Int(fQdcOff, fEventData);
479 fPb = (Double_t)((word >> 16) & 0xFFF);
480 fCher = (Double_t)((word ) & 0xFFF);
481
482 //cout << fCher << " " << fPb << endl;
483 return 1;
6511d650 484
0eedb47d 485}
486//____________________________________________________________________________
987ba9a3 487AliTRDrawStream *AliTRDtestBeam::GetTRDrawStream()
6511d650 488{
489 //
490 // Get the TRD raw stream
491 //
0eedb47d 492
987ba9a3 493 // needs AliTRDrawStream
0eedb47d 494 //cout << "Chamber reader:" << (Int_t)(fEventData+fDdlOff) << " " << fDdlSize << endl;
495 //int ifout = open("dump.dat", O_WRONLY | O_TRUNC | O_CREAT);
496 //write(ifout, (void*)(fEventData+fDdlOff+16), fDdlSize);
497 //close(ifout);
498
499 AliRawReaderMemory *reader = new AliRawReaderMemory((UChar_t*)(fEventData+fDdlOff), (UInt_t)fDdlSize);
500 reader->SetEquipmentID(1024);
501 reader->ReadHeader();
987ba9a3 502 //AliTRDrawStream::RawBufferMissAligned(kTRUE);
c29b11ff 503
987ba9a3 504 AliTRDrawStream::SetNoErrorWarning();
505 AliTRDrawStream::SetExtraWordsFix();
506 AliTRDrawStream::AllowCorruptedData();
f5fd37e8 507
987ba9a3 508 AliTRDrawStream *tb = new AliTRDrawStream(reader);
b32c70da 509 //tb->Init();
0eedb47d 510 return tb;
511 /*
512 return
513
508c0804 514 AliEawReaderMemory *rmem = data->GetRawReader();
0eedb47d 515 rmem->ReadHeader();
516
987ba9a3 517 AliTRDrawStream tb(rmem);
0eedb47d 518 tb.Init();
987ba9a3 519 AliTRDrawStream::SupressWarnings(kTRUE);
0eedb47d 520
521 */
522}
0eedb47d 523
6511d650 524//____________________________________________________________________________
00d52803 525Int_t AliTRDtestBeam::Int(Int_t i, Char_t *start) const
6511d650 526{
527 //
528 // ?????
529 //
0eedb47d 530
531 bool swap = kFALSE;
532
533 if(swap) {
534 char *q=(char*)(start+i);
535 char p[] = {q[3], q[2], q[1], q[0]};
536 return *((int*) p);
537 } else return *((int*)(start+i));
6511d650 538
0eedb47d 539}
540
541//____________________________________________________________________________