ROOT logo
#ifndef _HBEAM_H_
#define _HBEAM_H_

//#define isPluto


/////////////////////////////////////////////////////////////////////
//
// Beam line transport
//
// Author:  T. Hennino (core), J. Biernat , J.Markert
//          Pluto implementation by IF
// Written: 14.11.2013
//
/////////////////////////////////////////////////////////////////////



//_HADES_CLASS_DESCRIPTION
///////////////////////////////////////////////////////////////////////////////
//
// HBeam
//
//    This class simulates a beam line. Beam particles (HBeamParticle) will be produced using
//    a certain beam profile and momentum smearing. The particle will be tranported through
//    the beam line using matrices for the single beam elements (HBeamElement). The particles
//    cand be registered in detectors along the beam line.
//
//    HBeam pionbeam;
//    pionbeam.setBeam           (HPhysicsConstants::pid("pi-"),1.3,60,60,0.0,0.0); // id, total mom [GeV], beamtube x and y size, xoff,yoff
//    pionbeam.SetBeamResolution (0.01,0.05,0.06);                                  // dpx [rad],dpy [rad],dpz [frac]  [+-] white random
//    pionbeam.setBeamProfile    (.05,0.0);                                         // sigma [mm], flatradius [mm]
//                                                                                  // sigma>0 will select the beam profile (gausian + flat top)
//                                                                                  // sigma==0, flatradius>0 will give an extended beam spot without gaussian borders
//                                                                                  // sigma==0,flatradius==0 spot like beam
//    if(!pionbeam.initBeamLine  ("pibeam_set6_mod.data",32)) return;               // transform input file and target element number
//    pionbeam.addDetector("det1", -17092.6,2,50.,50.);                             // [mm] relative to HADES 0,0,0    cutype (0,1,2), xcut[mm],ycut[mm]
//    pionbeam.addDetector("det2",  -5400.0,2,50.,50.);                             // [mm] relative to HADES 0,0,0    cutype (0,1,2), xcut[mm],ycut[mm]
//    pionbeam.addDetector("plane", -1300.0,1,60.,60.);                             // [mm] relative to HADES 0,0,0    cutype (0,1,2), xcut[mm],ycut[mm]
//
//    pionbeam.printBeamLine(kTRUE);          // kTRUE : print transform matrices in addition to name and distance
//    pionbeam.printDetectors();
//    pionbeam.printBeamProfile();
//
//
//    //-----------------------------------------------------------
//    // create particles
//    vector<HBeamParticle>& vhistory = pionbeam.newParticle();  // beam particle at : beam,det1,det2,plane.target
//    // check if particle was in acceptance (using the particle at the end of transport)
//    if(vhistory[vhistory.size()-1].fAccepted) {   ......   }
//
//    // get access to all beam elements and detectors (positions,accptance,statistics ...)
//    vector<HBeamElement>& elements  = pionbeam.getElements();
//    vector<HBeamElement>& detectors = pionbeam.getDetectors();
//
///////////////////////////////////////////////////////////////////////////////

#include "TString.h"
#include "TVector3.h"
#include "TRandom.h"
#include "TF1.h"
#include "TMath.h"
#include "TArrayD.h"
#include "TObject.h"

#include <sstream>
#include <iostream>
#include <iomanip>
#include <fstream>

using namespace std;


class HBeamElement : public TObject {

public:

    TString  fName;          // name of the lement : default  = Element[i]
    Double_t fDistance;      // position along beam line
    Int_t    fId;            // element number
    Double_t Tij [5][5];     // first order transform
    Double_t Tijk[5][5][5];  // second order transform
    Double_t fout[5];        // particle state at element [x,tx,y,ty,dp]
    Double_t fxCut;          // cut ond x for apperture
    Double_t fyCut;          // cut ond y for apperture
    Int_t    fCutType;       // 0: no cut ,1 : radial cut (default), 2: square cut
    // square cut uses x,y , radial cut only x!
    Long64_t fCtAll ;        // count all beam particle transport
    Long64_t fCtFail;        // count if particle is out of apperture
    Bool_t   fAccepted;      // remembers the last result of check()

    HBeamElement(TString n="not_set",Double_t dist = 0, Int_t id=-1)
    {
	fName     = n;
	fDistance = dist;
	fId       = id;
	fCutType  = 1;
	fxCut     = 60;
	fyCut     = 60;
	fCtAll    = 0;
	fCtFail   = 0;
	fAccepted = kTRUE;
	clear();
    }

    void setElement(TString name,Int_t cuttype,Double_t xcut,Double_t ycut){
	if(name.CompareTo("")!=0) fName = name;
	fCutType = cuttype;
	fxCut = xcut;
        fyCut = ycut;
    }

    void clear()
    {
	for(Int_t i = 0; i < 5; i ++) {
	    fout[i] = 0;
	    for(Int_t j = 0; j < 5; j ++){
		Tij[i][j] =0;
		for(Int_t k = 0; k < 5; k ++){
		    Tijk[i][j][k] = 0;
		}
	    }
	}
    }

    void printTij()
    {
	for(Int_t i = 0; i < 5; i ++) {     // row
	    for(Int_t j = 0; j < 5; j ++){  // column
		cout<<setw(12)<<Tij[i][j]<<" "<<flush;

	    }
	    cout<<endl;
	}
    }

    void printTijk()
    {
	for(Int_t i = 0; i < 5; i ++) { // block
	    for(Int_t k = 0; k < 5; k ++){ // line
		for(Int_t j = 0; j <= k; j ++){    // values
		    cout<<setw(12)<<Tijk[i][j][k]<<" "<<flush;
		}
		cout<<endl;
	    }
	    cout<<endl;
	}
    }


