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