MyTree.C

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2010 by Olivier Stezowski                               *
 *   stezow(AT)ipnl.in2p3.fr                                               *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/


#include "MyTree.h"
#include "TSystem.h"

ClassImp(DoMyKeyTree);

DoMyKeyTree::DoMyKeyTree(const char *name, const char *title, TTree *tree) : 
        Watcher(name,title), TTreeBuilder(name,title,tree)
{
        if ( GetTree() ) { // just add the branch for that watcher
                SetBranches();
        }
}

Bool_t DoMyKeyTree::SetTrigger( ADF::DFTrigger *trigger ) 
{
        // Action once a the Frame is passed to this 
        fFrame = trigger->GetInputSharedFP();
        
        // unset the frame to be treated
        if ( fFrame == 0x0 ) 
                return false;
        
        // trigger and Frame set,
        fKey = dynamic_cast<const AgataKey *>(fFrame->GetFrame()->GetKey());
        
        if ( fKey == 0x0 ) 
                { return false; }
        else
                { return true; }
}

void DoMyKeyTree::SetBranches()
{
        GetTree()->Branch("EventNumber", &fEventNumber, "evnumber/i");
        GetTree()->Branch("TimeStamp", &fTimestamp,   "timestamp/l");
}

void DoMyKeyTree::Exec(Option_t * /*option*/) 
{
        // for debugging purpose
        if ( gDebug > 3 ) 
                printf("DoMyKeyTree::Exec is called \n");
        
        // be sure the frame has been set properly
        if ( fFrame == 0x0 || !fFrame->IsValid() ) 
                { SetLastExecStatus(1); return; }
        
        // get event # and timestamp from the key and fill the tree
        fEventNumber = fKey->GetEventNumber();
        fTimestamp   = fKey->GetTimeStamp();
        
        // fill the tree only if this is the owner of the tree ... otherwise by TreeMaster
        FillTree();
}

ClassImp(DoMyRawAncillaryTree);

DoMyRawAncillaryTree::DoMyRawAncillaryTree(const char *name, const char *title, TTree *tree) : 
        RancLegnaroWatcher(name,title), TTreeBuilder(name,title,tree)
{
        if ( GetTree() ) { // just add the branch for that watcher
                SetBranches();
        }
}

void DoMyRawAncillaryTree::SetBranches()
{
        GetTree()->Branch("lenrawBuf", &lenrawBuf, "lenrawBuf/I");
        GetTree()->Branch("rawBuf",     rawBuf,    "rawBuf[lenrawBuf]/F");
        
        GetTree()->Branch("TSanc",     &fTimestamp,    "TSanc/l");
        
}

DoMyRawAncillaryTree::~DoMyRawAncillaryTree()
{
}

void DoMyRawAncillaryTree::Exec(Option_t *option) 
{
        if ( gDebug > 3 ) 
                printf("DoMyRawAncillaryTree::Exec is called \n");
        
        // be sure the frame has been set properly
        if ( fFrame == 0x0 || !fFrame->IsValid() ) 
                { SetLastExecStatus(1); return; }
        
        fTimestamp = ((AgataKey *)fFrame->GetFrame()->GetKey())->GetTimeStamp();
        
        // filling the raw buffer is done in the base class 
        RancLegnaroWatcher::Exec(option);
        // fill the tree only if this is the owner of the tree ... otherwise by TreeMaster
        FillTree();
}

ClassImp(DoMyAncillaryTree);

DoMyAncillaryTree::DoMyAncillaryTree(const char *name, const char *title, TTree *tree) : 
        BaseDanteWatcher(name,title), TTreeBuilder(name,title,tree)
{
        if ( GetTree() ) { // just add the branch for that watcher
                SetBranches();

        }
}

DoMyAncillaryTree::~DoMyAncillaryTree()
{
}

void DoMyAncillaryTree::SetBranches()
{
         GetTree()->Branch("DR",&DR,"DR/F");
         GetTree()->Branch("DTheta",&DTheta,"DTheta/F");
         GetTree()->Branch("DPhi",&DPhi,"DPhi/F"); 
}

