PrismaWatchers.C

Go to the documentation of this file.

#include "PrismaWatchers.h"
#include "MetaWatchers.h"
#include "AgataGeometryTransformer.h"

#include "TMath.h"

#include <iomanip>

//includes for libPrisma.so or any other library needed
#include "RandomGenerator.hh"

namespace  {
        const int   pM = 0; //monitor
        const int   pA = 2; //MCP:x,y 
        const int   pR = 6; //MWPPC:tof,xleft,xright,xcathode,y 
        const int   pI = 4; //4 IC sections:A,B,C,D
        const int   pS = 4; //4 Side sections:A,B,C,D
}

ClassImp(BasePrismaWatcher);

BasePrismaWatcher::BasePrismaWatcher(const char *name, const char *title): 
        RancLegnaroWatcher(name,title), 
        fPM(0x0),
        fCurrentSize(0),
        fDetData()
{
} 

BasePrismaWatcher::~BasePrismaWatcher()
{
        if ( fPM )
                delete fPM;
        
        for (size_t i = 0; i < fDetData.size(); i++) {
                if ( fDetData [i] ) {
                        delete fDetData[i]; fDetData[i] = 0x0; 
                }
        }
}

void BasePrismaWatcher::Configure(const char *pmanager, const char *plut) 
{
        if ( fPM )
        { delete fPM; fPM = 0x0; }
        
        if ( InitLUT(plut) ) {
                fPM = new prismaManager(pmanager); Random::randomize();
        }
}

Bool_t BasePrismaWatcher::InitLUT(const char *filename)
{
        FILE * conf_file;
        if( (conf_file = fopen(filename,"r")) == 0x0) {
                std::cout << " Error, cannot open " << filename << std::endl; return false;
        }
        
        for (Int_t ii = 0; ii< gMaxSlotsForPrisma*gMaxChannels; ii++) {
                Int_t Conv, Ch, Cont;
                if(fscanf(conf_file,"%d %d %d",&Conv,&Ch,&Cont) != 3) {
                        cout << "Error reading line " << ii << " of " << filename << endl;
                        fclose(conf_file);
                        
                        return false;
                }
                LUT[Conv][Ch]= Cont;
                std::cout << "LUT[ " << Conv << " ][ " << std::setw(2) <<  Ch << " ] = " << std::setw(3) << LUT[Conv][Ch] << std::endl;
        }       
        
        std::cout << "Prisma LUT " << filename << " is correctl initialised " << std::endl; return true;
}

detData *BasePrismaWatcher::GetNextDetData(Int_t ndet, Int_t npar)
{
        detData *data;
        
        // stack not large enough, so add one element
        if ( fCurrentSize == fDetData.size() ) {
                data = new detData (ndet,npar);
                fDetData.push_back( data );
        }
        else { data = fDetData[fCurrentSize++]; data->reset_detData(ndet,npar); }
        
        return data;
}

