GSI/TrackedWatchers.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.             *
 ***************************************************************************/

// Watcher definition
#ifndef _TrackedWatchers
#include "TrackedWatchers.h"
#endif

// #include "GenericFrame.h"
#include "AgataKeyFactory.h"

#include "MetaWatchers.h"
#include "TVector3.h"

#ifndef ROOT_TGeoMatrix
#include "TGeoMatrix.h"
#endif

#ifndef Gw_AgataGeometryTransformer
#include "AgataGeometryTransformer.h"
#endif

using namespace ADF;

ClassImp(TrackedWatcher);

Bool_t TrackedWatcher::SetTrigger( ADF::DFTrigger *trigger ) 
{
        Bool_t ok = Watcher::GetFromTrigger(trigger, "data:tracked", fFrame);
        return ok;
}

void TrackedWatcher::SetAgataOffset(Double_t offset)
{
        fZOffset = offset;
}

Double_t TrackedWatcher::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);
        
        Double_t x = hit->GetX() - px;
        Double_t y = hit->GetY() - py;
        Double_t z = hit->GetZ() + fZOffset - 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;
}

ClassImp(Coinc2D);

Coinc2D::Coinc2D(const char *name, const char *title, TDirectory *sp_dir, TDirectory *tag_dir):
        TrackedWatcher(name,title,sp_dir,tag_dir),
        fVertexBuilder(0x0),
    fTrackedEnergy(0x0),
    fTrackedEnergyDopler(0x0),
        fFold(0x0),
        fGxG(0x0),
        fGxG_DC(0x0),
        fPos3D(0x0)
{
        // Vertex used for doppler correction, last registered
        fVertexBuilder = VertexBuilder::theCurrentVertexBuilder();
        
        // Now spectra
    fTrackedEnergy = MakeTH1<TH1F>("TrackedEnergy","Tracked Energy",6000,0,6000);
        
    fTrackedEnergyDopler = MakeTH1<TH1F>("TrackedEnergyDopler","Tracked Energy with Doppler correction",6000,0,6000);
    TagOn(fTrackedEnergyDopler);
        
        fFold    = MakeTH1<TH1F>("Fold","Fold distributions",10,0,10);
        TagOn(fFold);
                
        fGxG     = MakeTH2<TH2F>("GxG","Gamma Gamma matrix",2000,0,1000,2000,0,1000);
        fGxG_DC  = MakeTH2<TH2F>("GxG_DC","Gamma Gamma matrix with Doppler Correction",2000,0,1000,2000,0,1000);

    fPos3D =  MakeTH3<TH3F>("Pos3D","Positions of Tracked Gamma-rays ;x;y;z",200,-200,200,200,-200,200,200,-400,0);
}

Coinc2D::~Coinc2D()
{
// no need to delete the histograms, done because they have been added to the pool
}

void Coinc2D::DoCanvas(TCanvas *c, Option_t *o)
{
        TString opt = o;
        
        if ( opt.Contains("") ) {
                c->Divide(1,2,0.001,0.001);
                c->cd(1);
        fTrackedEnergy->Draw(o);
                TPad *pad = new TPad("GxG_FOLD", "Pad for Fold distribution",0.7,0.6,1,1);
                pad->Draw();
                pad->cd();
                fFold->Draw();
                c->cd(2);
        fTrackedEnergyDopler->Draw(o);
        }
}

void Coinc2D::Exec(Option_t */*option*/) 
{
// Frame by Frame action
        // be sure the frame has been set properly
        if ( fFrame == 0x0 || !fFrame->IsValid() ) 
                { SetLastExecStatus(1); return; }
//      printf("... is called \n");
                
        if ( gDebug > 3 ) 
                printf("Coinc2D::::Exec is called \n");
                
        // to get the data part of the frame
        const GammaTrackedInterface *data =     
                GetCstDataPointer<GammaTrackedInterface>(fFrame); 
        VertexInterface *vertex = fVertexBuilder->GetVertex();
  //  vertex->SetDirection(0.,-1.,1.);
  //  vertex->SetPosition(0.,0.,-50.);
  //  vertex->SetBeta(0.02);
        
        // now fill histograms
        UInt_t fold = data->GetNbGamma();
        
        fFold->Fill( fold );
        for (UShort_t i = 0u; i < fold; i++) {
                
                const TrackedHit *gamma1 = data->GetGamma(i); 
                
                fPos3D->Fill(gamma1->GetX(),gamma1->GetY(),gamma1->GetZ());
//              cout << gamma1->GetHit(1)->GetX() << " " << gamma1->GetHit(1)->GetY() << " " << gamma1->GetHit(1)->GetZ() << endl;

                Double_t edc1 = DoDopplerCorrection(gamma1, vertex);
                Double_t e1   = gamma1->GetE();
                
        fTrackedEnergy->Fill(e1);
                if ( vertex->GetBeta() > 0.001 )  
            fTrackedEnergyDopler->Fill(edc1);
                        
                for (UShort_t j = i+1u; j < fold; j++) {
                        
                        const TrackedHit *gamma2 = data->GetGamma(j); 
                        
                        Double_t edc2 = DoDopplerCorrection(gamma2, vertex);    
                        Double_t e2   = gamma2->GetE();
                        
                        fGxG->Fill(e1,e2);
                        fGxG->Fill(e2,e1);      
                        if ( vertex->GetBeta() > 0.001 )  {
                                fGxG_DC->Fill(edc1,edc2);
                                fGxG_DC->Fill(edc2,edc1);
                        }                               
                }
        }
}

