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