void BasePrismaWatcher::Exec(Option_t * option)
{       
        // cannot extract data, manager is not set properly
        if ( fPM == 0x0 )
                return;
        
        // event
        Int_t numDet[gmaxTypes]; detData* aData;
        
        // to load data fron ancillary frame to rawBuf
        SetLastExecStatus(0);
        RancLegnaroWatcher::Exec(option);
        if ( GetLastExecStatus() != 0 ) 
                return;
        
        // fill the map with the current data obtained using the LUT and rawBuf
        ::memset(theMap, 0, sizeof(theMap));
        
        for (int ii= 0; ii < gMaxSlotsForPrisma; ii++) { 
                for (int jj=0; jj < gMaxChannels; jj++) {
                        //std::cout << " LUT[" << ii << "]["<< jj<<"]= " <<  LUT[ii][jj] << std::endl;
                        if ( LUT[ii][jj] < 0) 
                                continue;
                        theMap[ LUT[ii][jj] ] =  rawBuf[ii*gMaxChannels+jj]; // Converter[ii].GetValue(jj);
                }
        }
        
        // count the number of data (parameter) in the event for each type of detector
        int nParam[gmaxTypes][gMaxSections]; ::memset(nParam,0,sizeof(int)*gmaxTypes*gMaxSections);
        
        for (int zz=0; zz<gmaxTypes; zz++) {
                numDet[zz]=0;
                for (int jj=0; jj<gMaxSections; jj++){
                        for (int ii=0; ii<6; ii++) {
                                if ( theMap[60*zz+6*jj+ii] > 250.f  ){
                                        nParam[zz][jj]++;
                                }
                        }
                }
        }
        
        for (int zz=0; zz<gmaxTypes; zz++){
                for (int jj=0; jj<gMaxSections; jj++){
                        if( nParam[zz][jj] > 0 )
                                numDet[zz]++;
                }
        }
        
        // Reset stack of data
        ResetDetDataStack();
        
        
        theMcpData.clear();
        
        for (int jj=0; jj<gMaxSections; jj++ ) {
                if( nParam[0][jj] < 1 ) 
                        continue;
                aData = GetNextDetData(jj,pA); theMcpData.push_back( aData );
                for(int ii=0; ii<pA; ii++ )
                        if (theMap[6*jj+ii] > 0)
                                aData->set( ii, theMap[6*jj+ii] + Random::uniform() - 0.5 );
                        else 
                                aData->set( ii, 0 );                            
        }       
        
        
        thePpacData.clear();
        
        //  for(int jj=0; jj<numDet[1]; jj++ ) {
        for(int jj=0; jj<gMaxSections; jj++ ) {
                if( nParam[1][jj] < 1 ) 
                        continue;
                aData = GetNextDetData(jj,pR); thePpacData.push_back( aData );
                for(int ii=0; ii<pR; ii++ )
                        if ( theMap[60 + 6*jj+ii] > 0 )
                                aData->set( ii, theMap[60 + 6*jj+ii] + Random::uniform() - 0.5 );
                        else 
                                aData->set( ii, 0 );    
        }
        
        /* for(int jj=0; jj<numDet[1]; jj++ ) {
         aData = thePpacData[jj];
         std::cout << "numm = " << numDet[1] << " " << nParam[1][aData->num_det()] << std::endl;
    for(int ii=0; ii<pR; ii++ )
         std::cout << "numDet= " <<aData->num_det() << " det_data= "<<aData->get(ii) << std::endl;
         } */
        
        
        theIcData.clear();
        
        //  for(int jj=0; jj<numDet[2]; jj++ ) {
        for(int jj=0; jj<gMaxSections; jj++ ) {
                if( nParam[2][jj] < 1 ) 
                        continue;
                aData = GetNextDetData(jj,pI); theIcData.push_back( aData );
                for(int ii=0; ii<pI; ii++ )
                        if (theMap[120 + 6*jj+ii]>0)
                                aData->set( ii, theMap[120 + 6*jj+ii] + Random::uniform() - 0.5 );
                        else 
                                aData->set( ii, 0 );    
        }
        
        
        theSideData.clear();
        
        //  for( int jj=0; jj<numDet[3]; jj++ ) {
        for(int jj=0; jj<gMaxSections; jj++ ) {
                if( nParam[3][jj] < 1 ) 
                        continue;
                aData = GetNextDetData(jj,pS); theSideData.push_back( aData );
                for(int ii=0; ii<pS; ii++ )
                        if (theMap[180 + 6*jj+ii]>0)
                                aData->set( ii, theMap[180 + 6*jj+ii] + Random::uniform() - 0.5 );
                        else 
                                aData->set( ii, 0 );    
        }
        
#if 0 
        
        cout << " Testing new event " << endl;
        for(int ii=0; ii<(int)theMcpData.size(); ii++ ) {
                aData = theMcpData[ii];
                cout << endl << "New MCP " << (*aData) << endl;
        }
        for(int ii=0; ii<(int)thePpacData.size(); ii++ ) {
                aData = thePpacData[ii];
                cout << endl << "New PPAC " << (*aData) << endl;
        }
        
        for(int ii=0; ii<(int)theIcData.size(); ii++ ) {
                aData = theIcData[ii];
                cout << endl << "New IC " << (*aData) << endl;
        }  
        
        for(int ii=0; ii<(int)theSideData.size(); ii++ ) {
                aData = theSideData[ii];
                cout << endl << "New Side " << (*aData) << endl;
        }
        
#endif  
        
        fPM->new_event( theMcpData, thePpacData, theIcData, theSideData );
}