//-------------------------------------------------------------------------------------

ClassImp(EveTrack);

Float_t EveTrack::fgMetric = 10.;

EveTrack::EveTrack(const char *name, const char *title, TDirectory *sp_dir, TDirectory *tag_dir): 
        TrackedWatcher(name, title, sp_dir, tag_dir),
        fMaxCycle(100),
        fPathFile(Gw::AgataEventDisplay::GetDefaultAgataPath())
{
  AgataEventDisplay::Instance()->SetTrackStyle("Cone");

  fMeanRadius = MakeTH1<TH1F>("MeanRadius","Mean radius of the tracks (cm)", 500, -1, 34);  
  fFold  = MakeTH1<TH1I>("Fold","Number of Tracks per event", 22, -1, 20) ;
  fHitEnergy = MakeTH1<TH1F>("HitEnergy","Energy of hits", 1000,-1,1000.) ;
  fTrackEnergy   = MakeTH1<TH1F>("TrackEnergy","Energy of tracks", 2000,0.,2000.) ;
}

EveTrack::~EveTrack()
{ 

}

void EveTrack::SetActive(Bool_t active) 
{
  TTask::SetActive(active); 
  AgataEventDisplay::Instance()->Reset();
} 

void EveTrack::DoCanvas(TCanvas* c, Option_t* /*o*/)
{
  c->Divide(2,2,0.001,0.001);
  
  c->cd(1); 
  fMeanRadius->Draw();  
  
  c->cd(2); 
  fFold->Draw();
  
  c->cd(3); 
  fHitEnergy->Draw();   
  
  c->cd(4); 
  fTrackEnergy->Draw(); 
}

void EveTrack::ShowEve(const Char_t* agataPathFile)
{
   AgataEventDisplay::SetDefaultAgataPath(agataPathFile);
   AgataEventDisplay::Instance()->BuildDefaultGeometry(false, false);
   AgataEventDisplay::Instance()->ShowDisplay();
}

void EveTrack::Exec(Option_t */*option*/) 
{
  // Frame by Frame action
  
  // be sure the frame has been set properly
        // be sure the frame has been set properly
        if ( fFrame == 0x0 || !fFrame->IsValid() ) 
                { SetLastExecStatus(1); return; }
  
  if ( gDebug > 3 ) 
    printf("EveTrack::::Exec is called \n");
  
  // to get the data part of the frame
  const GammaTrackedInterface *data = GetCstDataPointer<GammaTrackedInterface>(fFrame); 

  // now fill display container
  UInt_t fold = data->GetNbGamma();
  
  if (fold == 0)
    return;
  
  fFold->Fill(fold);

  AgataEventContainer* agataContainer = AgataEventDisplay::Instance()->GetEventContainer(); 
  
  if (agataContainer->GetNofEvents() > fMaxCycle) {
    TTask::SetActive(false);
  }
  
  if ( gDebug > 3 )   
          std::cout << "TrackHits " << std::endl;
  
  if ( gDebug > 3 ) 
          std::cout << "Fold " << fold << std::endl;

  for (UShort_t i = 0u; i < fold; ++i) {
        
    const TrackedHit* gamma = (TrackedHit*)data->GetGamma(i); 
    TrackHit* hit = agataContainer->NewTrackHit();
    fTrackEnergy->Fill(gamma->GetE());
    
    hit->SetX(gamma->GetX()/fgMetric);
    hit->SetY(gamma->GetY()/fgMetric);
    hit->SetZ(gamma->GetZ()/fgMetric);
    hit->SetE(gamma->GetE());
    
    for (Int_t j = 0; j < gamma->GetNbHits(); ++j) {
      Hit* gammaHit = gamma->GetHit(j);
      StdHit* shit = hit->NewHit();

      Float_t radius = TMath::Sqrt(gammaHit->GetX()*gammaHit->GetX() + gammaHit->GetY()*gammaHit->GetY() + 
                                   gammaHit->GetZ()*gammaHit->GetZ());
      fMeanRadius->Fill(radius/fgMetric);
      fHitEnergy->Fill(gammaHit->GetE());
         if ( gDebug > 3 ) {
                 std::cout << "TrackHit Pos " << gammaHit->GetX()/fgMetric << " " <<  gammaHit->GetY()/fgMetric 
                 <<  " " << gammaHit->GetZ()/fgMetric << " " << gammaHit->GetE() << std::endl;
          }

      shit->SetX(gammaHit->GetX()/fgMetric);
      shit->SetY(gammaHit->GetY()/fgMetric);
      shit->SetZ(gammaHit->GetZ()/fgMetric);
      shit->SetE(gammaHit->GetE());
    }
  }
  agataContainer->FillTracks();
}