    void print(Bool_t printAll = kFALSE)
    {
	// printAll = kTRUE will print transforms in addition
	cout<<"---------------------------------"<<endl;
	cout<<"ID  : "<<setw(3)<<fId<<" , "<<fName.Data()<<endl;
	cout<<"distance   : "<<setw(12)<<fDistance<<endl;
	cout<<"CutType    : "<<setw(3)<<fCutType<<" , x : "<<setw(12)<<fxCut<<", y : "<<setw(12)<<fyCut<<endl;
	cout<<"Acceptance : "<<setw(8)<<fCtAll<<" , Failed : "<<setw(8)<<fCtFail<<" ,  "<<setw(8)<<(fCtAll > 0?  (fCtFail/(Double_t)fCtAll)*100: 100)<<" %"<<endl;
        if(printAll){
	    cout<<"first  order transform:"<<endl;
	    printTij();
	    cout<<"second order transform:"<<endl;
	    printTijk();
	}
    }

    Bool_t isInAcceptance()
    {
	// checks if the particle is inside acceptance
	// for this element. to be called after beam transport (for elements)
	// and after createDetectorHits() for detectors
	if(fCutType == 0 || fCutType > 2) return kTRUE;
	if(fCutType == 1 ){
	    Double_t radius = TMath::Sqrt(fout[0] * fout[0] + fout[2] * fout[2]) ;
	    if(radius <  fxCut )   return kTRUE;
	    else                   return kFALSE;
	}
	if(fCutType == 2 ){
	    if(TMath::Abs(fout[0]) < fxCut && TMath::Abs(fout[2]) < fyCut ) return kTRUE;
	    else                                                            return kFALSE;
	}
	return kTRUE;
    }

    void check()
    {
	// checks the acceptance and fill the stats counters
	// and the fAccepted flag
	if(!isInAcceptance()) {
	    fCtFail++;
	    fAccepted = kFALSE;
	} else  fAccepted = kTRUE;

	fCtAll++;
    }

    void  toLinearArray(TArrayD& all)
    {
	// store in linearized fashion  i*5+j  and  i*25+j*5+k fID,fDistance,fxCut,fyCut,fCuttype
        // last element is the element id

	all.Set(5*5 + 5*5*5 + 5);  // 1st,2nd, fID,fDistance,fxCut,fyCut,fCuttype

	for(Int_t i = 0; i < 5; i ++) {
	    for(Int_t j = 0; j < 5; j ++){
		all[i*5+j]=Tij[i][j];
	    }
	}
	for(Int_t i = 0; i < 5; i ++) {
	    for(Int_t j = 0; j < 5; j ++){
		for(Int_t k = 0; k < 5; k ++){
		    all[25+i*25+j*5+k]=Tijk[i][j][k];
		}
	    }
	}
        all[150+0] = fId;
        all[150+1] = fDistance;
        all[150+2] = fxCut;
        all[150+3] = fyCut;
        all[150+4] = fCutType;
    }

    void  fromLinearArray(TArrayD& all)
    {
	// store in linearized fashion  i*5+j  and  i*25+j*5+k   fID,fDistance,fxCut,fyCut,fCuttype
        // last element is the element id

	for(Int_t i = 0; i < 5; i ++) {
	    for(Int_t j = 0; j < 5; j ++){
                Tij[i][j]=all[i*5+j];
	    }
	}
	for(Int_t i = 0; i < 5; i ++) {
	    for(Int_t j = 0; j < 5; j ++){
		for(Int_t k = 0; k < 5; k ++){
		    Tijk[i][j][k]=all[25+i*25+j*5+k];
		}
	    }
	}
        fId       = all[150+0];
        fDistance = all[150+1];
        fxCut     = all[150+2];
        fyCut     = all[150+3];
        fCutType  = all[150+4];
    }
#ifndef isPluto
     ClassDef(HBeamElement,1)
#endif
};


class HBeamParticle : public TObject {

public:
    TVector3 fP;                 // momentum vector
    TVector3 fPos;               // position vector
    TVector3 fOffset;            // beam could be displaced : profile function makes use of it
    TVector3 fMomRes;            // momentum resolution

    Double_t fBeamMom;           // nominal beam momentum
    Double_t fBeamMomSmeared;    // smeared beam momentum
    Int_t    fPid;               // type of beam particle
    Int_t    fDetnum;            // detector number where the transported particle was registered
    TString  fName;              // name of detecor where particle was registered
    Double_t fbeamtube_size_x;   // x general acceptance of the beam tube
    Double_t fbeamtube_size_y;   // y general acceptance of the beam tube
    Bool_t   fAccepted;          // whether the particle was reaching the detector without violation acceptance

    HBeamParticle(){
	clear();
    }

    void clear (){
	fP     .SetXYZ(0.,0.,0.);
	fPos   .SetXYZ(0.,0.,0.);
	fOffset.SetXYZ(0.,0.,0.);
	fMomRes.SetXYZ(0.,0.,0.);
	fBeamMom         = -1.;
	fBeamMomSmeared  = -1.;
	fPid             = -10;
	fDetnum          = -1;
	fbeamtube_size_x = -1.;
	fbeamtube_size_y = -1.;
	fAccepted        = kTRUE;
    }


    void print(Bool_t printAll = kFALSE){
	// printAll=kTRUE print in addition momentum vector
	cout<<"---------------------------------"<<endl;
	cout<<setw(12)<<fName.Data()<<" acc : "<<fAccepted<<" , ID "<<fPid<<" , Beam Mom : "<<setw(12)<<fBeamMom<<endl;
	cout<<"Beam tube [x,y]       : "<<setw(12)<<fbeamtube_size_x <<" , "<<setw(12)<<fbeamtube_size_y <<endl;
	cout<<"P smear   [x,y,z]     : "<<setw(12)<<fMomRes.X()<<" , "<<setw(12)<<fMomRes .Y()<<" , "<<setw(12)<<fMomRes.Z()<<endl;
	cout<<"Position  [x,y,z]     : "<<setw(12)<<fPos   .X()<<" , "<<setw(12)<<fPos    .Y()<<" , "<<setw(12)<<fPos   .Z()<<endl;
	cout<<"Offset    [x,y,z]     : "<<setw(12)<<fOffset.X()<<" , "<<setw(12)<<fOffset .Y()<<" , "<<setw(12)<<fOffset.Z()<<endl;

	if(printAll){
	    cout<<"P         [x,y,z,tot] : "<<setw(12)<<fP     .X()<<" , "<<setw(12)<<fP      .Y()<<" , "<<setw(12)<<fP     .Z()<<" , "<<setw(12)<<fP  .Mag()<<endl;
	    cout<<"Detector ID           : "<<setw(12)<<fDetnum    <<endl;
	}


    }
#ifndef isPluto
     ClassDef(HBeamParticle,1)
#endif

};




