ROOT logo
//*--Author A.Rustamov
//
//Input position vector,output field vector
//In LAB coordinate system
//
//               |Y             
//               |
//        ***************
//         *     |     * 
//          *    |    *
//           *   |   *
//            *  |  *
//             * | *
//              ***
//               | 
//   -------------             
//   X
//
// UNITS: TESLA,cm, GeV for track momentum,
// polarity is reverse compared to HYDRA (defined by Eve ..)
// switching is done automatically

#include "hedfield.h"
#include "hades.h"
#include "hruntimedb.h"
#include "hmdctrackgfieldpar.h"

#include "TEveVSDStructs.h"
#include "TSystem.h"
#include "TMath.h"
#include <math.h>
#include <unistd.h>
#include <iostream>
#include <iomanip>
#include <fstream>

using namespace std;

ClassImp(HEDField)

HEDField::HEDField(const Char_t *name, const Char_t *title,Double_t fPolarity)
    : TNamed(name,title) , TEveMagField()
{

    //----------------------------------
    // eve
    gEDField=this;
    fFieldConstant = kFALSE;
    //----------------------------------

    fpol  = -fPolarity;
    dconv = 0.0001; // 0.0001 [Tesla]

    p_tzfl  = new Double_t[708400]; //161*176*25
    p_tpfl  = new Double_t[708400];
    p_trfl  = new Double_t[708400];
    nfz     = 161;
    nfr     = 176;
    nfp     = 25;
    nfz_nfr = nfz*nfr;

    one_sixtyth = 1.0 / 60.0;
    
}

HEDField::~HEDField()
{
    delete [] p_tzfl;
    delete [] p_trfl;
    delete [] p_tpfl;
}
void HEDField::Clear()
{
}

void HEDField::ReadAscii(TString infile)
{

    if(gSystem->AccessPathName(infile.Data()) != 0){
	Error("ReadAscii()","Could not find fieldmap %s !",infile.Data());
        exit(1);
    }

    ifstream input;
    input.open(infile.Data(),ios::in);
    while(!input.eof())
    {

	input>>zflmin ;
	input>>zflmax ;
	input>>zfldel;
	input>>rflmin ;
	input>>rflmax ;
	input>>rfldel ;
	input>>pflmin ;
	input>>pflmax ;
	input>>pfldel ;
	for (Int_t i1 = 0; i1 < nfp; i1 ++){
	    for (Int_t i2 = 0; i2 < nfr; i2 ++){
		for (Int_t i3 = 0;i3 < nfz; i3 ++){input>>p_tzfl[i3 + i2 * nfz + i1 * nfz * nfr];}
		for (Int_t i3 = 0;i3 < nfz; i3 ++){input>>p_trfl[i3 + i2 * nfz + i1 * nfz * nfr];}
		for (Int_t i3 = 0;i3 < nfz; i3 ++){input>>p_tpfl[i3 + i2 * nfz + i1 * nfz * nfr];}
	    }
	}
    }
    cout<<"HEDField::READING FIELDMAP FINISHED!"<<endl;
    input.close();
    step1z = 1.0 / zfldel;
    step1r = 1.0 / rfldel;
    step1p = 1.0 / pfldel;
}

void HEDField::WriteAscii(TString outfile)
{
    ofstream output;
    output.open(outfile.Data(),ios::out);

    output<<zflmin<<endl;
    output<<zflmax<<endl;
    output<<zfldel<<endl;
    output<<rflmin<<endl;
    output<<rflmax<<endl;
    output<<rfldel<<endl;
    output<<pflmin<<endl;
    output<<pflmax<<endl;
    output<<pfldel<<endl;
    Int_t ct = 0;
    for (Int_t i1 = 0; i1 < nfp; i1 ++){
	for (Int_t i2 = 0; i2 < nfr; i2 ++){
	    for (Int_t i3 = 0;i3 < nfz; i3 ++){ct++; output<<" "<<setw(13)<<p_tzfl[i3 + i2 * nfz + i1 * nfz * nfr]; if(ct%10 == 0) output<<endl; }
	    for (Int_t i3 = 0;i3 < nfz; i3 ++){ct++; output<<" "<<setw(13)<<p_trfl[i3 + i2 * nfz + i1 * nfz * nfr]; if(ct%10 == 0) output<<endl; }
	    for (Int_t i3 = 0;i3 < nfz; i3 ++){ct++; output<<" "<<setw(13)<<p_tpfl[i3 + i2 * nfz + i1 * nfz * nfr]; if(ct%10 == 0) output<<endl; }
	}
    }
    cout<<"HEDField::WRITING FIELDMAP FINISHED!"<<endl;
    output.close();
}
	  	  
