1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 Revision 1.2 2000/07/12 08:56:25 fca
19 Coding convention correction and warning removal
21 Revision 1.1 2000/07/11 18:24:59 fca
22 Coding convention corrections + few minor bug fixes
28 #include "AliMagFDM.h"
34 //________________________________________
35 AliMagFDM::AliMagFDM(const char *name, const char *title, const Int_t integ,
36 const Int_t map, const Float_t factor, const Float_t fmax)
37 : AliMagF(name,title,integ,map,factor,fmax)
41 // Standard constructor for the Dipole field
45 printf("Field Map for Muon Arm from IP till muon filter %s created: map= %d, factor= %f, file=%s\n",fName.Data(),map,factor,fTitle.Data());
49 //________________________________________
51 void AliMagFDM::Field(Float_t *xfi, Float_t *b)
54 // Main routine to compute the field in a point
56 static const Double_t keps=0.1E-06;
57 static const Double_t kpi2=.6283185E+01;
58 static const Double_t kone=1;
60 static const Int_t kiip=33;
61 static const Int_t kmiip=0;
62 static const Int_t kliip=0;
64 static const Int_t kiic=0;
65 static const Int_t kmiic=0;
66 static const Int_t kliic=0;
68 static const Double_t kfdZbg=502.92; // Start of Map using in z
69 static const Double_t kfdZL3=600; // Beginning of L3 door in z
80 Double_t zp1, zp2,xp1,xp2,yp1,yp2;
81 Double_t zz1, zz2,yy1,yy2,x2,x1;
83 // --- start the map fiel from z = 502.92 cm ---
89 // printf("x[0] %f,x[1] %f,x[2] %f\n",x[0],x[1],x[2]);
91 Double_t rr=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
95 if ( (-700<x[2] && x[2]<=kfdZbg &&
96 (x[0]*x[0]+(x[1]+30)*(x[1]+30))< 560*560)
97 || (kfdZbg<x[2] && x[2]<=kfdZL3 && rr>=rPmax*100) )
103 xL3[1]=(x[1]+30)/100;
106 Double_t xminn=xL3[2]*fdAx1+fdCx1;
107 Double_t xmaxx=xL3[2]*fdAx2+fdCx2;
108 Double_t zCmin,zCmax,yCmin,yCmax;
115 if ((kfdZbg/100<xL3[2] && xL3[2]<zCmin && r0<rPmax) || ((zCmin<=xL3[2] && xL3[2] <= zCmax ) && (yCmin<=xL3[1] && xL3[1]<= yCmax) && (xminn <= xL3[0] && xL3[0] <= xmaxx)))
119 if (xL3[2]<zCmin && r0<rPmax)
121 //--------------------- Polar part ----------------------
126 Double_t alp1,alp2,alp3;
127 Double_t zpz,rp,fip,cphi;
135 FZ(&zpz, fdZp ,&fdZpdl,&kpi,&kp0,&zp1 ,&zp2,&fdZpl) ;
139 ph0=TMath::ACos(cphi);
140 if (xL3[0]< 0) {ph0=kpi2 - ph0;}
143 FZ(&fip,fdPhi,&fdPhid ,&mpi,&mp0, &xp1,&xp2,&fdPhin);
154 bint[0]=(zp1*fdBpx[kp0][0][0] + zp2*fdBpx[kp0+1][0][0])*10;
155 bint[1]=(zp1*fdBpy[kp0][0][0] + zp2*fdBpy[kp0+1][0][0])*10;
156 bint[2]=(zp1*fdBpz[kp0][0][0] + zp2*fdBpz[kp0+1][0][0])*10;
160 alp2= fdB[0][0][mp0]*yyp + fdB[0][1][mp0]*xL3[0];
161 alp3= fdB[1][0][mp0]*yyp + fdB[1][1][mp0]*xL3[0];
162 alp1= kone - alp2 - alp3;
164 for (jb=0; jb<3 ; jb++)
167 FRfuncBi(&iKvar,&zp1,&zp2,&alp1,&alp2,&alp3, &kp0,&mp0, &bbj);
175 FZ(&rp,fdR ,&fdRdel,&lpi,&lp0,&yp1,&yp2,&fdRn);
177 for (jb=0; jb<3 ; jb++)
180 FGfuncBi(&zp1,&zp2,&yp1,&yp2,&xp1,&xp2,&iKvar,&kp0,&lp0,&mp0,&bbj);
190 // fprintf(fitest,"------------- Freg2 run -------------\n");
195 //-------------- Cartensian part ------------------
199 Double_t xx1, xx2,dx, xxx ,xXl;
206 xx1 = fdAx1*xL3[2] + fdCx1;
207 xx2 = fdAx2*xL3[2] + fdCx2;
210 FZ(&zzc, fdZc ,&fdZdel, &kci,&k0, &zz1, &zz2, &fdZl);
213 FZ(&yyc, fdY , &fdYdel,&lci, &l0, &yy1, &yy2,&fdYl);
221 if(xL3[0]<(xx1+m0*dx) || xL3[0] >(xx1+(m0+1)*dx))
224 printf(" m0 %d, m0+1 %d\n",m0,m0+1);
227 x2=(xL3[0]-( xx1+m0*dx))/dx;
230 for (jb=3; jb<6; jb++)
233 FGfuncBi(&zz1,&zz2,&yy1,&yy2,&x1,&x2,&iKvar,&k0, &l0, &m0, &bbj);
234 bint[jb-3] = bbj*10 ;
241 // fprintf(fitest,"------------ Freg1 run -----------------\n");
245 printf("Unknown map of Dipole region %d\n",fMap);
250 //This is the ZDC part
251 Float_t rad2=x[0]*x[0]+x[1]*x[1];
255 // Separator Dipole D2
256 if(x[2]<kD2END) b[1]=kFDIP;
257 } else if(x[2]>kD1BEG) {
259 // Separator Dipole D1
260 if(x[2]<kD1END) b[1]=-kFDIP;
264 // First quadrupole of inner triplet de-focussing in x-direction
273 } else if(x[2]>kZ3BEG) {
280 } else if(x[2]>kZ2BEG) {
287 } else if(x[2]>kZ1BEG) {
294 } else if(x[2]>kCORBEG) {
296 // Corrector dipole (because of dimuon arm)
312 //_________________________________________
314 void AliMagFDM::FZ(Double_t *u, Float_t *Ar, Float_t *du,Int_t *ki,Int_t *kf,Double_t *a1,Double_t *a2 ,Int_t *nu)
318 // Z component of the field
320 static const Double_t kone=1;
323 Double_t ddu,delu,ar;
338 for(l=ikj; l<nk; l++)
344 *a2=(temp-Ar[l])/(Ar[l+1]-Ar[l]);
351 /*-------------FRfuncBi----------------*/
353 void AliMagFDM::FRfuncBi(Int_t *kai,Double_t *za1, Double_t *za2, Double_t *al1, Double_t *al2, Double_t *al3, Int_t *ka, Int_t *ma, Double_t *ba)
356 // This method needs to be commented
358 Double_t fa11,fa12,fa13;
359 Double_t fa21,fa22,fa23;
363 Double_t zaa1,zaa2,alf1,alf2,alf3;
376 fa11 = fdBpx[kaa][0][0];
377 fa12 = fdBpx[kaa][0][maa];
378 fa13 = fdBpx[kaa][0][maa+1];
379 fa21 = fdBpx[kaa+1][0][0];
380 fa22 = fdBpx[kaa+1][0][maa];
381 fa23 = fdBpx[kaa+1][0][maa+1];
384 fa11 = fdBpy[kaa][0][0];
385 fa12 = fdBpy[kaa][0][maa];
386 fa13 = fdBpy[kaa][0][maa+1];
387 fa21 = fdBpy[kaa+1][0][0];
388 fa22 = fdBpy[kaa+1][0][maa];
389 fa23 = fdBpy[kaa+1][0][maa+1];
392 fa11 = fdBpz[kaa][0][0];
393 fa12 = fdBpz[kaa][0][maa];
394 fa13 = fdBpz[kaa][0][maa+1];
395 fa21 = fdBpz[kaa+1][0][0];
396 fa22 = fdBpz[kaa+1][0][maa];
397 fa23 = fdBpz[kaa+1][0][maa+1];
400 Fatal("FRfuncBi","Invalid value of kaai %d\n",kaai);
403 faY1=alf1*fa11+alf2*fa12+alf3*fa13;
404 faY2=alf1*fa21+alf2*fa22+alf3*fa23;
405 bba = zaa1*faY1+zaa2*faY2;
411 /*----------- FGfuncBi------------*/
413 void AliMagFDM::FGfuncBi(Double_t *zz1,Double_t *zz2, Double_t *yy1,Double_t *yy2, Double_t *xx1,Double_t *xx2, Int_t *kvr, Int_t *kk, Int_t *ll, Int_t *mm, Double_t *bb)
416 // This method needs to be commented
418 Double_t fy1, fy2, ffy;
419 Double_t gy1,gy2,ggy;
420 Double_t z1,z2,y1,y2,x1,x2;
425 Double_t bf11,bf12,bf21,bf22;
426 Double_t bg11,bg12,bg21,bg22;
440 /*-----------------Polar part ------------------*/
445 bf12=fdBpx[k+1][l][m];
446 bf21=fdBpx[k+1][l+1][m];
447 bf22=fdBpx[k][l+1][m];
449 bg11=fdBpx[k][l][m+1];
450 bg12=fdBpx[k+1][l][m+1];
451 bg21=fdBpx[k+1][l+1][m+1];
452 bg22=fdBpx[k][l+1][m+1];
457 bf12=fdBpy[k+1][l][m];
458 bf21=fdBpy[k+1][l+1][m];
459 bf22=fdBpy[k][l+1][m];
461 bg11=fdBpy[k][l][m+1];
462 bg12=fdBpy[k+1][l][m+1];
463 bg21=fdBpy[k+1][l+1][m+1];
464 bg22=fdBpy[k][l+1][m+1];
469 bf12=fdBpz[k+1][l][m];
470 bf21=fdBpz[k+1][l+1][m];
471 bf22=fdBpz[k][l+1][m];
473 bg11=fdBpz[k][l][m+1];
474 bg12=fdBpz[k+1][l][m+1];
475 bg21=fdBpz[k+1][l+1][m+1];
476 bg22=fdBpz[k][l+1][m+1];
478 /*-----------------Cartensian part ---------------*/
482 bf12=fdBcx[k+1][l][m];
483 bf21=fdBcx[k+1][l+1][m];
484 bf22=fdBcx[k][l+1][m];
486 bg11=fdBcx[k][l][m+1];
487 bg12=fdBcx[k+1][l][m+1];
488 bg21=fdBcx[k+1][l+1][m+1];
489 bg22=fdBcx[k][l+1][m+1];
494 bf12=fdBcy[k+1][l][m];
495 bf21=fdBcy[k+1][l+1][m];
496 bf22=fdBcy[k][l+1][m];
498 bg11=fdBcy[k][l][m+1];
499 bg12=fdBcy[k+1][l][m+1];
500 bg21=fdBcy[k+1][l+1][m+1];
501 bg22=fdBcy[k][l+1][m+1];
506 bf12=fdBcz[k+1][l][m];
507 bf21=fdBcz[k+1][l+1][m];
508 bf22=fdBcz[k][l+1][m];
510 bg11=fdBcz[k][l][m+1];
511 bg12=fdBcz[k+1][l][m+1];
512 bg21=fdBcz[k+1][l+1][m+1];
513 bg22=fdBcz[k][l+1][m+1];
517 Fatal("FGfuncBi","Invalid value of kv %d\n",kv);
524 fy2=z2*bf21+z1* bf22;
528 gy1 = z1*bg11+z2*bg12;
529 gy2 = z2*bg21+z1*bg22;
530 ggy= y1*gy1 + y2*gy2;
537 //____________________________________________
539 void AliMagFDM::ReadField()
542 // Method to read the magnetic field from file
547 Float_t zzp, rr,phii;
548 Float_t zz, yy, bx,by,bz,bb;
551 printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
552 fname = gSystem->ExpandPathName(fTitle.Data());
553 magfile=fopen(fname,"r");
556 printf("Cartensian part\n");
562 fscanf(magfile,"%d %d %d ",&fdYl, &fdXl, &fdZl);
564 printf("fdYl %d, fdXl %d, fdZl %d\n",fdYl, fdXl, fdZl);
566 for (ik=0; ik<fdZl; ik++)
569 fscanf(magfile, " %e ", &zz);
574 for (ik=0; ik<fdYl; ik++)
576 fscanf(magfile, " %e ", &yy);
580 for (ik=0; ik<81; ik++)
582 printf("fdZc %e,fdY %e\n", fdZc[ik],fdY[ik]);
585 fscanf(magfile," %e %e %e %e %e %e %e %e %e %e %e ", &fdYdel,&fdXdel,&fdZdel,&fdZmax,&fdZmin,&fdYmax,&fdYmin,&fdAx1,&fdCx1,&fdAx2,&fdCx2);
587 printf("fdYdel %e, fdXdel %e, fdZdel %e\n",fdYdel,fdXdel,fdZdel);
588 printf("fdZmax %e, fdZmin %e, fdYmax %e,fdYmin %e\n",fdZmax,fdZmin,fdYmax,fdYmin);
589 printf("fdAx1 %e, fdCx1 %e, fdAx2 %e, fdCx %e\n",fdAx1,fdCx1,fdAx2,fdCx2);
591 for (il=0; il<44; il++) {
592 for (im=0; im<81; im++) {
593 for (ik=0; ik<81; ik++) {
595 fscanf(magfile, " %e ", &by);
596 fdBcy[ik][im][il]=by;
601 for (il=0; il<44; il++) {
602 for (im=0; im<81; im++) {
603 for (ik=0; ik<81; ik++) {
605 fscanf(magfile, " %e ", &bx);
606 fdBcx[ik][im][il]=bx;
611 for (il=0; il<44; il++) {
612 for (im=0; im<81; im++) {
613 for (ik=0; ik<81; ik++) {
615 fscanf(magfile, " %e ", &bz);
616 fdBcz[ik][im][il]=bz;
620 //---------------------- Polar part ---------------------------------
622 printf("Polar part\n");
623 fscanf(magfile,"%d %d %d ", &fdZpl, &fdRn, &fdPhin);
624 printf("fdZpl %d, fdRn %d, fdPhin %d\n",fdZpl,fdRn,fdPhin);
626 printf(" fdZp array\n");
628 for (ik=0; ik<51; ik++)
630 fscanf(magfile, " %e ", &zzp);
632 printf(" %e\n",fdZp[ik]);
635 printf(" fdR array\n");
637 for (ik=0; ik<10; ik++)
639 fscanf(magfile, " %e ", &rr);
641 printf(" %e\n",fdR[ik]);
644 // printf("fdPhi array\n");
646 for (il=0; il<33; il++)
648 fscanf(magfile, " %e ", &phii);
650 // printf(" %e\n",fdPhi[il]);
653 fscanf(magfile," %e %e %e %e %e %e %e ",&fdZpdl,&fdPhid,&fdRdel,&fdZpmx,&fdZpmn,&fdRmax, &fdRmin);
655 printf("fdZpdl %e, fdPhid %e, fdRdel %e, fdZpmx %e, fdZpmn %e,fdRmax %e,fdRmin %e \n", fdZpdl,fdPhid, fdRdel,fdZpmx, fdZpmn,fdRmax, fdRmin);
658 for (il=0; il<33; il++) {
659 for (im=0; im<10; im++) {
660 for (ik=0; ik<51; ik++) {
661 fscanf(magfile, " %e ", &by);
662 fdBpy[ik][im][il]=by;
667 for (il=0; il<33; il++) {
668 for (im=0; im<10; im++) {
669 for (ik=0; ik<51; ik++) {
670 fscanf(magfile, " %e ", &bx);
671 fdBpx[ik][im][il]=bx;
677 for (il=0; il<33; il++) {
678 for (im=0; im<10; im++) {
679 for (ik=0; ik<51; ik++) {
680 fscanf(magfile, " %e ", &bz);
681 fdBpz[ik][im][il]=bz;
687 for (il=0; il<32; il++) {
688 for (im=0; im<2; im++) {
689 for (ik=0; ik<2; ik++) {
690 fscanf(magfile, " %e ", &bb);
697 printf("File %s not found !\n",fTitle.Data());
701 //________________________________