class HBeam : public TObject {



private:

    static Double_t profile(Double_t *x, Double_t *par)
    {
	Double_t xx    = x[0];
	Double_t c     = 1;
	Double_t mean  = 0;
	Double_t sigma = par[0];
	Double_t flat  = par[1];
	Double_t result = 0;

	if(xx < mean - flat) {         // left gaus
	    result = c * TMath::Gaus(xx,mean - flat,sigma);
	} else if(xx > mean + flat) {  // right gaus
	    result = c * TMath::Gaus(xx,mean + flat,sigma);
	} else {                       // flat top
	    result = c;
	}

	return result;
    }



    HBeamParticle fBeam;                 // beam particle
    TF1*          fProfile;              //-> beam profile

    vector <HBeamElement> fdetectors;    // vector of all detectors
    vector <HBeamElement>  felements;    // vector of all beam elements
    vector <HBeamParticle>  fhistory;    //! vector of beam particle at all detectors
    Int_t                fnum_target;    // index of elemnt target
    Double_t            ftoPluto  [5];   // units -> transport ->Pluto
    Double_t            ffromPluto[5];   // units -> Pluto     ->transport
    Bool_t setTargetElement(UInt_t n)
    {
	if (n < 1 || n >= felements.size()){
	    cout<<"HBeam::setTargetElement() : specifified target number outside range!"<<endl;
	    return kFALSE;
	}
	fnum_target = n;
	Double_t shift = felements[fnum_target].fDistance;
	for (UInt_t i = 0; i < felements.size(); i ++){
	    felements[i].fDistance -= shift;
	}
	return kTRUE;
    }

public:

    HBeam()
    {
	fnum_target = -1;
	fProfile = new TF1("fprofile",profile,-20,20,2);
	fProfile->SetParNames("sigma","flatradius");
	setBeamProfile(0.05,0);
	setBeamResolution(0.01,0.05,0.06); //
	ftoPluto[0] = 10.;      // x  cm   -> mm
	ftoPluto[1] =  0.001;   // tx mrad -> rad
	ftoPluto[2] = 10.;      // y  cm   -> mm
	ftoPluto[3] =  0.001;   // ty mrad -> rad
	ftoPluto[4] =  0.01;    // dp %    -> frac

	for(Int_t i = 0; i < 5; i ++) ffromPluto[i] = 1./ftoPluto[i];
    }

    ~HBeam()
    {
	if(fProfile) delete fProfile;
    }
    vector<HBeamElement>&    getDetectors()     { return fdetectors ;  }
    vector<HBeamElement>&    getElements()      { return felements ;   }
    vector<HBeamParticle>&   getHistory()       { return fhistory ;    }

    void                     setTargetElementOnly(Int_t id)        { fnum_target=id;      }

    Int_t addDetector(TString name, Double_t distance,Int_t cutType,Double_t xcut=50,Double_t ycut=50)
    {
	// adds a detector "name" at place "distance" [in mm]
	// (distance must be negative if placed before the target)
	// cutType (0,1,2)(no,radial,square) xcut/ycut [mm]

	HBeamElement det(name,distance,fdetectors.size()+1);
	det.fxCut    = TMath::Abs(xcut);
	det.fyCut    = TMath::Abs(ycut);
	det.fCutType = cutType;

	fdetectors.push_back(det);
	return fdetectors.size();
    }

    void setBeam(Int_t id,
		 Double_t p,
		 Double_t dbeamtube_size_x,Double_t dbeamtube_size_y,
		 Double_t xoff = 0 ,Double_t yoff = 0)
    {
	fBeam.fBeamMom = TMath::Abs(p);
	fBeam.fPid     = id;
	fBeam.fOffset.SetXYZ(xoff,yoff,0.);
	fBeam.fbeamtube_size_x = TMath::Abs(dbeamtube_size_x);
	fBeam.fbeamtube_size_y = TMath::Abs(dbeamtube_size_y);
    }

    void setBeamResolution(Double_t dpx = 0.001,Double_t dpy = 0.005,Double_t dpz = 0.06)
    {
	fBeam.fMomRes.SetXYZ(TMath::Abs(dpx),TMath::Abs(dpy),TMath::Abs(dpz));
    }

    void setBeamProfile(Double_t sigma_beam = 0.05,Double_t flat_radius = 0)
    {
	fProfile->SetParameter(0,TMath::Abs(sigma_beam));
	fProfile->SetParameter(1,TMath::Abs(flat_radius));
    }


    void createBeam(HBeamParticle& part)
    {
	//----------------------------------------
	// beam profile
	// x-y symmetric profile
	Double_t x = 0;
	if      ( fProfile->GetParameter(0) != 0)                                  x  = fProfile->GetRandom();
	else if ( fProfile->GetParameter(0) == 0  && fProfile->GetParameter(1)!=0) x  = 2 * (gRandom->Rndm() - 0.5) * fProfile->GetParameter(3);
	else                                                                       x  = 0;
	Double_t phi = gRandom->Rndm()*360.;
	TVector3 p1(x,0.,felements[fnum_target].fDistance);
	p1.RotateZ(phi*TMath::DegToRad());
	p1 += fBeam.fOffset;  // beam displacement
	part.fPos += p1;
	// beam divergence
	Double_t p  = fBeam.fBeamMom + fBeam.fBeamMom * (gRandom->Rndm() - 0.5) * 2*fBeam.fMomRes.Z();
	Double_t px = p * (gRandom->Rndm() - 0.5) * 2*fBeam.fMomRes.X();
	Double_t py = p * (gRandom->Rndm() - 0.5) * 2*fBeam.fMomRes.Y();
	Double_t pz = TMath::Sqrt(p*p - px*px - py*py);

	part.fP.SetXYZ(px,py,pz);
	part.fBeamMomSmeared = p;
	part.fName   = "beam";
	part.fDetnum = -1;
	//----------------------------------------
    }

