---

CalDisplay.cxx

Go to the documentation of this file.

#include "bfemDisplay/CalDisplay.h"
#include "bfemDisplay/Enum.h"

#include "reconRootData/ReconHeader.h"

CalDisplay::CalDisplay( TNode *parent, TObjArray *peds ) {

  m_peds = peds;

  m_Gammas = 0;

  // !Note: All measurements in cm.  Numbers are half what they are in reality.  
  // !This is to accomadate the fact that TBrik accept half lengths as input variables
  // !Coordinate system centere in the center of the tower frame.

  m_CalTkrdistance  = 8.61;
  m_CalHeightOffset = -21.68;
  m_CalWidth        = 15.525 ; 
  m_CalHeight       = 13.07;
  m_LogWidth        = 1.513;
  m_LogHeight       = 1.307;
  m_LogLength       = 15.13;
  m_NumCalLayers    = 8;
  m_CalHitLength    = 1.5;



  m_CalNodeList = 0;
  m_CalHitList = 0;


  // !Create representation for the CAL.

  parent->cd();

  TBRIK* calFrame = new TBRIK( "calFrame", "calFrame", "void", m_CalWidth, m_CalWidth, m_CalHeight );

  calFrame->SetLineColor( 1 );

  m_calFrameNode  = new TNode( "calFrameNode", "calFrameNode", "calFrame", 0.0, 0.0, m_CalHeightOffset );

  m_calFrameNode->cd();

  // Top-level node for hits on logs parallel to Y-axis ( measured coordinate is X ).
  m_AllHitsX = new TNode( "AllHitsX", "AllHitsX", "calFrame" );
  m_AllHitsX->SetParent( m_calFrameNode );

  // Top-level node for hits on logs parallel to X-axis ( measured coordinate is Y ).
  m_calFrameNode->cd();
  m_AllHitsY = new TNode( "AllHitsY", "AllHitsY", "calFrame" );
  m_AllHitsY->SetParent( m_calFrameNode );


}

void CalDisplay::DrawEvents( DigiEvent *event, Recon *recon, Bool_t isRealData  ) {


  TNNode *newnode;
  THits *calLog;

  m_event = event;
  m_recon = recon;


  // !Draw displays the current event

  // !Delete the old hit node and hit arrays.

  if ( m_CalNodeList ) m_CalNodeList->Delete();
  m_CalNodeList = new TObjArray();

  if ( m_CalHitList ) m_CalHitList->Delete();
  m_CalHitList = new TObjArray();

 
  Char_t* hit_name=0; // name of elements in CalHitList

  Char_t* cal_name=0; // name of elements in CalShapeList 


  Char_t* view;

  Int_t nhit = 0;

  Int_t nCAL = m_event->getCalDigi()->GetEntries();

  for (Int_t ihit=0; ihit < nCAL; ihit++) {

    CalLog *hit = (CalLog*)m_event->getCalDigi()->At(ihit);

    LogId *log = (LogId*)hit->getLogId();

    //    Int_t id = log->getId(); // Log's ID. now from 0 to 39 (before 0 - 79)

    Float_t eTotal, e0, e1;

    Int_t layer = log->getLayer();

    Int_t column = log->getColumn();

    Int_t xyDir = log->getXY(); 

    Int_t iped = layer*10+column;

    cal_name = getCalName(iped);

    hit_name = cal_name;


    //printf(" id %i -- iped %i layer %i column %i XY %i\n",id,iped,layer,column,xyDir);


    e0 = getLogEnergy( hit, iped, 0, isRealData );

    e1 = getLogEnergy( hit, iped, 1, isRealData );


    if (e1<0 || e0<0) continue;

    eTotal = (e0 + e1)/2.0;

    // Cut on Total Energy: Display if Et>1.3 MeV (~0.1 MIP)

    if (eTotal > m_cal_lev1) {



      // !Figure out position of Log      


      Float_t position = -(e1 - e0)/(e1 + e0) * m_LogLength;

      Float_t off = (m_CalTkrdistance + m_LogHeight) + m_CalHeightOffset;

      Float_t z = -off - 2*(Float_t)layer*m_LogHeight;

      //Float_t m_Log = m_LogWidth*eTotal*100.;
      
      //if (eTotal > 0.01) m_Log = m_LogWidth; //from 100 MeV on the size of the square rapresenting 
                                             //the energy in the log is fixed and equal to the log width
      //

      Float_t x, y;

      if (xyDir == 0){ // Coord. system according to doc: BalloonTestV12_20010415.ps pag.9

        view = "Y"; // xyDir = 0 measures Y

        y = -m_CalWidth + (1.0 + 2.0*(Float_t)column)*m_LogWidth;

        x = position;

        calLog = new THits(cal_name,cal_name,"CsI",m_CalHitLength,m_LogWidth,m_LogHeight);

        SetLogColor( calLog, eTotal );

        m_AllHitsY->cd();

        newnode = new TNNode(cal_name,cal_name,cal_name,x,y,z);

        newnode->setObject( calLog );
        newnode->SetParent( m_AllHitsY );
        m_CalNodeList->Add( newnode );

      }

      else {

        view = "X"; // xyDir = 1 measures X

        x = -m_CalWidth + (1.0 + 2.0*(Float_t)column)*m_LogWidth;

        y = position;

        calLog = new THits(cal_name,cal_name,"CsI",m_LogWidth,m_CalHitLength,m_LogHeight);

        SetLogColor( calLog, eTotal );

        m_AllHitsX->cd();

        newnode = new TNNode(cal_name,cal_name,cal_name,x,y,z);

        newnode->setObject( calLog );
        newnode->SetParent( m_AllHitsX );
        m_CalNodeList->Add( newnode );

      }


      // !Build-up list of CAL Hits

      /*      m_CalHitList->Add( calLog );

      ((THits*)m_CalHitList->At(nhit))->SetX(x);

      ((THits*)m_CalHitList->At(nhit))->SetY(y);

      ((THits*)m_CalHitList->At(nhit))->SetZ(z+m_CalHeightOffset);

      ((THits*)m_CalHitList->At(nhit))->SetEn(eTotal*1000);

      ((THits*)m_CalHitList->At(nhit))->SetLayer(layer);

      ((THits*)m_CalHitList->At(nhit))->SetView(view);

      ((THits*)m_CalHitList->At(nhit))->SetElement(column);*/

      calLog->SetX(x);

      calLog->SetY(y);

      calLog->SetZ(z+m_CalHeightOffset);

      calLog->SetEn(eTotal*1000);

      calLog->SetLayer(layer);

      calLog->SetView(view);

      calLog->SetElement(column);

      m_CalHitList->Add( calLog );



      // printf("log %i calname %s, layer %i column %i energy %f\n",id,cal_name,layer,column,eTotal);


      nhit++;

    }

  }
  
  if ( m_recon ) drawCalRecon();

}