void DoMyAncillaryTree::Exec(Option_t *option) 
{
        if ( gDebug > 3 ) 
                printf("DoAncillaryTree::Exec is called \n");
        
        // be sure the frame has been set properly
        if ( fFrame == 0x0 || !fFrame->IsValid() ) 
                { SetLastExecStatus(1); return; }
        
        // filling the raw buffer is done in the base class 
        RancLegnaroWatcher::Exec(option);       
        
        // Compute from the Recoil velocity from Dante
        /*
         TVector3 vtmp(fDante[i].Xr,fDante[i].Yr,fDante[i].Zr); 
         DR = vtmp.Mag();
         DTheta = vtmp.Theta();
         DPhi = vtmp.Phi();     */
        // fill the tree only if this is the owner of the tree ... otherwise by TreeMaster
        FillTree();
}

ClassImp(DoMyTrackTree);

DoMyTrackTree::DoMyTrackTree(const char *name, const char *title, TTree *tree) : 
        TrackedWatcher(name,title), 
        TTreeBuilder(name,title,tree),
        fVertexBuilder(0x0)
{
        // Vertex used for doppler correction, last registered
        fVertexBuilder = VertexBuilder::theCurrentVertexBuilder();
        
        if ( GetTree() ) { // just add the branch for that watcher
                SetBranches();
        }
}

void DoMyTrackTree::SetBranches()
{
        GetTree()->Branch("number_of_gammas",&number_of_gammas,"number_of_gammas/I");
        GetTree()->Branch("gammaE", gammaE, "gammaE[number_of_gammas]/F");
        GetTree()->Branch("gammaEDC", gammaEDC, "gammaEDC[number_of_gammas]/F");        
        GetTree()->Branch("gammaX1",gammaX1,"gammaX1[number_of_gammas]/F");
        GetTree()->Branch("gammaY1",gammaY1,"gammaY1[number_of_gammas]/F");
        GetTree()->Branch("gammaZ1",gammaZ1,"gammaZ1[number_of_gammas]/F");
        GetTree()->Branch("gammaX2",gammaX2,"gammaX2[number_of_gammas]/F");
        GetTree()->Branch("gammaY2",gammaY2,"gammaY2[number_of_gammas]/F");
        GetTree()->Branch("gammaZ2",gammaZ2,"gammaZ2[number_of_gammas]/F"); 
        
        GetTree()->Branch("TStrack",     &fTimestamp,    "TStrack/l");
}

DoMyTrackTree::~DoMyTrackTree()
{
}

void DoMyTrackTree::Exec(Option_t * /*option*/) 
{       
        // be sure the frame has been set properly
        if ( fFrame == 0x0 || !fFrame->IsValid() ) 
                { SetLastExecStatus(1); return; }
        
        fTimestamp = ((AgataKey *)fFrame->GetFrame()->GetKey())->GetTimeStamp();

        // to get the data part of the frame
        const GammaTrackedInterface *data = 
                GetCstDataPointer<GammaTrackedInterface>(fFrame); 
        VertexInterface *vertex = fVertexBuilder->GetVertex();
        
        // now fill the tree
        number_of_gammas = data->GetNbGamma();
        for (UShort_t i = 0u; i < number_of_gammas; i++) {
                
                const TrackedHit *gamma1 = data->GetGamma(i); 
                
                gammaE[i]    = gamma1->GetE();
                gammaEDC[i]  = DoDopplerCorrection(gamma1, vertex);
                
                gammaX1[i] = gamma1->GetX();
                gammaY1[i] = gamma1->GetY();
                gammaZ1[i] = gamma1->GetZ();
                
                if ( gamma1->GetHit(1) ) {
                        gammaX2[i] = gamma1->GetHit(1)->GetX();
                        gammaY2[i] = gamma1->GetHit(1)->GetY();
                        gammaZ2[i] = gamma1->GetHit(1)->GetZ(); 
                }
                else {
                        gammaX2[i] = gamma1->GetX();
                        gammaY2[i] = gamma1->GetY();
                        gammaZ2[i] = gamma1->GetZ();                    
                }
        }
        FillTree();
}