void HEDField::CalcField(Double_t *xv,Double_t *btos) const
{
    Double_t xloc,yloc,zloc,rhog,phigd,phil,bz,br,bp;
    Double_t m1,m2,m3;
    Double_t dc,ds,a,b,c,ac,ab,bc,abc,w1,w2,w3,w4,w5,w6,w7,w8;
    Double_t spol;
    Double_t mult;
    Int_t ifz,ifr,ifp,ifr1,ifp1;
    Int_t i1,i2,i3,i4,i5,i6,i7,i8;
    Double_t phig;
    zloc = xv[2];

    if(zloc < zflmin || zloc > zflmax) {
	btos[0] = 0.; btos[1] = 0.;btos[2] = 0.;
	return;
    }
    xloc = xv[0];
    yloc = xv[1];
    rhog = sqrt(xloc*xloc+yloc*yloc);

    if(rhog > rflmax){
	btos[0] = 0.;
	btos[1] = 0.;
	btos[2] = 0.;
	return;
    }
    if(rhog > 0.01){

	phig = atan2(yloc,xloc);
	if(phig < 0.) phig = phig + TMath::Pi() * 2;
	phigd = phig * TMath::RadToDeg();
	dc = TMath::Cos(phig);
	ds = TMath::Sin(phig);

    } else {
	rhog  = 0.;
	phigd = 0.;
	dc    = 0.;
	ds    = 0.;
    }

    phil=fmod(phigd,60.);

    if(phil > 30.0){
	spol = -1.0;
	phil = 60.0 - phil;
    } else {
	spol = 1.0;
    }


    m1 = (zloc - zflmin) * step1z;
    m2 = (rhog - rflmin) * step1r;
    m3 = (phil - pflmin) * step1p;

    ifz = (Int_t)(m1);
    ifr = (Int_t)(m2);
    ifp = (Int_t)(m3);

    a = m1 - ifz;
    b = m2 - ifr;
    c = m3 - ifp;

    ifz ++;
    ifr  = (ifr + 1) * nfz;
    ifp  = (ifp + 1) * nfz_nfr;
    ifr1 = ifr - nfz;
    ifp1 = ifp - nfz_nfr;

    ab  = a * b;
    ac  = a * c;
    bc  = b * c;
    abc = a * bc;
    w1  = 1.0 - a - b + ab - c + ac + bc - abc;
    w2  = b - ab - bc + abc;
    w3  = bc - abc;
    w4  = c - ac - bc + abc;
    w5  = a - ab - ac + abc;
    w6  = ab - abc;
    w7  = abc;
    w8  = ac - abc;

    i1 = ifz -1+ ifr1 + ifp1;
    i2 = ifz -1+ ifr  + ifp1;
    i3 = ifz -1+ ifr  + ifp;
    i4 = ifz -1+ ifr1 + ifp;
    i5 = ifz   + ifr1 + ifp1;
    i6 = ifz   + ifr  + ifp1;
    i7 = ifz   + ifr  + ifp;
    i8 = ifz   + ifr1 + ifp;

    bz =  w1 * p_tzfl[i1] + w2 * p_tzfl[i2]
	+ w3 * p_tzfl[i3] + w4 * p_tzfl[i4]
	+ w5 * p_tzfl[i5] + w6 * p_tzfl[i6]
	+ w7 * p_tzfl[i7] + w8 * p_tzfl[i8];

    br =  w1 * p_trfl[i1] + w2 * p_trfl[i2]
	+ w3 * p_trfl[i3] + w4 * p_trfl[i4]
	+ w5 * p_trfl[i5] + w6 * p_trfl[i6]
	+ w7 * p_trfl[i7] + w8 * p_trfl[i8];

    bp =  w1 * p_tpfl[i1] + w2 * p_tpfl[i2]
	+ w3 * p_tpfl[i3] + w4 * p_tpfl[i4]
	+ w5 * p_tpfl[i5] + w6 * p_tpfl[i6]
	+ w7 * p_tpfl[i7] + w8 * p_tpfl[i8];



    mult = fpol * dconv;
    bz = -bz * spol * mult;
    br = -br * spol * mult;
    bp = -bp * mult;

    btos[0] = (dc * br - ds * bp);
    btos[1] = (ds * br + dc * bp);
    btos[2] = bz;

}

TEveVector HEDField::GetField(Float_t x, Float_t y, Float_t z) const
{

    Double_t xv[3] = {x,y,z}; // cm input from GEANT/EVENTDISPLAY
    Double_t b [3];
    CalcField(xv,b);

    TEveVector val((Float_t)b[0],(Float_t)b[1],(Float_t)b[2]);

    return val;

}