void CalDisplay::SetLogColor( THits *calLog, Float_t eTotal ) {
      

  // Define colour of logs based on the deposited energy

  if (eTotal > ( Float_t )m_cal_lev4) {

    calLog->SetLineColor(5);

  } else if (eTotal > ( Float_t )m_cal_lev3) {

    calLog->SetLineColor(2);

  } else if (eTotal > ( Float_t )m_cal_lev2) {

    calLog->SetLineColor(4); 

  } else if (eTotal > ( Float_t )m_cal_lev1) {

    calLog->SetLineColor(3); 

  }

}


void CalDisplay::setVisibility( Int_t vis ) {

  // Set visibility of tracker planes and hits.
  // Note: -2 = "node not drawn, sons not drawn", 2 = "node not drawn, sons drawn".

  switch( vis ) {

  case( kHideX ):
    m_AllHitsX->SetVisibility( -2 );
    m_AllHitsY->SetVisibility( 2 );
    break;

  case( kHideY ):
    m_AllHitsX->SetVisibility( 2 );
    m_AllHitsY->SetVisibility( -2 );
    break;

  case( kHideAll ):
    m_AllHitsX->SetVisibility( -2 );
    m_AllHitsY->SetVisibility( -2 );
    break;

  case( kHideNone ):
    m_AllHitsX->SetVisibility( 2 );
    m_AllHitsY->SetVisibility( 2 );
    break;

  }

}