ClassImp(RawPrisma);

RawPrisma::RawPrisma(const char *name, const char *title) : BasePrismaWatcher(name, title), 
        fFP( new TObjArray(gMaxSections) ),
        fLR( new TObjArray(gMaxSections) )
{       
        // histograms belongs to the pool
        fFP->SetOwner(false);
        fLR->SetOwner(false);
        
        TH2F *h;
        for (Int_t i = 0; i < gMaxSections; i++) {
                h = new TH2F(Form("TOF_Xfp%d",i),Form("TOF_Xfp%d",i), 512, 0, 4096, 512, 0, 4096);
                fFP->AddAt(h,i);
                AddToPool(h);
                h = new TH2F(Form("fLR%d",i),Form("fLR%d",i), 512, 0, 4096, 512, 0, 4096);
                fLR->AddAt(h,i);
                AddToPool(h);           
        }
} 

RawPrisma::~RawPrisma()
{
}

void RawPrisma::Exec(Option_t * option)
{       
        // No need to pass events to prismaManager, only the raw content is needed
        SetLastExecStatus(0);
        RancLegnaroWatcher::Exec(option);
        if ( GetLastExecStatus() != 0 ) 
                return;
        
        if ( lenrawBuf == 0 )
                return;
        
        //      fLR->At(0] ->Fill( Converter[2].GetValue( 3), Converter[2].GetValue( 0) + Converter[2].GetValue( 1));
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+3], rawBuf[2*gMaxChannels+0] + rawBuf[2*gMaxChannels+1] );   
        //      fLR->At(1] ->Fill( Converter[2].GetValue(23), Converter[2].GetValue( 2) + Converter[2].GetValue(22));
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+23], rawBuf[2*gMaxChannels+2] + rawBuf[2*gMaxChannels+22] );                                                 
        //      fLR->At(2] ->Fill( Converter[2].GetValue(24), Converter[2].GetValue( 4) + Converter[2].GetValue( 5));
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+24], rawBuf[2*gMaxChannels+4] + rawBuf[2*gMaxChannels+5] );                                                  
        //      fLR->At(3] ->Fill( Converter[2].GetValue(25), Converter[2].GetValue( 6) + Converter[2].GetValue( 7));
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+25], rawBuf[2*gMaxChannels+6] + rawBuf[2*gMaxChannels+7] );                                                  
        //      fLR->At(4] ->Fill( Converter[2].GetValue(26), Converter[2].GetValue( 8) + Converter[2].GetValue( 9));
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+26], rawBuf[2*gMaxChannels+8] + rawBuf[2*gMaxChannels+9] );                                                  
        //      fLR->At(5] ->Fill( Converter[2].GetValue(27), Converter[2].GetValue(10) + Converter[2].GetValue(11));
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+27], rawBuf[2*gMaxChannels+10] + rawBuf[2*gMaxChannels+11] );                                                        
        //      fLR->At(6] ->Fill( Converter[2].GetValue(28), Converter[2].GetValue(12) + Converter[2].GetValue(13));
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+28], rawBuf[2*gMaxChannels+12] + rawBuf[2*gMaxChannels+13] );                                                        
        //      fLR->At(7] ->Fill( Converter[2].GetValue(29), Converter[2].GetValue(14) + Converter[2].GetValue(15));  
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+29], rawBuf[2*gMaxChannels+14] + rawBuf[2*gMaxChannels+15] );                                                        
        //      fLR->At(8] ->Fill( Converter[2].GetValue(30), Converter[2].GetValue(16) + Converter[2].GetValue(17)); 
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+30], rawBuf[2*gMaxChannels+16] + rawBuf[2*gMaxChannels+17] );                                                        
        //      fLR->At(9] ->Fill( Converter[2].GetValue(31), Converter[2].GetValue(18) + Converter[2].GetValue(19)); 
        ((TH2F *)fLR->At(0)) ->Fill( rawBuf[2*gMaxChannels+31], rawBuf[2*gMaxChannels+18] + rawBuf[2*gMaxChannels+18] );                                                        
        
        ((TH2F *)fFP->At(0)) ->Fill( rawBuf[2*gMaxChannels+0]  - rawBuf[2*gMaxChannels+1],  rawBuf[22+0] );     
        ((TH2F *)fFP->At(1)) ->Fill( rawBuf[2*gMaxChannels+2]  - rawBuf[2*gMaxChannels+22], rawBuf[22+1] );                                                     
        ((TH2F *)fFP->At(2)) ->Fill( rawBuf[2*gMaxChannels+4]  - rawBuf[2*gMaxChannels+5],  rawBuf[22+2] );                                                     
        ((TH2F *)fFP->At(3)) ->Fill( rawBuf[2*gMaxChannels+6]  - rawBuf[2*gMaxChannels+7],  rawBuf[22+3] );                                                     
        ((TH2F *)fFP->At(4)) ->Fill( rawBuf[2*gMaxChannels+8]  - rawBuf[2*gMaxChannels+9],  rawBuf[22+4] );                                                     
        ((TH2F *)fFP->At(5)) ->Fill( rawBuf[2*gMaxChannels+10] - rawBuf[2*gMaxChannels+11], rawBuf[22+5] );                                                     
        ((TH2F *)fFP->At(6)) ->Fill( rawBuf[2*gMaxChannels+12] - rawBuf[2*gMaxChannels+13], rawBuf[22+6] );                                                     
        ((TH2F *)fFP->At(7)) ->Fill( rawBuf[2*gMaxChannels+14] - rawBuf[2*gMaxChannels+15], rawBuf[22+7] );                                                     
        ((TH2F *)fFP->At(8)) ->Fill( rawBuf[2*gMaxChannels+16] - rawBuf[2*gMaxChannels+17], rawBuf[22+8] );                                                     
        ((TH2F *)fFP->At(9)) ->Fill( rawBuf[2*gMaxChannels+18] - rawBuf[2*gMaxChannels+18], rawBuf[22+9] );                             
}