TEveVector  HEDField::GetField(const TEveVector& v) const { return HEDField::GetField(v.fX,v.fY,v.fZ);};


void HEDField::Streamer(TBuffer &R__b) {
  // Stream an object of class HEDField.
  // The data in the 1-dimensional arrays are streamed into/from the old
  // 3-dimensional arrays (no change of class version)
  UInt_t R__s, R__c;
  if (R__b.IsReading()) {
    Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
    TNamed::Streamer(R__b);
    R__b >> nfz;
    R__b >> nfr;
    R__b >> nfp;
    R__b >> zflmin;
    R__b >> zflmax;
    R__b >> zfldel;
    R__b >> rflmin;
    R__b >> rflmax;
    R__b >> rfldel;
    R__b >> pflmin;
    R__b >> pflmax;
    R__b >> pfldel;
    R__b >> fpol;
      
    Int_t n;
    R__b >> n;
    for (Int_t i3 = 0; i3 < nfz; i3 ++)
    for (Int_t i2 = 0; i2 < nfr; i2 ++)  
    for (Int_t i1 = 0; i1 < nfp; i1 ++) 
      R__b >> p_tzfl[i3+i2*nfz+i1*nfz*nfr];

    R__b >> n;
    for (Int_t i3 = 0; i3 < nfz; i3 ++)
    for (Int_t i2 = 0; i2 < nfr; i2 ++)  
    for (Int_t i1 = 0; i1 < nfp; i1 ++) 
      R__b >> p_trfl[i3+i2*nfz+i1*nfz*nfr];

    R__b >> n;
    for (Int_t i3 = 0; i3 < nfz; i3 ++)
    for (Int_t i2 = 0; i2 < nfr; i2 ++)  
    for (Int_t i1 = 0; i1 < nfp; i1 ++) 
      R__b >> p_tpfl[i3+i2*nfz+i1*nfz*nfr];

    step1z=1.0/zfldel;
    step1r=1.0/rfldel;
    step1p=1.0/pfldel;

    fFieldConstant = kFALSE;

    R__b.CheckByteCount(R__s, R__c, HEDField::IsA());
      
  } else {
    R__c = R__b.WriteVersion(HEDField::IsA(), kTRUE);
    TNamed::Streamer(R__b);
    R__b << nfz;
    R__b << nfr;
    R__b << nfp;
    R__b << zflmin;
    R__b << zflmax;
    R__b << zfldel;
    R__b << rflmin;
    R__b << rflmax;
    R__b << rfldel;
    R__b << pflmin;
    R__b << pflmax;
    R__b << pfldel;
    R__b << fpol;

    Int_t n=nfz*nfr*nfp;
    R__b << n;
    for (Int_t i3 = 0; i3 < nfz; i3 ++)
    for (Int_t i2 = 0; i2 < nfr; i2 ++)  
    for (Int_t i1 = 0; i1 < nfp; i1 ++) 
      R__b << p_tzfl[i3+i2*nfz+i1*nfz*nfr];

    R__b << n;
    for (Int_t i3 = 0; i3 < nfz; i3 ++)
    for (Int_t i2 = 0; i2 < nfr; i2 ++)  
    for (Int_t i1 = 0; i1 < nfp; i1 ++) 
      R__b << p_trfl[i3+i2*nfz+i1*nfz*nfr];

    R__b << n;
    for (Int_t i3 = 0; i3 < nfz; i3 ++)
    for (Int_t i2 = 0; i2 < nfr; i2 ++)  
    for (Int_t i1 = 0; i1 < nfp; i1 ++) 
      R__b << p_tpfl[i3+i2*nfz+i1*nfz*nfr];

    R__b.SetByteCount(R__c, kTRUE);
  }
}