Float_t CalDisplay::getLogEnergy( CalLog* hit, Int_t id, Int_t end, Bool_t isRealData ){

    Float_t eMax = -1000.; Int_t bestRange = -1; Int_t index;

    //Float_t gain[4] = {0.000065, 0.00022, 0.0037, 0.022};     //Berrie's gain (in Gev)
    Float_t gain[4] = {0.075, 0.3, 4.3, 25};       //Eric's gain Sept. 11 2001  - in MeV

    Int_t hitADC;                                        

    for (Int_t adc = 0; adc <= 1; adc++) {

       if (!end) hitADC = 4095 - hit->getAdc(hit->POS,adc);
       else hitADC = 4095 - hit->getAdc(hit->NEG,adc);

        if (hitADC > 3900) continue;

        if (hitADC > eMax) {

            eMax = hitADC;

            bestRange = adc;

        }       

    }

    if(!end) index = 2 * bestRange;

    else index = 2 * bestRange + 1;

    // Correct for unexpected ordering in ROOT ped file (Berrie's peds.)
    //if ( !end && bestRange==1 ) index = 1;
    //if ( end && !bestRange ) index = 2;

    TVector *vv;

    Float_t eHit;

    if ( bestRange > -1 ) { 
      
      if( isRealData ) { // raw real data !!!! 
        
        if ( eMax && m_peds ) {
          
          vv = (TVector *)m_peds->At(index);
          
          eHit = eMax-vv->operator()(id);
          
        }
        
        else eHit = 0;
      }
      
      else eHit = eMax; // MC data !!!

      // Otherwise just multiply by gain
      
      //Float_t eHit = eMax*gains[bestRange];
      
      return eHit*gain[ bestRange ];

    }
    
    return 0.0;
    
}




Char_t *CalDisplay::getCalName(Int_t ihit){

  Char_t cs[40];

   TString *CalName = new TString("CAL_log");

   sprintf(cs,"%d",ihit);

   CalName->Append(cs);   

   return (Char_t*)CalName->Data();

}



void CalDisplay::drawCalRecon() {

  Double_t x, y, z, anglex, angley,anglez;

  //  ReconHeader* recHeader;

  //  m_recon->Dump();

  //  recHeader = m_recon->getReconFlags(); //check way when switch back from analysis to display 
  //  if uncommented this line gives bad pointer to class recHeader.

  // Indeed this is a problem!  recHeader part re-commented by Nick, 28-09-2001 for compiled version.

//    if  ( recHeader->HasCalNoEnergyRecon()) {
//      cout << " ===> No energy reconstruction in CAL !! " << endl;
//      return;
//    }
 
//    m_Gammas = 0;

//    if ( recHeader->HasGamma() ) {
//      cout << "\n*********** GAMMA!!! ************\n";
//      m_Gammas = 1;
//    }
  
  CalRecon* calRecon = m_recon->getCalRecon();


  TObjArray *m_CalClusterList = calRecon->getCalClusters();
  
  Int_t numClus = m_CalClusterList->GetEntries();
  
 
  
 // Create a new list of reconstructed hits.

  //if ( m_CalReconHitList != 0 ) delete m_CalReconHitList;
  
  //m_CalReconHitList = new TObjArray();

  cout << " Num. of clusters in CAL: " << numClus << endl; 

  for ( Int_t iclus=0; iclus<numClus; ++iclus ) {
   
   // Loop over all clusters for this event.

   CalCluster *cluster = (CalCluster*) m_CalClusterList->At(iclus);

   m_eTot = (cluster->getEsum())*pow(10,-9);

   m_eCorr = (cluster->getEcorr())*pow(10,-9);

   m_eLeak = (cluster->getEleak())*pow(10,-9);
  
   cout << "Etot:" << m_eLeak << endl;
   
   const TArrayF* m_eneLayer = &cluster->getEneLayer();

   Int_t nHitsE = m_eneLayer->GetSize();
   
   const TArrayF* m_xposLayer = &cluster->getXLayer();

   //Int_t nhitsX = m_xposLayer->GetSize();
   
   const TArrayF* m_yposLayer = &cluster->getYLayer();

   //Int_t nhitsY = m_yposLayer->GetSize();


   for ( Int_t ihit=0; ihit<nHitsE; ++ihit ) {
   
     Float_t xpos = m_xposLayer->At(ihit);
     Float_t ypos = m_yposLayer->At(ihit);
     Float_t energy = m_eneLayer->At(ihit);
     printf( " lay: %d -- X: %f   Y: %f  Energy: %f \n",ihit,xpos,ypos,energy);

   }
   TVector3* position = cluster->getPosition();

      x = position->X();

      y = position->Y();

      z = position->Z();

   TVector3* direction = cluster->getDirection();
   

      anglex = direction->X();

      angley = direction->Y();

      anglez = direction->Z();
      
      cout << " X:" << x << " anglex: " << anglex << "  Y:" << y << " angley:" 
           << angley << " Z:" << z << " anglez; " << anglez << endl; 
  }

}

---