ClassImp(DisplayPrisma);

DisplayPrisma::DisplayPrisma(const char *name, const char *title): BasePrismaWatcher(name, title)
{
        /*      Ic_Range        = new TH2F("Ic_Range", "Ic_Range", 1024,0,2000,1024,0,2000); 
         AddToPool(Ic_Range); */
        IcABvsABCD       = new TH2F("Ic_ABvsABCD", "E_DeltaE_AB", 1024,0,12000,2048,0,8000); 
        AddToPool(IcABvsABCD);
        IcAvsABCD       = new TH2F("Ic_AvsABCD", "E_DeltaE_A", 1024,0,12000,2048,0,8000); 
        AddToPool(IcAvsABCD);
        TOF_D           = new TH2F("ToF(ns)vsX_FP(mm)", "ToFvsx_FP", 1024,0,1024,4096,0,4096);
        AddToPool(TOF_D);
        PRISMA_MCPcal   = new TH2F("MCPcal", "MCP calibrated", 1024,-50,50, 1024,-50,50); 
        AddToPool(PRISMA_MCPcal);
        PRISMA_X_FP     = new TH1I("PRISMA_Xfc","X focal plane in PRISMA",4096,0,4096);
        AddToPool(PRISMA_X_FP);
        PRISMA_Vel      = new TH1I("PRISMA_Vel","Velocity in PRISMA",8000,0,0.25);
        AddToPool(PRISMA_Vel);
        Range_Energy   = new TH2F("Range_Energy", "Range_Energy", 1024,0,4096, 1024,0,20000); 
        AddToPool(Range_Energy);
        Vel_Theta   = new TH2F("Vel_Theta", "Vel_Theta", 1023,1,1024, 1024,0,0.2); 
        AddToPool(Vel_Theta);
        a_over_q    = new TH2F("a_over_q", "a_over_q", 1024,0,1200, 4096,0,450); 
        AddToPool(a_over_q);
        Energy_RBeta = new TH2F("Energy_RBeta", "Energy_RBeta", 1024,0,4096, 1024,0,20000); 
        AddToPool(Energy_RBeta);
        PRISMA_Mass = new TH1F("PRISMA_Mas", "PRISMA_Mas", 4096,0,4096); 
        AddToPool(PRISMA_Mass);
} 