    Bool_t transport(HBeamParticle& part)
    {

	//----------------------------------------
	// coordinates  x [cm], px/pz [mrad],y [mm], py/pz [mrad], dp to cental p [%]
	//
	//
	//   x_i =  sum (j=0,5) T_ij*x_j + sum (j=0,5)(k=j,5)T_ijk*x_j*x_k
	//
	//
	//

	Double_t state   [5] = { 0,  0, 0,  0, 0};
	Double_t stateold[5] = {
	    part.fPos.X()                        *ffromPluto[0],
	    TMath::ATan2(part.fP.X(),part.fP.Z())*ffromPluto[1],
	    part.fPos.Y()                        *ffromPluto[2],
	    TMath::ATan2(part.fP.Y(),part.fP.Z())*ffromPluto[3],
	    ((part.fBeamMomSmeared-fBeam.fBeamMom)/fBeam.fBeamMom)*ffromPluto[4]
	};

	for(UInt_t el = 0; el < felements.size(); el++){
	    for(Int_t i = 0; i < 5; i ++){ // state vars
		// first order Term
		for(Int_t j = 0; j < 5; j ++){ //
		    state[i] += felements[el].Tij[i][j] * stateold[j];
		}
		// second order Term
		for(Int_t j = 0; j < 5; j ++){ //
		    for(Int_t k = j; k < 5; k ++){ //
			state[i] += felements[el].Tijk[i][j][k] * stateold[j] * stateold[k];
		    }
		}
	    }

	    // propagate and reset
	    for(Int_t i = 0; i < 5; i ++) {
		felements[el].fout[i] = state[i] * ftoPluto[i];
		state[i] = 0.;
	    }
	    felements[el].fout[4] = ((part.fBeamMomSmeared-fBeam.fBeamMom)/fBeam.fBeamMom) ; // keep dp
	    felements[el].check(); // stats and accepted flag
	}

	return kTRUE;
    }

    Bool_t createDetectorHits(HBeamParticle& part)
    {
	//calc particles at detector places
	Double_t state   [5] = { 0,  0, 0,  0, 0};
	Double_t p,px,py,pz;
	Bool_t   accepted    = kTRUE;
	Bool_t   acceptedDet = kTRUE;

	for (UInt_t i = 0; i < fdetectors.size(); i ++) {
	    accepted = kTRUE;
	    for (UInt_t k = 0; k < felements.size() - 1; k ++) {
		if(!felements[k].fAccepted) {
		    accepted = kFALSE;
		}
		if (felements[k].fDistance <= fdetectors[i].fDistance &&
		    fdetectors[i].fDistance < felements[k+1].fDistance) {

		    Double_t frac = (fdetectors[i].fDistance - felements[k].fDistance) / (felements[k+1].fDistance - felements[k].fDistance);

		    for(Int_t j = 0; j < 5; j ++){
			state[j] = felements[k].fout[j] + frac*(felements[k+1].fout[j] - felements[k].fout[j]);
		    }

		    p  = fBeam.fBeamMom * (1. + state[4]);
		    pz = p / sqrt(1.+ tan(state[1])*tan(state[1]) + tan(state[3])*tan(state[3]));
		    px = pz * tan(state[1]);
		    py = pz * tan(state[3]);

		    part.fP  .SetXYZ(px,py,pz);
		    part.fPos.SetXYZ(state[0],state[2],fdetectors[i].fDistance);
		    part.fDetnum   = fdetectors[i].fId;
		    part.fName     = fdetectors[i].fName;
		    part.fAccepted = accepted;

		    fdetectors[i].fout[0] = state[0] ;
		    fdetectors[i].fout[2] = state[2] ;
		    fdetectors[i].check();

		    if(!fdetectors[i].fAccepted) {
			part.fAccepted = kFALSE;
			acceptedDet    = kFALSE;
		    }
		    fhistory.push_back(part);
		}
	    }
	}
	if(!felements[felements.size()-1].fAccepted) accepted = kFALSE;

	return  accepted&acceptedDet;
    }

    vector <HBeamParticle>& newParticle()
    {
	HBeamParticle part(fBeam);
	fhistory.clear();

	createBeam(part);         // beam profile + smearing

	fhistory.push_back(part);

	transport(part);          // transport particle through beamline

	Bool_t accepted = createDetectorHits(part); // fill detectors

	// calc particle at 0,0,0 ideal target
	Double_t* out = &felements[fnum_target].fout[0];
	Double_t p,px,py,pz;
	p  = fBeam.fBeamMom * (1. + out[4]);
	pz = p / sqrt(1.+ tan(out[1])*tan(out[1]) + tan(out[3])*tan(out[3]));
	px = pz * tan(out[1]);
	py = pz * tan(out[3]);

	part.fP  .SetXYZ(px,py,pz);
	part.fPos.SetXYZ(out[0],out[2],felements[fnum_target].fDistance);
	part.fDetnum = -1;
	part.fName   = felements[fnum_target].fName;
	part.fAccepted = accepted;
	fhistory.push_back(part);

	return fhistory;
    }