ClassImp(DoMyHitTree);

DoMyHitTree::DoMyHitTree(const char *name, const char *title, TTree *tree) : 
  EventPSAWatcher(name,title),
  TTreeBuilder(name,title,tree),
        fVertexBuilder(0x0),
        fTrans(0x0)
{
        fVertexBuilder = VertexBuilder::theCurrentVertexBuilder();
        if ( fVertexBuilder == 0x0 )
                std::cout << "Should not be NULL " << std::endl;

  fTrans = new AgataGeometryTransformer();
        if ( gSystem->AccessPathName("CrystalPositionLookUpTable") ) { // does not exist
                if ( gSystem->AccessPathName("CrystalPositionLookUpTable.dat") ) {
                        fLog << error << " CrystalPositionLookUpTable required " << nline; 
                }
                else {
                        fLog << info << " Ge positions extracted from CrystalPositionLookUpTable.dat " << nline; 
                        fTrans->ReadTransformations("CrystalPositionLookUpTable.dat");
                }
        }
        else {
                fLog << info << " Ge positions extracted from CrystalPositionLookUpTable.dat " << nline; 
                fTrans->ReadTransformations("CrystalPositionLookUpTable");
        }

        if ( GetTree() ) { // just add the branch for that watcher
        SetBranches();
        }

          // create spetra
  fIdSpectra = MakeTH1<TH1F>("hit_ID_spectra","Spectra of the ID ",100,0,100);
  fSumSpectra = MakeTH1<TH1F>("SumSpectra","Sum of all core energy ",4096,0,4096);
        fhitAgata  = MakeTH3<TH3F>("hitDW","Tracking from TaDSSD to WallD ;x;y;z",100,-200.,200.,100,-200.,200.,200,-4000.,4000.);


}

DoMyHitTree::~DoMyHitTree()
{
}

void DoMyHitTree::SetBranches()
{
         GetTree()->Branch("number_of_hits",&number_of_hits,"number_of_hits/I");
         GetTree()->Branch("mult_gamma",&mult_gamma,"mult_gamma/I");
         GetTree()->Branch("hitE",hitE,"hitE[number_of_hits]/F");
         GetTree()->Branch("hitX",hitX,"hitX[number_of_hits]/F");
         GetTree()->Branch("hitY",hitY,"hitY[number_of_hits]/F");
         GetTree()->Branch("hitZ",hitZ,"hitZ[number_of_hits]/F");
         GetTree()->Branch("hit_per_cryst",hit_per_cryst,"hit_per_cryst[number_of_hits]/I"); 
         GetTree()->Branch("hitId",hitId,"hitId[number_of_hits]/I");
         GetTree()->Branch("hitSg",hitSg,"hitSg[number_of_hits]/I"); 
   GetTree()->Branch("coreE0",coreE0,"coreE0[number_of_hits]/F");
   GetTree()->Branch("coreE1",coreE1,"coreE1[number_of_hits]/F");
 
         GetTree()->Branch("MhitX",MhitX,"MhitX[number_of_hits]/F");
         GetTree()->Branch("MhitY",MhitY,"MhitY[number_of_hits]/F");
         GetTree()->Branch("MhitZ",MhitZ,"MhitZ[number_of_hits]/F");
   GetTree()->Branch("velDop",velDop,"velDop[number_of_hits]/F");
   GetTree()->Branch("coreDE0",coreDE0,"coreDE0[number_of_hits]/F");
   GetTree()->Branch("coreDE1",coreDE1,"coreDE1[number_of_hits]/F");
        
         GetTree()->Branch("TSHit",     &fTimestamp,    "TSHit/l");
}