void DisplayPrisma::DoCanvas(TCanvas *c, Option_t * /*o*/)
{
        //      TString opt = o; 
        
        c->Divide(5,2,0.0001,0.0001);
        c->cd(1);
        IcABvsABCD->Draw("col");
        c->cd(2);
        IcAvsABCD->Draw("col");
        c->cd(3);
        TOF_D->Draw("col");
        c->cd(4);
        PRISMA_MCPcal->Draw("col");
        c->cd(5);
        PRISMA_X_FP->Draw();
        c->cd(6);
        Range_Energy->Draw("col");
        c->cd(7);
        PRISMA_Vel->Draw("");
        c->cd(8);
        a_over_q->Draw("col");
        c->cd(9);
        Energy_RBeta->Draw("col");
        c->cd(10);
        PRISMA_Mass->Draw();    
        /*
         if ( opt.Contains("cor") ) {
         cas = 1;
         }
         switch( cas ) {
         case 1:
         c->Divide(2,5,0.0001,0.0001);
         for (Int_t i = 0; i < CrystalInterface::kNbCores; i++ ) {
         c->cd(i+1);
         fCor->At(i)->Draw();   
         }
         break;
         default:
         // show segments
         c->Divide(6,6,0.0001,0.0001);
         for (Int_t i = 0; i < CrystalInterface::kNbSegments; i++ ) {
         c->cd(i+1);
         fSeg->At(i)->Draw();   
         }      
         }
         */
        
}


void DisplayPrisma::Exec(Option_t *option) 
{
        const double  c_light = 29.97295 * cm/ns;
        
        // fill the array of floats and compute the PRISMA physics
        SetLastExecStatus(0u);
        BasePrismaWatcher::Exec(option);
        if ( GetLastExecStatus() != 0 ) 
                return;
        
        if ( fPM->length() > 0. )  {
                
                Vector3D Vel = fPM->velocity()/c_light; 
                
                //              Ic_Range      -> Fill(fPM->range()/mm,fPM->ic_energy());
                IcABvsABCD    -> Fill(fPM->ic_energy(),fPM->ic_energy_AB_DE());
                //IcAvsABCD     -> Fill(fPM->ic_energy(),fPM->ic_energy_A_DE());
                TOF_D         -> Fill(fPM->x_fp() /mm, 10*fPM->tof()/ns);
                PRISMA_MCPcal -> Fill(fPM->mcp_x()/mm , fPM->mcp_y()/mm) ;
                PRISMA_X_FP   -> Fill( fPM->x_fp() /mm) ;
                PRISMA_Vel    -> Fill( Vel.rho() ) ;  
                Range_Energy->Fill(fPM->range(), fPM->ic_energy()) ;
                Vel_Theta->Fill(10.*fPM->theta_c()/degree, Vel.rho() );
                a_over_q->Fill(fPM->x_fp()/mm, fPM->a_over_q_uncal());
                Energy_RBeta->Fill(fPM->r_beta(),fPM->ic_energy());
                PRISMA_Mass->Fill(fPM->Mass());
        }
        
}

ClassImp(EvePrisma);
EvePrisma::EvePrisma(const char *name, const char *title)
          : BasePrismaWatcher(name, title), 
            fMaxCycle(100),
            fConfFileName(Gw::AgataEventDisplay::GetDefaultPrismaConf()),
                        fPathFile(Gw::AgataEventDisplay::GetDefaultAgataPath())
{
        fTypeName =  AgataEventDisplay::Instance()->RegisterBranch("Prisma");
        AgataEventDisplay::Instance()->SetTrackStyle("Rectangle", fTypeName);
} 

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

void EvePrisma::ShowEve(const Char_t* prismaConfFile, const Char_t* agataPathFile)
{
    Gw::AgataEventDisplay::SetDefaultPrismaConf(prismaConfFile);
    Gw::AgataEventDisplay::SetDefaultAgataPath(agataPathFile);
        Gw::AgataEventDisplay::Instance()->BuildDefaultGeometry(false, true);
        Gw::AgataEventDisplay::Instance()->ShowDisplay();
}