    void printBeamLine(Bool_t trans=kFALSE)
    {
	cout<<"############################################################################"<<endl;
	cout<<"BEAMLINE ELEMENTS :"<<endl;
	for(UInt_t i = 0; i < felements.size(); i ++){
	    felements[i].print(trans);
	}
	cout<<"############################################################################"<<endl;
    }
    void printBeamProfile()
    {
	cout<<"############################################################################"<<endl;
	cout<<"BEAM PROFILE :"<<endl;
	for(Int_t i = 0 ; i < fProfile->GetNpar();i ++){
	    cout<<setw(20)<<fProfile->GetParName(i)<<" : "<<setw(12)<<fProfile->GetParameter(i)<<endl;
	}
	fBeam.print(kFALSE);
	cout<<"############################################################################"<<endl;
    }
    void printDetectors()
    {
	cout<<"############################################################################"<<endl;
	cout<<"DETECTORS :"<<endl;
	for(UInt_t i = 0; i < fdetectors.size(); i ++){
	    fdetectors[i].print(kFALSE);
	}
	cout<<"############################################################################"<<endl;
    }

    Bool_t initBeamLine(TString filename,Int_t targetElementNum,Bool_t debug=kFALSE)
    {
	//opens the ASCII description of the beam line


	cout<<"HBeam:: reading input!"<<endl;
	Int_t nelem = 0;



	std::ifstream input;
	input.open(filename.Data());

	if(input.fail()){
	    cout<<"HBeam:: could not open input!"<<endl;
	    return kFALSE;
	}


	Char_t str[1000];
	Int_t  nread = 1000;
	Bool_t readfile = kTRUE;
	//----------------------------
	// local dummy vars
	Double_t aa;
	Int_t ind1;
	Int_t ind23;
	Double_t distance;
	//----------------------------


	while (readfile && !input.fail()) {

	    input>>distance;
	    if(input.peek()!='\n') input.ignore(nread,'\n');
	    if(input.fail() && nelem > 0 ) continue;
	    if(input.fail() ){ cout<<"HBeam:: "<<nelem<<" failed reading distance"<<endl; return kFALSE;}
	    distance *= 1000;  // [m] -> [mm]
	    if(debug) cout<<nelem<<" ------------------------"<<distance<<endl;
	    for(Int_t i = 0; i < 5; i++) {
		input.ignore(nread,'\n');
	    }
	    if(input.fail() ) {
		cout<<"HBeam:: "<<nelem<<" failed after first block "<<endl;
		return kFALSE;
	    }
	    input.getline(str,nread,'\n');  // *TRANSFORM 1*
	    if(debug) cout<<"HBeam:: "<<"1st "<<str<<endl;

	    HBeamElement e(Form("Element[%i]",nelem),distance,nelem);
	    // first order transform
	    Int_t ind = 0;
	    for(Int_t i = 0; i < 6; i++){ // lines
		if(i == 4) {
		    input>>aa>>aa>>aa>>aa>>aa>>aa;
		    continue;
		}
		if(i > 4) ind = i-1;
		else      ind = i;
		input>>e.Tij[ind][0]>>e.Tij[ind][1]>>e.Tij[ind][2]>>e.Tij[ind][3]>>aa>>e.Tij[ind][4];
	    }
	    if(debug) e.printTij();

	    if(input.fail() ) {
		cout<<"HBeam:: "<<nelem<<" failed after first order "<<endl;
		return kFALSE;
	    }
	    if(input.peek() =='\n') input.ignore(nread,'\n');
	    input.getline(str,nread,'\n');   // *2ND ORDER TRANSFORM*

	    if(debug) cout<<"HBeam:: "<<"2nd : "<<str<<endl;
	    //C
	    //C     second order transform
	    //C

	    for(Int_t i = 0; i < 5; i++){ // blocks i
		input>>ind1>>ind23>>e.Tijk[i][0][0];
		input>>ind1>>ind23>>e.Tijk[i][0][1]>>ind1>>ind23>>e.Tijk[i][1][1];
		input>>ind1>>ind23>>e.Tijk[i][0][2]>>ind1>>ind23>>e.Tijk[i][1][2]>>ind1>>ind23>>e.Tijk[i][2][2];
		input>>ind1>>ind23>>e.Tijk[i][0][3]>>ind1>>ind23>>e.Tijk[i][1][3]>>ind1>>ind23>>e.Tijk[i][2][3]>>ind1>>ind23>>e.Tijk[i][3][3];
		input>>ind1>>ind23>>aa             >>ind1>>ind23>>aa             >>ind1>>ind23>>aa             >>ind1>>ind23>>aa             >>ind1>>ind23>>aa;
		input>>ind1>>ind23>>e.Tijk[i][0][4]>>ind1>>ind23>>e.Tijk[i][1][4]>>ind1>>ind23>>e.Tijk[i][2][4]>>ind1>>ind23>>e.Tijk[i][3][4]>>ind1>>ind23>>aa>>ind1>>ind23>>e.Tijk[i][4][4];
	    }

	    if(debug) e.printTijk();

	    if(input.fail() ) {
		cout<<"HBeam:: "<<nelem<<" failed after second order "<<endl;
		return kFALSE;
	    }

	    felements.push_back(e);
	    input.ignore(nread,'\n');

	    if(input.fail() ) return kFALSE;

	    if (!input.eof())
		nelem++;
	    else
		readfile = kFALSE;

	} //end nelem

	input.close();
	cout<<"HBeam:: finished reading input!"<<endl;
	return setTargetElement(targetElementNum);

    }



#ifndef isPluto
    ClassDef(HBeam,1)
#endif

};