HEDField *gEDField;

 hedfield.cc:1
 hedfield.cc:2
 hedfield.cc:3
 hedfield.cc:4
 hedfield.cc:5
 hedfield.cc:6
 hedfield.cc:7
 hedfield.cc:8
 hedfield.cc:9
 hedfield.cc:10
 hedfield.cc:11
 hedfield.cc:12
 hedfield.cc:13
 hedfield.cc:14
 hedfield.cc:15
 hedfield.cc:16
 hedfield.cc:17
 hedfield.cc:18
 hedfield.cc:19
 hedfield.cc:20
 hedfield.cc:21
 hedfield.cc:22
 hedfield.cc:23
 hedfield.cc:24
 hedfield.cc:25
 hedfield.cc:26
 hedfield.cc:27
 hedfield.cc:28
 hedfield.cc:29
 hedfield.cc:30
 hedfield.cc:31
 hedfield.cc:32
 hedfield.cc:33
 hedfield.cc:34
 hedfield.cc:35
 hedfield.cc:36
 hedfield.cc:37
 hedfield.cc:38
 hedfield.cc:39
 hedfield.cc:40
 hedfield.cc:41
 hedfield.cc:42
 hedfield.cc:43
 hedfield.cc:44
 hedfield.cc:45
 hedfield.cc:46
 hedfield.cc:47
 hedfield.cc:48
 hedfield.cc:49
 hedfield.cc:50
 hedfield.cc:51
 hedfield.cc:52
 hedfield.cc:53
 hedfield.cc:54
 hedfield.cc:55
 hedfield.cc:56
 hedfield.cc:57
 hedfield.cc:58
 hedfield.cc:59
 hedfield.cc:60
 hedfield.cc:61
 hedfield.cc:62
 hedfield.cc:63
 hedfield.cc:64
 hedfield.cc:65
 hedfield.cc:66
 hedfield.cc:67
 hedfield.cc:68
 hedfield.cc:69
 hedfield.cc:70
 hedfield.cc:71
 hedfield.cc:72
 hedfield.cc:73
 hedfield.cc:74
 hedfield.cc:75
 hedfield.cc:76
 hedfield.cc:77
 hedfield.cc:78
 hedfield.cc:79
 hedfield.cc:80
 hedfield.cc:81
 hedfield.cc:82
 hedfield.cc:83
 hedfield.cc:84
 hedfield.cc:85
 hedfield.cc:86
 hedfield.cc:87
 hedfield.cc:88
 hedfield.cc:89
 hedfield.cc:90
 hedfield.cc:91
 hedfield.cc:92
 hedfield.cc:93
 hedfield.cc:94
 hedfield.cc:95
 hedfield.cc:96
 hedfield.cc:97
 hedfield.cc:98
 hedfield.cc:99
 hedfield.cc:100
 hedfield.cc:101
 hedfield.cc:102
 hedfield.cc:103
 hedfield.cc:104
 hedfield.cc:105
 hedfield.cc:106
 hedfield.cc:107
 hedfield.cc:108
 hedfield.cc:109
 hedfield.cc:110
 hedfield.cc:111
 hedfield.cc:112
 hedfield.cc:113
 hedfield.cc:114
 hedfield.cc:115
 hedfield.cc:116
 hedfield.cc:117
 hedfield.cc:118
 hedfield.cc:119
 hedfield.cc:120
 hedfield.cc:121
 hedfield.cc:122
 hedfield.cc:123
 hedfield.cc:124
 hedfield.cc:125
 hedfield.cc:126
 hedfield.cc:127
 hedfield.cc:128
 hedfield.cc:129
 hedfield.cc:130
 hedfield.cc:131
 hedfield.cc:132
 hedfield.cc:133
 hedfield.cc:134
 hedfield.cc:135
 hedfield.cc:136
 hedfield.cc:137
 hedfield.cc:138
 hedfield.cc:139
 hedfield.cc:140
 hedfield.cc:141
 hedfield.cc:142
 hedfield.cc:143
 hedfield.cc:144
 hedfield.cc:145
 hedfield.cc:146
 hedfield.cc:147
 hedfield.cc:148
 hedfield.cc:149
 hedfield.cc:150
 hedfield.cc:151
 hedfield.cc:152
 hedfield.cc:153
 hedfield.cc:154
 hedfield.cc:155
 hedfield.cc:156
 hedfield.cc:157
 hedfield.cc:158
 hedfield.cc:159
 hedfield.cc:160
 hedfield.cc:161
 hedfield.cc:162
 hedfield.cc:163
 hedfield.cc:164
 hedfield.cc:165
 hedfield.cc:166
 hedfield.cc:167
 hedfield.cc:168
 hedfield.cc:169
 hedfield.cc:170
 hedfield.cc:171
 hedfield.cc:172
 hedfield.cc:173
 hedfield.cc:174
 hedfield.cc:175
 hedfield.cc:176
 hedfield.cc:177
 hedfield.cc:178
 hedfield.cc:179
 hedfield.cc:180
 hedfield.cc:181
 hedfield.cc:182
 hedfield.cc:183
 hedfield.cc:184
 hedfield.cc:185
 hedfield.cc:186
 hedfield.cc:187
 hedfield.cc:188
 hedfield.cc:189
 hedfield.cc:190
 hedfield.cc:191
 hedfield.cc:192
 hedfield.cc:193
 hedfield.cc:194
 hedfield.cc:195
 hedfield.cc:196
 hedfield.cc:197
 hedfield.cc:198
 hedfield.cc:199
 hedfield.cc:200
 hedfield.cc:201
 hedfield.cc:202
 hedfield.cc:203
 hedfield.cc:204
 hedfield.cc:205
 hedfield.cc:206
 hedfield.cc:207
 hedfield.cc:208
 hedfield.cc:209
 hedfield.cc:210
 hedfield.cc:211
 hedfield.cc:212
 hedfield.cc:213
 hedfield.cc:214
 hedfield.cc:215
 hedfield.cc:216
 hedfield.cc:217
 hedfield.cc:218
 hedfield.cc:219
 hedfield.cc:220
 hedfield.cc:221
 hedfield.cc:222
 hedfield.cc:223
 hedfield.cc:224
 hedfield.cc:225
 hedfield.cc:226
 hedfield.cc:227
 hedfield.cc:228
 hedfield.cc:229
 hedfield.cc:230
 hedfield.cc:231
 hedfield.cc:232
 hedfield.cc:233
 hedfield.cc:234
 hedfield.cc:235
 hedfield.cc:236
 hedfield.cc:237
 hedfield.cc:238
 hedfield.cc:239
 hedfield.cc:240
 hedfield.cc:241
 hedfield.cc:242
 hedfield.cc:243
 hedfield.cc:244
 hedfield.cc:245
 hedfield.cc:246
 hedfield.cc:247
 hedfield.cc:248
 hedfield.cc:249
 hedfield.cc:250
 hedfield.cc:251
 hedfield.cc:252
 hedfield.cc:253
 hedfield.cc:254
 hedfield.cc:255
 hedfield.cc:256
 hedfield.cc:257
 hedfield.cc:258
 hedfield.cc:259
 hedfield.cc:260
 hedfield.cc:261
 hedfield.cc:262
 hedfield.cc:263
 hedfield.cc:264
 hedfield.cc:265
 hedfield.cc:266
 hedfield.cc:267
 hedfield.cc:268
 hedfield.cc:269
 hedfield.cc:270
 hedfield.cc:271
 hedfield.cc:272
 hedfield.cc:273
 hedfield.cc:274
 hedfield.cc:275
 hedfield.cc:276
 hedfield.cc:277
 hedfield.cc:278
 hedfield.cc:279
 hedfield.cc:280
 hedfield.cc:281
 hedfield.cc:282
 hedfield.cc:283
 hedfield.cc:284
 hedfield.cc:285
 hedfield.cc:286
 hedfield.cc:287
 hedfield.cc:288
 hedfield.cc:289
 hedfield.cc:290
 hedfield.cc:291
 hedfield.cc:292
 hedfield.cc:293
 hedfield.cc:294
 hedfield.cc:295
 hedfield.cc:296
 hedfield.cc:297
 hedfield.cc:298
 hedfield.cc:299
 hedfield.cc:300
 hedfield.cc:301
 hedfield.cc:302
 hedfield.cc:303
 hedfield.cc:304
 hedfield.cc:305
 hedfield.cc:306
 hedfield.cc:307
 hedfield.cc:308
 hedfield.cc:309
 hedfield.cc:310
 hedfield.cc:311
 hedfield.cc:312
 hedfield.cc:313
 hedfield.cc:314
 hedfield.cc:315
 hedfield.cc:316
 hedfield.cc:317
 hedfield.cc:318
 hedfield.cc:319
 hedfield.cc:320
 hedfield.cc:321
 hedfield.cc:322
 hedfield.cc:323
 hedfield.cc:324
 hedfield.cc:325
 hedfield.cc:326
 hedfield.cc:327
 hedfield.cc:328
 hedfield.cc:329
 hedfield.cc:330
 hedfield.cc:331
 hedfield.cc:332
 hedfield.cc:333
 hedfield.cc:334
 hedfield.cc:335
 hedfield.cc:336
 hedfield.cc:337
 hedfield.cc:338
 hedfield.cc:339
 hedfield.cc:340
 hedfield.cc:341
 hedfield.cc:342
 hedfield.cc:343
 hedfield.cc:344
 hedfield.cc:345
 hedfield.cc:346
 hedfield.cc:347
 hedfield.cc:348
 hedfield.cc:349
 hedfield.cc:350
 hedfield.cc:351
 hedfield.cc:352
 hedfield.cc:353
 hedfield.cc:354
 hedfield.cc:355
 hedfield.cc:356
 hedfield.cc:357
 hedfield.cc:358
 hedfield.cc:359
 hedfield.cc:360
 hedfield.cc:361
 hedfield.cc:362
 hedfield.cc:363
 hedfield.cc:364