TRD/AliTRDtestBeam.cxx

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