void EvePrisma::Exec(Option_t *option) 
{
        const double  c_light = 29.97295 * cm/ns;
        //static Int_t count = 0;
        
        // fill the array of floats and compute the PRISMA physics
        SetLastExecStatus(0);
        BasePrismaWatcher::Exec(option);
        if ( GetLastExecStatus() != 0 ) 
                return;
        
        Vector3D Vel         = fPM->velocity()/c_light; 
        Vector3D posQpoleIn  = fPM->getPosQuadrupoleIn();
        Vector3D posQpoleOut = fPM->getPosQuadrupoleOut();
        Vector3D posDipIn    = fPM->getPosDipoleIn();
        Vector3D posDipOut   = fPM->getPosDipoleOut();
        Vector3D posFP       = fPM->getPosFocalPlane();
        Vector3D center      = fPM->getCenter();
        Float_t radius      =  fPM->radius();
        
        if ( fPM->length() > 0. )  {
                
                AgataEventContainer* agataContainer = Gw::AgataEventDisplay::Instance()->GetEventContainer(); 
                TrackHit* hit = agataContainer->NewTrackHit(fTypeName);
                StdHit* thit = 0x0;
                
                // MCP
                thit = hit->NewHit();
                thit->SetX(-fPM->mcp_x()/cm);
                thit->SetY( fPM->mcp_y()/cm);
                thit->SetZ(25);
                thit->SetE(fPM->energy_over_beta()+100);
                
                // QPole
                thit = hit->NewHit();
                thit->SetX(-posQpoleIn.X()/cm);
                thit->SetY( posQpoleIn.Z()/cm);
                thit->SetZ( posQpoleIn.Y()/cm);
                thit->SetE(fPM->energy_over_beta()+100);
                
                thit = hit->NewHit();
                thit->SetX(-posQpoleOut.X()/cm);
                thit->SetY( posQpoleOut.Z()/cm);
                thit->SetZ( posQpoleOut.Y()/cm);
                thit->SetE(fPM->energy_over_beta()+100);
                
                // Dipole in
                thit = hit->NewHit();
                thit->SetX(-posDipIn.X()/cm);
                thit->SetY( posDipIn.Z()/cm);
                thit->SetZ( posDipIn.Y()/cm);
                thit->SetE(fPM->energy_over_beta()+100);
                
                // In Dipole
                Float_t phi1 = TMath::Abs(TMath::ATan((posDipIn.X()-center.X())/(posDipIn.Y()-center.Y())));
                Float_t phi2 = TMath::Abs(TMath::ATan((posDipOut.X()-center.X())/(posDipOut.Y()-center.Y())));
                Float_t dphi = TMath::Abs(phi2-phi1)/10.;
                
                for (Int_t i = 1; i < 10; ++i) {
                        thit = hit->NewHit();
                        thit->SetX(-(center.X()/cm - radius*TMath::Sin(phi1 - float(i)*dphi)/cm));
                        thit->SetY( posDipIn.Z()/cm);
                        thit->SetZ( center.Y()/cm  + radius*TMath::Cos(phi1 - float(i)*dphi)/cm);
                        thit->SetE(fPM->energy_over_beta()+100);
                }
                
                // Dipole Out
                thit = hit->NewHit();
                thit->SetX(-posDipOut.X()/cm);
                thit->SetY( posDipIn.Z()/cm);
                thit->SetZ( posDipOut.Y()/cm);
                thit->SetE(fPM->energy_over_beta()+100);
                
                // Focal Plane
                thit = hit->NewHit();
                thit->SetX(-posFP.X()/cm);
                thit->SetY( posFP.Z()/cm);
                thit->SetZ( posFP.Y()/cm);
                thit->SetE(fPM->energy_over_beta()+100);
                
                // Ionization Chamber
                if (fPM->path() > 0) {
                        Float_t phi = TMath::ATan(-(posFP.X()-posDipOut.X())/(posFP.Y()-posDipOut.Y()));
                        Float_t x   = -posFP.X()/cm + fPM->path()*TMath::Sin(phi)/10; // path in mm !
                        Float_t z   =  posFP.Y()/cm + fPM->path()*TMath::Cos(phi)/10;
                        thit = hit->NewHit();
                        thit->SetX(x);
                        thit->SetY(0.);
                        thit->SetZ(z);
                        thit->SetE(fPM->energy_over_beta()+100);
                }
                //        printf("%f %f %f %f\n", phi1*TMath::RadToDeg(), phi2*TMath::RadToDeg(),   posDipOut.X()/cm, center.X()/cm - radius*TMath::Sin(phi2)/cm);
                //        printf("%f %f %f %f\n\n", phi1*TMath::RadToDeg(), phi2*TMath::RadToDeg(), posDipOut.Y()/cm, center.Y()/cm + radius*TMath::Cos(phi2)/cm);
                
                agataContainer->FillTracks(fTypeName);
                
                if (agataContainer->GetNofEvents() > fMaxCycle) {
                        TTask::SetActive(false);
                }
        }
        
}

