Analysis.cxx 12.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
/*****************************************************************************
 * Copyright (C) 2009-2016    this file is part of the NPTool Project        *
 *                                                                           *
 * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
 * For the list of contributors see $NPTOOL/Licence/Contributors             *
 *****************************************************************************/

/*****************************************************************************
 * Original Author: Morfouace Pierre  contact address: morfouace@ganil.fr    *
 *                                                                           *
 * Creation Date  : April 2018                                               *
 * Last update    :                                                          *
 *---------------------------------------------------------------------------*
 * Decription:                                                               *
 *  This class describe  Actar analysis project                              *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * Comment:                                                                  *
 *                                                                           *
 *****************************************************************************/

#include<iostream>

using namespace std;
#include"Analysis.h"
#include"NPAnalysisFactory.h"
#include"NPDetectorManager.h"
#include"RootOutput.h"
#include"RootInput.h"


////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
////////////////////////////////////////////////////////////////////////////////
Analysis::~Analysis(){
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
    Actar= (TActarPhysics*) m_DetectorManager->GetDetector("Actar");
    ReactionConditions = new TReactionConditions();
43
    
44 45 46 47
    Actar->ReadAnalysisConfig();
    if(Actar->GetRansacStatus()){
        Actar->SetRansacParameter("./configs/RansacConfig.dat");
    }
48 49 50 51
    else if(Actar->GetClusterStatus()){
        Actar->SetClusterParameter("./configs/ClusterConfig.dat");
    }
    
52 53 54 55 56
    DriftVelocity = Actar->GetDriftVelocity();
    PadSizeX = Actar->GetPadSizeX();
    PadSizeY = Actar->GetPadSizeY();
    NumberOfPadsX = Actar->GetNumberOfPadsX();
    NumberOfPadsY = Actar->GetNumberOfPadsY();
57
    
58 59 60
    EnergyLoss_1H = NPL::EnergyLoss("./EnergyLossTable/proton_iC4H10_6.24151e+07_295.G4table","G4Table",100);
    EnergyLoss_18O = NPL::EnergyLoss("./EnergyLossTable/O18_iC4H10_6.24151e+07_295.G4table","G4Table",100);
    TheReaction = new NPL::Reaction("18O(p,p)18O@59");
61
    
62 63 64 65 66 67 68
    InitInputBranch();
    InitOutputBranch();
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
    ReInitValue();
69
    
70 71
    LightName="";
    if(ReactionConditions->GetParticleMultiplicity()>0){
72 73 74 75
        InitXVertex = ReactionConditions->GetVertexPositionZ();
        InitTheta3 = ReactionConditions->GetTheta(0);
        InitE3 = ReactionConditions->GetKineticEnergy(0);
        LightName = ReactionConditions->GetParticleName(0);
76 77
    }
    int TrackMult = Actar->GetTrackMult();
78
    
79 80
    TVector3 vX = TVector3(1,0,0);
    TVector3 aTrack, vB;
Morfouace's avatar
Morfouace committed
81
    if(TrackMult>1){
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
        vTrack = Actar->GetTracks();
        double scalarproduct=0;
        int BeamTrack=0;
        for(unsigned int i=0; i<TrackMult; i++){
            TVector3 vtest = TVector3(vTrack[i].GetDirectionVector().X(),vTrack[i].GetDirectionVector().Y(),vTrack[i].GetDirectionVector().Z());
            TVector3 vunit = vtest.Unit();
            double scalar = abs(vunit.Dot(vX));
            vScalar.push_back(scalar);
            //cout << scalar << endl;
            //cout << scalarproduct << endl;
            if(scalar>scalarproduct){
                BeamTrack=i;
                scalarproduct=scalar;
            }
        }
97
        
98 99 100 101
        double XBeam = vTrack[BeamTrack].GetDirectionVector().X();
        double YBeam = vTrack[BeamTrack].GetDirectionVector().Y();
        double ZBeam = vTrack[BeamTrack].GetDirectionVector().Z();
        TVector3 vBeam = TVector3(XBeam,YBeam,ZBeam);
102
        
103 104 105 106
        double XBeamPoint = vTrack[BeamTrack].GetXh();
        double YBeamPoint = vTrack[BeamTrack].GetYh();
        double ZBeamPoint = vTrack[BeamTrack].GetZh();
        TVector3 vBeamPos = TVector3(XBeamPoint,YBeamPoint,ZBeamPoint);
107
        
108 109
        //vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam*DriftVelocity);
        vB = TVector3(XBeam*PadSizeX, YBeam*PadSizeY,ZBeam);
110
        BeamAngle = (vX.Angle(vB))*180/TMath::Pi();
111
        
112 113 114 115 116
        for(unsigned int i=0; i<TrackMult; i++){
            if(i!=BeamTrack){
                double Xdir = vTrack[i].GetDirectionVector().X();
                double Ydir = vTrack[i].GetDirectionVector().Y();
                double Zdir = vTrack[i].GetDirectionVector().Z();
117
                
118 119
                double vertex_x = vTrack[i].GetVertexPostion(vBeam,vBeamPos).X()*PadSizeX;
                double vertex_y = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Y()*PadSizeY;
120 121
                //double vertex_z = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z()*DriftVelocity;
                double vertex_z = vTrack[i].GetVertexPostion(vBeam,vBeamPos).Z();
122
                
123 124
                //aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir*DriftVelocity);
                aTrack = TVector3(Xdir*PadSizeX, Ydir*PadSizeY, Zdir);
125 126 127
                double angle = vX.Angle(aTrack)*180/TMath::Pi();
                //double angle = vB.Angle(aTrack)*180/TMath::Pi();
                if(angle>90) angle = 180-angle;
128
                
129 130 131 132 133 134
                double x1 = vTrack[i].GetXm()*PadSizeX;
                double x2 = vTrack[i].GetXh()*PadSizeX;
                double y1 = vTrack[i].GetYm()*PadSizeY-0.5*NumberOfPadsY*PadSizeY;
                double y2 = vTrack[i].GetYh()*PadSizeY-0.5*NumberOfPadsY*PadSizeY;
                //double z1 = -(vTrack[i].GetZm()-256)*DriftVelocity;
                //double z2 = -(vTrack[i].GetZh()-256)*DriftVelocity;
135 136 137 138
                //double z1 = vTrack[i].GetZm()*DriftVelocity;
                double z1 = vTrack[i].GetZm();
                //double z2 = vTrack[i].GetZh()*DriftVelocity;
                double z2 = vTrack[i].GetZh();
139 140
                
                
141 142 143
                if(vertex_x>0 && vertex_x<256){
                    double LengthInGas = fSiDistanceX - vertex_x;
                    for(unsigned int k=0; k<Actar->Si_E.size(); k++){
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
                        XVertex.push_back(vertex_x);
                        YVertex.push_back(vertex_y);
                        ZVertex.push_back(vertex_z);
                        
                        GetMayaSiHitPosition(x1,x2,y1,y2,z1,z2);
                        ESi.push_back(Actar->Si_E[k]);
                        SiNumber.push_back(Actar->Si_Number[k]);
                        
                        DE.push_back(vTrack[i].GetPartialCharge(108,128)/(20./cos(angle*TMath::Pi()/180)));
                        double E3;
                        
                        if(LightName=="proton")E3 = EnergyLoss_1H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
                        if(LightName=="deuteron")E3 = EnergyLoss_2H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
                        if(LightName=="triton")E3 = EnergyLoss_3H.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
                        if(LightName=="He3")E3 = EnergyLoss_3He.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
                        if(LightName=="alpha")E3 = EnergyLoss_4He.EvaluateInitialEnergy(Actar->Si_E[k]*MeV,LengthInGas*mm,angle*TMath::Pi()/180);
                        //double BeamEnergy = EnergyLoss_18O.Slow(59.4*MeV,(vertex_x[i]+60)*mm, BeamAngle*TMath::Pi()/180);
                        double Energy_beam = EnergyLoss_18O.Slow(59.4*MeV,(vertex_x+60)*mm, 0);
                        BeamEnergy.push_back(Energy_beam);
                        TheReaction->SetBeamEnergy(Energy_beam);
                        ELab.push_back(E3);
                        ThetaLab.push_back(angle);
                        TheReaction->SetNuclei3(E3,angle*TMath::Pi()/180);
                        Ex.push_back(TheReaction->GetExcitation4());
                        ThetaCM.push_back(TheReaction->GetThetaCM()*180./TMath::Pi());
169 170 171 172 173 174 175 176 177 178 179
                    }
                }
            }
        }
    }
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::GetMayaSiHitPosition(double xm, double xh, double ym, double yh, double zm, double zh)
{
    double X1, X2, Y1, Y2, Z1, Z2;
180
    
181
    if(xm>xh){
182 183 184 185 186 187 188
        X1 = xh;
        Y1 = yh;
        Z1 = zh;
        
        X2 = xm;
        Y2 = ym;
        Z2 = zm;
189 190
    }
    else if(xh>xm){
191 192 193 194 195 196 197
        X1 = xm;
        Y1 = ym;
        Z1 = zm;
        
        X2 = xh;
        Y2 = yh;
        Z2 = zh;
198
    }
199
    
200 201
    double l, L, t;
    double zf, yf;
202
    
203
    if(fSiDistanceX>X2){
204 205 206 207 208 209
        L = fSiDistanceX-X2;
        l = X2 - X1;
        t = (l+L)/l;
        
        zf = Z1 + (Z2-Z1)*t;
        yf = Y1 + (Y2-Y1)*t;
210 211
    }
    else if(fSiDistanceX<X2){
212 213 214 215 216 217
        L = X2 - fSiDistanceX;
        l = fSiDistanceX - X1;
        t = (l+L)/l;
        
        zf = Z1 + (Z2-Z1)/t;
        yf = Y1 + (Y2-Y1)/t;
218
    }
219
    
220 221 222 223 224 225 226 227 228 229
    SiPosY.push_back(yf);
    SiPosZ.push_back(zf);
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch() {
230 231 232
    RootInput::getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true);
    RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*",true);
    RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
    RootOutput::getInstance()->GetTree()->Branch("DE",&DE);
    RootOutput::getInstance()->GetTree()->Branch("SiNumber",&SiNumber);
    RootOutput::getInstance()->GetTree()->Branch("ESi",&ESi);
    RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab);
    RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
    RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex);
    RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM);
    RootOutput::getInstance()->GetTree()->Branch("vScalar",&vScalar);
    RootOutput::getInstance()->GetTree()->Branch("XVertex",&XVertex);
    RootOutput::getInstance()->GetTree()->Branch("YVertex",&YVertex);
    RootOutput::getInstance()->GetTree()->Branch("ZVertex",&ZVertex);
    RootOutput::getInstance()->GetTree()->Branch("BeamAngle",&BeamAngle,"BeamAngle/D");
    RootOutput::getInstance()->GetTree()->Branch("BeamEnergy",&BeamEnergy);
    RootOutput::getInstance()->GetTree()->Branch("SiPosY",&SiPosY);
    RootOutput::getInstance()->GetTree()->Branch("SiPosZ",&SiPosZ);
    RootOutput::getInstance()->GetTree()->Branch("InitXVertex",&InitXVertex,"InitXVertex/D");
    RootOutput::getInstance()->GetTree()->Branch("InitE3",&InitE3,"InitE3/D");
    RootOutput::getInstance()->GetTree()->Branch("InitTheta3",&InitTheta3,"InitTheta3/D");
254
    
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
    DE.clear();
    SiNumber.clear();
    ESi.clear();
    ELab.clear();
    ThetaLab.clear();
    Ex.clear();
    ThetaCM.clear();
    vScalar.clear();
    XVertex.clear();
    YVertex.clear();
    ZVertex.clear();
    SiPosY.clear();
    SiPosZ.clear();
    BeamEnergy.clear();
273
    
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
    BeamAngle=-1000;
    InitE3=-1000;
    InitTheta3=-1000;
    InitXVertex=-1000;
}

////////////////////////////////////////////////////////////////////////////////
//            Construct Method to be pass to the DetectorFactory              //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
    return (NPL::VAnalysis*) new Analysis();
}

////////////////////////////////////////////////////////////////////////////////
//            Registering the construct method to the factory                 //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
    class proxy{
    public:
        proxy(){
            NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
        }
    };
297
    
298 299
    proxy p;
}