void DoMyHitTree::Exec(Option_t * /*option*/) 
{
//  cout << "Exec MyhitTree " << endl;
        // be sure the frame has been set properly
        if ( fFrame == 0x0 || !fFrame->IsValid() ) 
                { cout << "non valid " << endl; SetLastExecStatus(1u); return; }
  if ( !fFrame->IsValid()){  //For empty event:data:psa purpose
     number_of_hits =0;
     mult_gamma =0;
     for(int i=0; i<MaxHits ; i++) {//empty array coreE0 coreE1
       coreE0[i]=-1;
       coreE1[i]=-1;
       velDop[i]=-1;
       coreDE0[i]=-1;
       coreDE1[i]=-1;
    }
    return;
  } 

        fTimestamp = ((AgataKey *)fFrame->GetFrame()->GetKey())->GetTimeStamp();

  Double_t max_e = 0.0;
        Double_t xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal, xt,yt,zt; TrackedHitImp<EHit,EHit> gamma;
        // the vertex used for doppler correction
        VertexInterface *vertex = 0x0;
        vertex = fVertexBuilder->GetVertex();

 
  //Multiplicity gamma in case if more than one frame per event
  mult_gamma=0;

  // Get beta0 for doppler correction
  Float_t beta0 =-1;
  // Nb_of hits per crystals
  Int_t nb_hits =-1;
  if ( vertex ) {
    beta0=vertex->GetBeta();
  }

  //Initialisation:
  //--Init total number of hit per event
  number_of_hits = 0;

        for(UInt_t i = 0u; i < GetNbFramePSA(); ++i) {
    mult_gamma++;
                
    PSAInterface *data = GetDataPSA(i); 
    max_e = 0.0;
    nb_hits =0;

    //Fill spectra
    fIdSpectra->Fill(data->GetUID());
    fSumSpectra->Fill(data->GetE());

      

   // First loop on the hit
   // to get over the hits with the max of the energy, and set the energy to the core
   // and the number of hits per crystal
                for (Int_t j = 0; j < data->GetNbHits(); ++j) {
                        const PSAHit* pHit = (PSAHit*)data->GetHit(j); 

      // Now look for the hit with the highest energy, taken for Doppler correction
     // Kind of tracking in crystal
     //
      //Get hit position
                        pHit->GetXYZ(xLocal, yLocal, zLocal);
      //To the AGATA global ref
                        fTrans->Local2Global(data->GetUID(), xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal);


                        if ( pHit->GetE() > max_e ) {
                                gamma.SetXYZ(-yGlobal, -xGlobal, zGlobal); //x,y,z->-y,-x,-z
                                max_e = pHit->GetE();
                        }

    // High gain core as gamma energy   
                gamma.SetE(data->GetE());

    //Increment number of hit per crystal:
    nb_hits++;
    }


                //Second loop on the hit, for all treatment and gamma
                for (Int_t j = 0; j < data->GetNbHits(); ++j) {
                        
                        const PSAHit* pHit = (PSAHit*)data->GetHit(j); 
                        
                        hitE[number_of_hits]   = pHit->GetE();
                        hitX[number_of_hits]   = pHit->GetX();
                        hitY[number_of_hits]   = pHit->GetY();
                        hitZ[number_of_hits]   = pHit->GetZ();
      hit_per_cryst[number_of_hits] = nb_hits;
                        hitId[number_of_hits]  = data->GetUID() ;
                        hitSg[number_of_hits]  = pHit->GetID();
                  coreE0[number_of_hits] = data->GetE();
      coreE1[number_of_hits] = data->GetE(1u);

      //Get hit position
                        pHit->GetXYZ(xLocal, yLocal, zLocal);
      //To the AGATA global ref
                        fTrans->Local2Global(data->GetUID(), xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal);

      // Set the Position of the hit in global ref of PRESPEC
      MhitX[number_of_hits]= yGlobal; // inversion X-Y and opposite direction
      MhitY[number_of_hits]= xGlobal; // inversion Y-X and opposite direction
      MhitZ[number_of_hits]= zGlobal; // z
      fhitAgata->Fill(yGlobal,xGlobal, zGlobal);



    if ( vertex &&beta0>0 && beta0<1) {
    
      // Get the proper LyccaVertec pointer
      vertex->GetPosition(xt,yt,zt);
 
      fhitAgata->Fill(xt,yt,zt);
    }

    //High gain core
        gamma.SetE(data->GetE());
    if ( vertex ) {
       Double_t edc = DoDopplerCorrection(&gamma,vertex);
       coreDE0[number_of_hits] = edc;            
       velDop[number_of_hits] = beta0;
    }
    else coreDE1[number_of_hits] = -1; //if there is no vertex, no energy


    // Low gain core        
                gamma.SetE(data->GetE(1u));
    if ( vertex ) {
       Double_t edc = DoDopplerCorrection(&gamma,vertex);
       coreDE1[number_of_hits] = edc;            
    }
    else coreDE1[number_of_hits] = -1; //if there is no vertex, no energy



    number_of_hits++;

        if ( number_of_hits == 500 ) {
                SetLastExecStatus(1u);
      cout << "merdaasse" <<endl;
                break;
    }
    

        } // j, hit loop

  
    if ( number_of_hits == 500 )
                fLog << error << " Number of Hits to hit in ! " << number_of_hits << nline;

    
        } // end loop over PSA frame
        FillTree();
  
}