ClassImp(DoPrismaDoppler);

DoPrismaDoppler::DoPrismaDoppler(const char *name, const char *title): BasePrismaWatcher(name, title), VertexBuilder()
{
        fVertexFrame = MainFrameFactory::theMainFactory().
                NewSharedFrame( FactoryItem("Agata","meta:vertex",Version(0,0)), FactoryItem("Agata","meta:vertex",Version(1,1)) );
        
        fIsToBeDeleted = true;  
} 

void DoPrismaDoppler::Exec(Option_t *option) 
{       
        const double  c_light = 29.97295 * cm/ns;
        
        // fill the array of floats and compute the PRISMA physics
        SetLastExecStatus(0);
        BasePrismaWatcher::Exec(option);
        if ( GetLastExecStatus() != 0 ) 
                return;
        
        // now get Velocity vector from prismaManager and pass it to VertexWatcher
        Vector3D Vel; VertexInterface *gdata = GetVertex();
        
        if ( fPM->length() > 0. && fPM->tof()/ns > 100 && fPM->tof()/ns < 300 )  
                //      if ( fPM->length() > 0. )  
                Vel = fPM->velocity()/c_light;          
        else 
                Vel = Vector3D(0,0,0);
        
        gdata->SetPosition(0,0,0);
        gdata->SetDirection(Vel.X(),Vel.Y(),Vel.Z());   
        gdata->SetBeta(Vel.rho()); 
}

ClassImp(MyPrismaTree);

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

void MyPrismaTree::SetBranches()
{
        GetTree()->Branch("AoverQ", &fAoverQ, "AoverQ/F");
        GetTree()->Branch("Radius", &fRadius, "Radius/F");
        GetTree()->Branch("Ic_AB_DE", &fIc_AB_DE, "AoverQ/F");
        GetTree()->Branch("Ic", &fIc, "Ic/F");
        GetTree()->Branch("x_fp", &fx_fp, "x_fp/F");
        GetTree()->Branch("TOF", &fTOF, "TOF/F");
        GetTree()->Branch("Range", &fRange, "Range/F");
        GetTree()->Branch("Ic_A_DE",&fIc_A_DE,"Ic_A_DE");
}

MyPrismaTree::~MyPrismaTree()
{
}

void MyPrismaTree::Exec(Option_t *option) 
{
        if ( gDebug > 3 ) 
                printf("MyPrismaTree::Exec is called \n");
        
        // fill the array of floats and compute the PRISMA physics
        SetLastExecStatus(0);
        BasePrismaWatcher::Exec(option);
        if ( GetLastExecStatus() != 0 ) 
                return;
        
        fAoverQ = fPM->a_over_q_uncal();
        fRadius = fPM->path(); 
        fIc_AB_DE = fPM->ic_energy_AB_DE();
        fIc = fPM->ic_energy();
        fx_fp = fPM->x_fp() /mm; 
        fTOF = fPM->tof()/ns; 
        fRange = fPM->range();
        fIc_A_DE = fPM->ic_energy_A_DE();
        
        // fill the tree only if this is the owner of the tree ... otherwise by TreeMaster
        FillTree();
}