#endif


 HBeam.h:1
 HBeam.h:2
 HBeam.h:3
 HBeam.h:4
 HBeam.h:5
 HBeam.h:6
 HBeam.h:7
 HBeam.h:8
 HBeam.h:9
 HBeam.h:10
 HBeam.h:11
 HBeam.h:12
 HBeam.h:13
 HBeam.h:14
 HBeam.h:15
 HBeam.h:16
 HBeam.h:17
 HBeam.h:18
 HBeam.h:19
 HBeam.h:20
 HBeam.h:21
 HBeam.h:22
 HBeam.h:23
 HBeam.h:24
 HBeam.h:25
 HBeam.h:26
 HBeam.h:27
 HBeam.h:28
 HBeam.h:29
 HBeam.h:30
 HBeam.h:31
 HBeam.h:32
 HBeam.h:33
 HBeam.h:34
 HBeam.h:35
 HBeam.h:36
 HBeam.h:37
 HBeam.h:38
 HBeam.h:39
 HBeam.h:40
 HBeam.h:41
 HBeam.h:42
 HBeam.h:43
 HBeam.h:44
 HBeam.h:45
 HBeam.h:46
 HBeam.h:47
 HBeam.h:48
 HBeam.h:49
 HBeam.h:50
 HBeam.h:51
 HBeam.h:52
 HBeam.h:53
 HBeam.h:54
 HBeam.h:55
 HBeam.h:56
 HBeam.h:57
 HBeam.h:58
 HBeam.h:59
 HBeam.h:60
 HBeam.h:61
 HBeam.h:62
 HBeam.h:63
 HBeam.h:64
 HBeam.h:65
 HBeam.h:66
 HBeam.h:67
 HBeam.h:68
 HBeam.h:69
 HBeam.h:70
 HBeam.h:71
 HBeam.h:72
 HBeam.h:73
 HBeam.h:74
 HBeam.h:75
 HBeam.h:76
 HBeam.h:77
 HBeam.h:78
 HBeam.h:79
 HBeam.h:80
 HBeam.h:81
 HBeam.h:82
 HBeam.h:83
 HBeam.h:84
 HBeam.h:85
 HBeam.h:86
 HBeam.h:87
 HBeam.h:88
 HBeam.h:89
 HBeam.h:90
 HBeam.h:91
 HBeam.h:92
 HBeam.h:93
 HBeam.h:94
 HBeam.h:95
 HBeam.h:96
 HBeam.h:97
 HBeam.h:98
 HBeam.h:99
 HBeam.h:100
 HBeam.h:101
 HBeam.h:102
 HBeam.h:103
 HBeam.h:104
 HBeam.h:105
 HBeam.h:106
 HBeam.h:107
 HBeam.h:108
 HBeam.h:109
 HBeam.h:110
 HBeam.h:111
 HBeam.h:112
 HBeam.h:113
 HBeam.h:114
 HBeam.h:115
 HBeam.h:116
 HBeam.h:117
 HBeam.h:118
 HBeam.h:119
 HBeam.h:120
 HBeam.h:121
 HBeam.h:122
 HBeam.h:123
 HBeam.h:124
 HBeam.h:125
 HBeam.h:126
 HBeam.h:127
 HBeam.h:128
 HBeam.h:129
 HBeam.h:130
 HBeam.h:131
 HBeam.h:132
 HBeam.h:133
 HBeam.h:134
 HBeam.h:135
 HBeam.h:136
 HBeam.h:137
 HBeam.h:138
 HBeam.h:139
 HBeam.h:140
 HBeam.h:141
 HBeam.h:142
 HBeam.h:143
 HBeam.h:144
 HBeam.h:145
 HBeam.h:146
 HBeam.h:147
 HBeam.h:148
 HBeam.h:149
 HBeam.h:150
 HBeam.h:151
 HBeam.h:152
 HBeam.h:153
 HBeam.h:154
 HBeam.h:155
 HBeam.h:156
 HBeam.h:157
 HBeam.h:158
 HBeam.h:159
 HBeam.h:160
 HBeam.h:161
 HBeam.h:162
 HBeam.h:163
 HBeam.h:164
 HBeam.h:165
 HBeam.h:166
 HBeam.h:167
 HBeam.h:168
 HBeam.h:169
 HBeam.h:170
 HBeam.h:171
 HBeam.h:172
 HBeam.h:173
 HBeam.h:174
 HBeam.h:175
 HBeam.h:176
 HBeam.h:177
 HBeam.h:178
 HBeam.h:179
 HBeam.h:180
 HBeam.h:181
 HBeam.h:182
 HBeam.h:183
 HBeam.h:184
 HBeam.h:185
 HBeam.h:186
 HBeam.h:187
 HBeam.h:188
 HBeam.h:189
 HBeam.h:190
 HBeam.h:191
 HBeam.h:192
 HBeam.h:193
 HBeam.h:194
 HBeam.h:195
 HBeam.h:196
 HBeam.h:197
 HBeam.h:198
 HBeam.h:199
 HBeam.h:200
 HBeam.h:201
 HBeam.h:202
 HBeam.h:203
 HBeam.h:204
 HBeam.h:205
 HBeam.h:206
 HBeam.h:207
 HBeam.h:208
 HBeam.h:209
 HBeam.h:210
 HBeam.h:211
 HBeam.h:212
 HBeam.h:213
 HBeam.h:214
 HBeam.h:215
 HBeam.h:216
 HBeam.h:217
 HBeam.h:218
 HBeam.h:219
 HBeam.h:220
 HBeam.h:221
 HBeam.h:222
 HBeam.h:223
 HBeam.h:224
 HBeam.h:225
 HBeam.h:226
 HBeam.h:227
 HBeam.h:228
 HBeam.h:229
 HBeam.h:230
 HBeam.h:231
 HBeam.h:232
 HBeam.h:233
 HBeam.h:234
 HBeam.h:235
 HBeam.h:236
 HBeam.h:237
 HBeam.h:238
 HBeam.h:239
 HBeam.h:240
 HBeam.h:241
 HBeam.h:242
 HBeam.h:243
 HBeam.h:244
 HBeam.h:245
 HBeam.h:246
 HBeam.h:247
 HBeam.h:248
 HBeam.h:249
 HBeam.h:250
 HBeam.h:251
 HBeam.h:252
 HBeam.h:253
 HBeam.h:254
 HBeam.h:255
 HBeam.h:256
 HBeam.h:257
 HBeam.h:258
 HBeam.h:259
 HBeam.h:260
 HBeam.h:261
 HBeam.h:262
 HBeam.h:263
 HBeam.h:264
 HBeam.h:265
 HBeam.h:266
 HBeam.h:267
 HBeam.h:268
 HBeam.h:269
 HBeam.h:270
 HBeam.h:271
 HBeam.h:272
 HBeam.h:273
 HBeam.h:274
 HBeam.h:275
 HBeam.h:276
 HBeam.h:277
 HBeam.h:278
 HBeam.h:279
 HBeam.h:280
 HBeam.h:281
 HBeam.h:282
 HBeam.h:283
 HBeam.h:284
 HBeam.h:285
 HBeam.h:286
 HBeam.h:287
 HBeam.h:288
 HBeam.h:289
 HBeam.h:290
 HBeam.h:291
 HBeam.h:292
 HBeam.h:293
 HBeam.h:294
 HBeam.h:295
 HBeam.h:296
 HBeam.h:297
 HBeam.h:298
 HBeam.h:299
 HBeam.h:300
 HBeam.h:301
 HBeam.h:302
 HBeam.h:303
 HBeam.h:304
 HBeam.h:305
 HBeam.h:306
 HBeam.h:307
 HBeam.h:308
 HBeam.h:309
 HBeam.h:310
 HBeam.h:311
 HBeam.h:312
 HBeam.h:313
 HBeam.h:314
 HBeam.h:315
 HBeam.h:316
 HBeam.h:317
 HBeam.h:318
 HBeam.h:319
 HBeam.h:320
 HBeam.h:321
 HBeam.h:322
 HBeam.h:323
 HBeam.h:324
 HBeam.h:325
 HBeam.h:326
 HBeam.h:327
 HBeam.h:328
 HBeam.h:329
 HBeam.h:330
 HBeam.h:331
 HBeam.h:332
 HBeam.h:333
 HBeam.h:334
 HBeam.h:335
 HBeam.h:336
 HBeam.h:337
 HBeam.h:338
 HBeam.h:339
 HBeam.h:340
 HBeam.h:341
 HBeam.h:342
 HBeam.h:343
 HBeam.h:344
 HBeam.h:345
 HBeam.h:346
 HBeam.h:347
 HBeam.h:348
 HBeam.h:349
 HBeam.h:350
 HBeam.h:351
 HBeam.h:352
 HBeam.h:353
 HBeam.h:354
 HBeam.h:355
 HBeam.h:356
 HBeam.h:357
 HBeam.h:358
 HBeam.h:359
 HBeam.h:360
 HBeam.h:361
 HBeam.h:362
 HBeam.h:363
 HBeam.h:364
 HBeam.h:365
 HBeam.h:366
 HBeam.h:367
 HBeam.h:368
 HBeam.h:369
 HBeam.h:370
 HBeam.h:371
 HBeam.h:372
 HBeam.h:373
 HBeam.h:374
 HBeam.h:375
 HBeam.h:376
 HBeam.h:377
 HBeam.h:378
 HBeam.h:379
 HBeam.h:380
 HBeam.h:381
 HBeam.h:382
 HBeam.h:383
 HBeam.h:384
 HBeam.h:385
 HBeam.h:386
 HBeam.h:387
 HBeam.h:388
 HBeam.h:389
 HBeam.h:390
 HBeam.h:391
 HBeam.h:392
 HBeam.h:393
 HBeam.h:394
 HBeam.h:395
 HBeam.h:396
 HBeam.h:397
 HBeam.h:398
 HBeam.h:399
 HBeam.h:400
 HBeam.h:401
 HBeam.h:402
 HBeam.h:403
 HBeam.h:404
 HBeam.h:405
 HBeam.h:406
 HBeam.h:407
 HBeam.h:408
 HBeam.h:409
 HBeam.h:410
 HBeam.h:411
 HBeam.h:412
 HBeam.h:413
 HBeam.h:414
 HBeam.h:415
 HBeam.h:416
 HBeam.h:417
 HBeam.h:418
 HBeam.h:419
 HBeam.h:420
 HBeam.h:421
 HBeam.h:422
 HBeam.h:423
 HBeam.h:424
 HBeam.h:425
 HBeam.h:426
 HBeam.h:427
 HBeam.h:428
 HBeam.h:429
 HBeam.h:430
 HBeam.h:431
 HBeam.h:432
 HBeam.h:433
 HBeam.h:434
 HBeam.h:435
 HBeam.h:436
 HBeam.h:437
 HBeam.h:438
 HBeam.h:439
 HBeam.h:440
 HBeam.h:441
 HBeam.h:442
 HBeam.h:443
 HBeam.h:444
 HBeam.h:445
 HBeam.h:446
 HBeam.h:447
 HBeam.h:448
 HBeam.h:449
 HBeam.h:450
 HBeam.h:451
 HBeam.h:452
 HBeam.h:453
 HBeam.h:454
 HBeam.h:455
 HBeam.h:456
 HBeam.h:457
 HBeam.h:458
 HBeam.h:459
 HBeam.h:460
 HBeam.h:461
 HBeam.h:462
 HBeam.h:463
 HBeam.h:464
 HBeam.h:465
 HBeam.h:466
 HBeam.h:467
 HBeam.h:468
 HBeam.h:469
 HBeam.h:470
 HBeam.h:471
 HBeam.h:472
 HBeam.h:473
 HBeam.h:474
 HBeam.h:475
 HBeam.h:476
 HBeam.h:477
 HBeam.h:478
 HBeam.h:479
 HBeam.h:480
 HBeam.h:481
 HBeam.h:482
 HBeam.h:483
 HBeam.h:484
 HBeam.h:485
 HBeam.h:486
 HBeam.h:487
 HBeam.h:488
 HBeam.h:489
 HBeam.h:490
 HBeam.h:491
 HBeam.h:492
 HBeam.h:493
 HBeam.h:494
 HBeam.h:495
 HBeam.h:496
 HBeam.h:497
 HBeam.h:498
 HBeam.h:499
 HBeam.h:500
 HBeam.h:501
 HBeam.h:502
 HBeam.h:503
 HBeam.h:504
 HBeam.h:505
 HBeam.h:506
 HBeam.h:507
 HBeam.h:508
 HBeam.h:509
 HBeam.h:510
 HBeam.h:511
 HBeam.h:512
 HBeam.h:513
 HBeam.h:514
 HBeam.h:515
 HBeam.h:516
 HBeam.h:517
 HBeam.h:518
 HBeam.h:519
 HBeam.h:520
 HBeam.h:521
 HBeam.h:522
 HBeam.h:523
 HBeam.h:524
 HBeam.h:525
 HBeam.h:526
 HBeam.h:527
 HBeam.h:528
 HBeam.h:529
 HBeam.h:530
 HBeam.h:531
 HBeam.h:532
 HBeam.h:533
 HBeam.h:534
 HBeam.h:535
 HBeam.h:536
 HBeam.h:537
 HBeam.h:538
 HBeam.h:539
 HBeam.h:540
 HBeam.h:541
 HBeam.h:542
 HBeam.h:543
 HBeam.h:544
 HBeam.h:545
 HBeam.h:546
 HBeam.h:547
 HBeam.h:548
 HBeam.h:549
 HBeam.h:550
 HBeam.h:551
 HBeam.h:552
 HBeam.h:553
 HBeam.h:554
 HBeam.h:555
 HBeam.h:556
 HBeam.h:557
 HBeam.h:558
 HBeam.h:559
 HBeam.h:560
 HBeam.h:561
 HBeam.h:562
 HBeam.h:563
 HBeam.h:564
 HBeam.h:565
 HBeam.h:566
 HBeam.h:567
 HBeam.h:568
 HBeam.h:569
 HBeam.h:570
 HBeam.h:571
 HBeam.h:572
 HBeam.h:573
 HBeam.h:574
 HBeam.h:575
 HBeam.h:576
 HBeam.h:577
 HBeam.h:578
 HBeam.h:579
 HBeam.h:580
 HBeam.h:581
 HBeam.h:582
 HBeam.h:583
 HBeam.h:584
 HBeam.h:585
 HBeam.h:586
 HBeam.h:587
 HBeam.h:588
 HBeam.h:589
 HBeam.h:590
 HBeam.h:591
 HBeam.h:592
 HBeam.h:593
 HBeam.h:594
 HBeam.h:595
 HBeam.h:596
 HBeam.h:597
 HBeam.h:598
 HBeam.h:599
 HBeam.h:600
 HBeam.h:601
 HBeam.h:602
 HBeam.h:603
 HBeam.h:604
 HBeam.h:605
 HBeam.h:606
 HBeam.h:607
 HBeam.h:608
 HBeam.h:609
 HBeam.h:610
 HBeam.h:611
 HBeam.h:612
 HBeam.h:613
 HBeam.h:614
 HBeam.h:615
 HBeam.h:616
 HBeam.h:617
 HBeam.h:618
 HBeam.h:619
 HBeam.h:620
 HBeam.h:621
 HBeam.h:622
 HBeam.h:623
 HBeam.h:624
 HBeam.h:625
 HBeam.h:626
 HBeam.h:627
 HBeam.h:628
 HBeam.h:629
 HBeam.h:630
 HBeam.h:631
 HBeam.h:632
 HBeam.h:633
 HBeam.h:634
 HBeam.h:635
 HBeam.h:636
 HBeam.h:637
 HBeam.h:638
 HBeam.h:639
 HBeam.h:640
 HBeam.h:641
 HBeam.h:642
 HBeam.h:643
 HBeam.h:644
 HBeam.h:645
 HBeam.h:646
 HBeam.h:647
 HBeam.h:648
 HBeam.h:649
 HBeam.h:650
 HBeam.h:651
 HBeam.h:652
 HBeam.h:653
 HBeam.h:654
 HBeam.h:655
 HBeam.h:656
 HBeam.h:657
 HBeam.h:658
 HBeam.h:659
 HBeam.h:660
 HBeam.h:661
 HBeam.h:662
 HBeam.h:663
 HBeam.h:664
 HBeam.h:665
 HBeam.h:666
 HBeam.h:667
 HBeam.h:668
 HBeam.h:669
 HBeam.h:670
 HBeam.h:671
 HBeam.h:672
 HBeam.h:673
 HBeam.h:674
 HBeam.h:675
 HBeam.h:676
 HBeam.h:677
 HBeam.h:678
 HBeam.h:679
 HBeam.h:680
 HBeam.h:681
 HBeam.h:682
 HBeam.h:683
 HBeam.h:684
 HBeam.h:685
 HBeam.h:686
 HBeam.h:687
 HBeam.h:688
 HBeam.h:689
 HBeam.h:690
 HBeam.h:691
 HBeam.h:692
 HBeam.h:693
 HBeam.h:694
 HBeam.h:695
 HBeam.h:696
 HBeam.h:697
 HBeam.h:698
 HBeam.h:699
 HBeam.h:700
 HBeam.h:701
 HBeam.h:702
 HBeam.h:703
 HBeam.h:704
 HBeam.h:705
 HBeam.h:706
 HBeam.h:707
 HBeam.h:708
 HBeam.h:709
 HBeam.h:710
 HBeam.h:711
 HBeam.h:712
 HBeam.h:713
 HBeam.h:714
 HBeam.h:715
 HBeam.h:716
 HBeam.h:717
 HBeam.h:718
 HBeam.h:719
 HBeam.h:720
 HBeam.h:721
 HBeam.h:722
 HBeam.h:723
 HBeam.h:724
 HBeam.h:725
 HBeam.h:726
 HBeam.h:727
 HBeam.h:728
 HBeam.h:729
 HBeam.h:730
 HBeam.h:731
 HBeam.h:732
 HBeam.h:733
 HBeam.h:734
 HBeam.h:735
 HBeam.h:736
 HBeam.h:737