Double_t DoMyHitTree::DoDopplerCorrection(const TrackedHit *hit, VertexInterface *vertex)
{
        Double_t e = hit->GetE(), beta = vertex->GetBeta(), gamma = vertex->GetGamma();
        
        Double_t px, py, pz, dx, dy, dz;
        vertex->GetPosition (px, py, pz);
        vertex->GetDirection(dx, dy, dz); 
    
//    std::cout << " " << px << " " << py << " " << pz <<std::endl;
//    std::cout << " " << dx << " " << dy << " " << dz <<std::endl;

        Double_t x = hit->GetX() - px;
        Double_t y = hit->GetY() - py;
        Double_t z = hit->GetZ() - pz;
        
        Double_t dd = x*x + y*y + z*z;
        if( !(dd == 0.0) ) {
                Double_t cosTheta = (x*dx + y*dy + z*dz)/::sqrt(dd);
                e = VertexInterface::DopplerCorrection(e,beta,gamma,cosTheta);
        }
        else e = -1.0;
        
        return e;
}

void DoMyHitTree::FillSpectraByUser()
{
   //Nothing yet.
}



ClassImp(DoMyCoincTree);

DoMyCoincTree::DoMyCoincTree(const char *name, const char *title, TTree *tree)
: EventWatcher(name,title), TTreeBuilder(name,title,tree)
{
        if ( GetTree() ) { // just add the branch for that watcher
                SetBranches();
        }
}

DoMyCoincTree::~DoMyCoincTree()
{
}

void DoMyCoincTree::SetBranches()
{
        // cout << " SetBranches\t"<< fTSCoinc << endl;
        GetTree()->Branch("TSDiff", &fTSCoinc, "TSDiff/F");
}

void DoMyCoincTree::Exec(Option_t *option) 
{
        // to get the main key
        AgataKey * mainkey = dynamic_cast<AgataKey *> (GetTrigger()->GetInputFrame(0u)->GetKey()), *subkey = 0x0;
        
        fTSCoinc = 0;
                
        // loop on subframes and compute the timestamp differences
        for(UInt_t which = 1u; which < GetTrigger()->GetNbInputFrame(); which++){
                
                SharedFP *fp = GetTrigger()->GetInputSharedFP(which);
                if ( fp == 0x0 || !fp->IsValid() )
                        continue;

                // the subkey
                subkey = dynamic_cast<AgataKey *> (GetTrigger()->GetInputFrame(which)->GetKey());
                break;
        }
        if ( subkey ) {
                fTSCoinc = subkey->GetTimeStamp() - mainkey->GetTimeStamp();
        }
   
   // fill the tree only if this is the owner of the tree ... otherwise by TreeMaster
   FillTree();
}