[u/mrichter/AliRoot.git] / STRUCT / EvalAbso.C

void EvalAbso()
{
// Material numbers
    enum {kC=6, kAl=9, kFe=10, kCu=11, kW=12, kPb=13,
	  kNiCuW=21, kVacuum=16, kAir=15, kConcrete=17,
	  kPolyCH2=18, kSteel=10, kInsulation=14, kPolyCc=20};	   

    AliABSO *pABSO  = (AliABSO*) gAlice->GetModule("ABSO");

    Int_t nl[2];
    nl[0] = pABSO->NumberOfLayers(0);
    nl[1] = pABSO->NumberOfLayers(1);    
    
    Int_t   mat[2][15],  mat[2][15];
    Float_t zmin[2][15], zmin[2][15], zmax[2][15], zmax[2][15];
    Int_t i, j;

    for (j=0; j<   2; j++) {
	for (i=0; i< nl[j]; i++) {
	    mat[j][i]  = pABSO->MaterialOfLayer(j,i)-1599;
	    zmax[j][i] = pABSO->ZPositionOfLayer(j,i);
	    if (i+1 < nl[j]) zmin[j][i+1] = pABSO->ZPositionOfLayer(j,i);
	}
	zmin[j][0]=0.;
    }
    Float_t l = zmax[0][nl[0]-1];


//

//
// 1. Limiting angular resolution in the Branson formalism
//
    Float_t f0,f1,f2;
    Float_t a, z, dens, radl, absl;
    char  name[21];

    for (j=0; j< 2; j++) {
	printf("\n                        A            Z           ZPos         DZ            dens        radL        absL");
	printf("\n==================================================================================");
	f0=f1=f2=0;
	for (i=0; i< nl[j]; i++) {
	    pABSO->AliGetMaterial(mat[j][i], name, a, z, dens, radl, absl);
	    Float_t dz = zmax[j][i]-zmin[j][i];
	    f0 += dz/radl;
	    f1 += dz/(2.*radl)*(2.*zmin[j][i]+dz);
	    f2 += dz/(3.*radl)*(3.*zmin[j][i]*zmin[j][i]+3.*dz*zmin[j][i]+dz*dz);	
	    printf("\n %3d %14s %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f",
		   i+1, name, a, z, (zmax[j][i]+zmin[j][i])/2., dz, dens, radl, absl);
	}
//  
	Float_t deltaThetaB = 0.0136  * TMath::Sqrt(f0-f1*f1/f2);
	Float_t lBranson    = f2/f1;
	printf("\n==================================================================================");
	printf("\n Branson plane:                              %12.3f (cm)", lBranson);
	printf("\n Limiting resolution (Branson)             : %12.3f (mrad)", deltaThetaB*1000.);
    }
//
//  2. Limiting Resolution as a function of thickness of extra Fe
//
    Float_t zmin0[15], zmax0[15];
    Float_t x[100], y[100];

    for (i=0; i< nl[0]; i++) {
	zmin0[i] = zmin[0][i];
	zmax0[i] = zmax[0][i];	
    }
    
    Float_t ds = (zmax[0][3]-zmin[0][2])/100.;
    Int_t j;
    
    for (j=0; j< 100; j++) {
	f0=f1=f2=0;
	Float_t zFe = zmin[0][2]+Float_t(j)*ds;
	zmax0[2]=zFe;
	zmin0[3]=zFe;
	for (i=0; i< nl[0]; i++) {
	    Float_t a, z, dens, radl, absl;
	    pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
	    Float_t dz = zmax0[i]-zmin0[i];
	    f0 += dz/radl;
	    f1 += dz/(2.*radl)*(2.*zmin0[i]+dz);
	    f2 += dz/(3.*radl)*(3.*zmin0[i]*zmin0[i]+3.*dz*zmin0[i]+dz*dz);	
	}
	Float_t deltaThetaB = 0.0136  * TMath::Sqrt(f0-f1*f1/f2);
	x[j] = zmax[0][3]-zFe;
	y[j] = deltaThetaB*1000.;
    }
    TGraph *sGraph = new TGraph(100, x, y);
    TCanvas *c = new TCanvas("c","  ",400,10,600,700);
    c->Divide(2,2);
    c->cd(1);
    sGraph->Draw("AC");
    c->Update();
//
// 3. Limiting Resolution as a function of density of Carbon
//
    Float_t dRho = (2.5-1.5)/100.;

    for (j=0; j< 100; j++) {
	f0=f1=f2=0;
	Float_t rho = 1.5 + Float_t(j)*dRho;
	for (i=0; i< nl[0]; i++) {
	    Float_t a, z, dens, radl, absl;
	    pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
	    if (i==1) radl*=1.75/rho;
	    Float_t dz = zmax[0][i]-zmin[0][i];
	    f0 += dz/radl;
	    f1 += dz/(2.*radl)*(2.*zmin[0][i]+dz);
	    f2 += dz/(3.*radl)*(3.*zmin[0][i]*zmin[0][i]+3.*dz*zmin[0][i]+dz*dz);	
	}
	deltaThetaB = 0.0136  * TMath::Sqrt(f0-f1*f1/f2);
	x[j] = rho;
	y[j] = deltaThetaB*1000.;
    }
    TGraph *dGraph = new TGraph(100, x, y);
    c->cd(2);
    dGraph->Draw("AC");
    c->Update();

//
//   4. Energy Loss
//
//
    TCanvas *c1 = new TCanvas("c1","  ",400,10,600,700);
    c1->Divide(2,2);
 
    const Float_t kAvo=0.60221367;

    char* tmp;
    tmp = new char[5];
    strncpy(tmp, "LOSS", 4);
    tmp[4]='\0';

    Int_t ixst;
    Float_t ekin[100], dedx[100], de[100], pcut[100], deRad[100];
    Float_t emin =   2;    
    Float_t emax = 200;
    Float_t eps = (emax-emin)/100.;
//  3 < theta < 9
    for (j=0; j< 100; j++) {
	ekin[j] = emin + Float_t(j)*eps;
	de[j]=0;
	deRad[j]=0;
    }
// all losses    
    for (i=0; i< nl[0]; i++) {
	((TGeant3*) gMC)->Gftmat(pABSO->GetMatId(mat[0][i]), 5, tmp, 100, ekin, dedx, pcut, ixst);
	for (j=0; j< 100; j++) de[j]+=dedx[j]*(zmax[0][i]-zmin[0][i])/1000.;
    }
// radative losses
    for (i=0; i< nl[0]; i++) {
	Float_t a, z, dens, radl, absl;
	pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
	for (Int_t j=0; j<100; j++) {
	    Float_t nor= kAvo*dens/a*(zmax[0][i]-zmin[0][i]);
	    deRad[j] += (((TGeant3*) gMC)->Gbrelm(z,ekin[j],1.e10))*nor;
	    deRad[j] += (((TGeant3*) gMC)->Gprelm(z,ekin[j],1.e10))*nor;
	}
    }
    
    TH2F *hr = new TH2F("hr1","Several graphs in the same pad",2,0,15,2,0,5);
    hr->SetXTitle("X title");
    hr->SetYTitle("Y title");
    c1->cd(1);
    hr->Draw();

    TGraph *eGraph1  = new TGraph(100, ekin, de);
    TGraph *eGraphR1 = new TGraph(100, ekin, deRad);

    eGraph1 ->Draw("LP");
    eGraphR1->Draw("LP");

    hr = new TH2F("hr2","Several graphs in the same pad",2,15,200,2,0,5);
    c1->cd(2);
    hr->Draw();
    eGraph1 ->Draw("LP");
    eGraphR1->Draw("LP");
    

//  2 < theta < 3
    for (j=0; j< 100; j++) {
	ekin[j] = emin + Float_t(j)*eps;
	de[j]=0;
	deRad[j]=0;
    }
// all losses    
    for (i=0; i< nl[1]; i++) {
	((TGeant3*) gMC)->Gftmat(pABSO->GetMatId(mat[1][i]), 5, tmp, 100, ekin, dedx, pcut, ixst);
	for (j=0; j< 100; j++) de[j]+=dedx[j]*(zmax[1][i]-zmin[1][i])/1000.;
    }
// radative losses
    for (i=0; i< nl[1]; i++) {
	Float_t a, z, dens, radl, absl;
	pABSO->AliGetMaterial(mat[1][i], name, a, z, dens, radl, absl);
	for (Int_t j=0; j<100; j++) {
	    Float_t nor= kAvo*dens/a*(zmax[1][i]-zmin[1][i]);
	    deRad[j] += (((TGeant3*) gMC)->Gbrelm(z,ekin[j],1.e10))*nor;
	    deRad[j] += (((TGeant3*) gMC)->Gprelm(z,ekin[j],1.e10))*nor;
	}
    }

    TH2F *hr2 = new TH2F("hr3","Several graphs in the same pad",2,0,15,2,0,5);
    hr2->SetXTitle("X title");
    hr2->SetYTitle("Y title");
    c1->cd(3);
    hr2->Draw();

    TGraph *eGraph2  = new TGraph(100, ekin, de);
    TGraph *eGraphR2 = new TGraph(100, ekin, deRad);

    eGraph2 ->Draw("LP");
    eGraphR2->Draw("LP");

    hr2 = new TH2F("hr4","Several graphs in the same pad",2,15,200,2,0,5);
    c1->cd(4);
    hr2->Draw();
    eGraph2 ->Draw("LP");
    eGraphR2->Draw("LP");


    c->Update();

}
Commit	Line	Data
fd2bbc37	1	void EvalAbso()
	2	{
	3	// Material numbers
	4	enum {kC=6, kAl=9, kFe=10, kCu=11, kW=12, kPb=13,
	5	kNiCuW=21, kVacuum=16, kAir=15, kConcrete=17,
	6	kPolyCH2=18, kSteel=10, kInsulation=14, kPolyCc=20};
	7
	8	AliABSO pABSO = (AliABSO) gAlice->GetModule("ABSO");
	9
	10	Int_t nl[2];
	11	nl[0] = pABSO->NumberOfLayers(0);
	12	nl[1] = pABSO->NumberOfLayers(1);
	13
	14	Int_t mat[2][15], mat[2][15];
	15	Float_t zmin[2][15], zmin[2][15], zmax[2][15], zmax[2][15];
	16	Int_t i, j;
	17
	18	for (j=0; j< 2; j++) {
	19	for (i=0; i< nl[j]; i++) {
	20	mat[j][i] = pABSO->MaterialOfLayer(j,i)-1599;
	21	zmax[j][i] = pABSO->ZPositionOfLayer(j,i);
	22	if (i+1 < nl[j]) zmin[j][i+1] = pABSO->ZPositionOfLayer(j,i);
	23	}
	24	zmin[j][0]=0.;
	25	}
	26	Float_t l = zmax[0][nl[0]-1];
	27
	28
	29	//
	30
	31	//
	32	// 1. Limiting angular resolution in the Branson formalism
	33	//
	34	Float_t f0,f1,f2;
	35	Float_t a, z, dens, radl, absl;
	36	char name[21];
	37
	38	for (j=0; j< 2; j++) {
92664bc8	39	printf("\n A Z ZPos DZ dens radL absL");
fd2bbc37	40	printf("\n==================================================================================");
	41	f0=f1=f2=0;
	42	for (i=0; i< nl[j]; i++) {
	43	pABSO->AliGetMaterial(mat[j][i], name, a, z, dens, radl, absl);
	44	Float_t dz = zmax[j][i]-zmin[j][i];
	45	f0 += dz/radl;
	46	f1 += dz/(2.radl)(2.*zmin[j][i]+dz);
	47	f2 += dz/(3.radl)(3.zmin[j][i]zmin[j][i]+3.dzzmin[j][i]+dz*dz);
92664bc8	48	printf("\n %3d %14s %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f",
92664bc8	49	i+1, name, a, z, (zmax[j][i]+zmin[j][i])/2., dz, dens, radl, absl);
fd2bbc37	50	}
	51	//
	52	Float_t deltaThetaB = 0.0136 * TMath::Sqrt(f0-f1*f1/f2);
	53	Float_t lBranson = f2/f1;
	54	printf("\n==================================================================================");
	55	printf("\n Branson plane: %12.3f (cm)", lBranson);
	56	printf("\n Limiting resolution (Branson) : %12.3f (mrad)", deltaThetaB*1000.);
	57	}
	58	//
	59	// 2. Limiting Resolution as a function of thickness of extra Fe
	60	//
	61	Float_t zmin0[15], zmax0[15];
	62	Float_t x[100], y[100];
	63
	64	for (i=0; i< nl[0]; i++) {
	65	zmin0[i] = zmin[0][i];
	66	zmax0[i] = zmax[0][i];
	67	}
	68
	69	Float_t ds = (zmax[0][3]-zmin[0][2])/100.;
	70	Int_t j;
	71
	72	for (j=0; j< 100; j++) {
	73	f0=f1=f2=0;
	74	Float_t zFe = zmin[0][2]+Float_t(j)*ds;
	75	zmax0[2]=zFe;
	76	zmin0[3]=zFe;
	77	for (i=0; i< nl[0]; i++) {
	78	Float_t a, z, dens, radl, absl;
	79	pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
	80	Float_t dz = zmax0[i]-zmin0[i];
	81	f0 += dz/radl;
	82	f1 += dz/(2.radl)(2.*zmin0[i]+dz);
	83	f2 += dz/(3.radl)(3.zmin0[i]zmin0[i]+3.dzzmin0[i]+dz*dz);
	84	}
	85	Float_t deltaThetaB = 0.0136 * TMath::Sqrt(f0-f1*f1/f2);
	86	x[j] = zmax[0][3]-zFe;
	87	y[j] = deltaThetaB*1000.;
	88	}
	89	TGraph *sGraph = new TGraph(100, x, y);
	90	TCanvas *c = new TCanvas("c"," ",400,10,600,700);
	91	c->Divide(2,2);
	92	c->cd(1);
	93	sGraph->Draw("AC");
	94	c->Update();
	95	//
	96	// 3. Limiting Resolution as a function of density of Carbon
	97	//
	98	Float_t dRho = (2.5-1.5)/100.;
	99
	100	for (j=0; j< 100; j++) {
	101	f0=f1=f2=0;
	102	Float_t rho = 1.5 + Float_t(j)*dRho;
	103	for (i=0; i< nl[0]; i++) {
	104	Float_t a, z, dens, radl, absl;
	105	pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
	106	if (i==1) radl*=1.75/rho;
	107	Float_t dz = zmax[0][i]-zmin[0][i];
	108	f0 += dz/radl;
	109	f1 += dz/(2.radl)(2.*zmin[0][i]+dz);
	110	f2 += dz/(3.radl)(3.zmin[0][i]zmin[0][i]+3.dzzmin[0][i]+dz*dz);
	111	}
	112	deltaThetaB = 0.0136 * TMath::Sqrt(f0-f1*f1/f2);
	113	x[j] = rho;
114	y[j] = deltaThetaB*1000.;
115	}
116	TGraph *dGraph = new TGraph(100, x, y);
117	c->cd(2);
118	dGraph->Draw("AC");
119	c->Update();
120
121	//
122	// 4. Energy Loss
123	//
124	//
125	TCanvas *c1 = new TCanvas("c1"," ",400,10,600,700);
126	c1->Divide(2,2);
127
128	const Float_t kAvo=0.60221367;
129
130	char* tmp;
131	tmp = new char[5];
132	strncpy(tmp, "LOSS", 4);
133	tmp[4]='\0';
134
135	Int_t ixst;
136	Float_t ekin[100], dedx[100], de[100], pcut[100], deRad[100];
137	Float_t emin = 2;
138	Float_t emax = 200;
139	Float_t eps = (emax-emin)/100.;
140	// 3 < theta < 9
141	for (j=0; j< 100; j++) {
142	ekin[j] = emin + Float_t(j)*eps;
143	de[j]=0;
144	deRad[j]=0;
145	}
146	// all losses
147	for (i=0; i< nl[0]; i++) {
148	((TGeant3*) gMC)->Gftmat(pABSO->GetMatId(mat[0][i]), 5, tmp, 100, ekin, dedx, pcut, ixst);
149	for (j=0; j< 100; j++) de[j]+=dedx[j]*(zmax[0][i]-zmin[0][i])/1000.;
150	}
151	// radative losses
152	for (i=0; i< nl[0]; i++) {
153	Float_t a, z, dens, radl, absl;
154	pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
155	for (Int_t j=0; j<100; j++) {
156	Float_t nor= kAvodens/a(zmax[0][i]-zmin[0][i]);
157	deRad[j] += (((TGeant3) gMC)->Gbrelm(z,ekin[j],1.e10))nor;
158	deRad[j] += (((TGeant3) gMC)->Gprelm(z,ekin[j],1.e10))nor;
159	}
160	}
161
162	TH2F *hr = new TH2F("hr1","Several graphs in the same pad",2,0,15,2,0,5);
163	hr->SetXTitle("X title");
164	hr->SetYTitle("Y title");
165	c1->cd(1);
166	hr->Draw();
167
168	TGraph *eGraph1 = new TGraph(100, ekin, de);
169	TGraph *eGraphR1 = new TGraph(100, ekin, deRad);
170
171	eGraph1 ->Draw("LP");
172	eGraphR1->Draw("LP");
173
174	hr = new TH2F("hr2","Several graphs in the same pad",2,15,200,2,0,5);
175	c1->cd(2);
176	hr->Draw();
177	eGraph1 ->Draw("LP");
178	eGraphR1->Draw("LP");
179
180
181	// 2 < theta < 3
182	for (j=0; j< 100; j++) {
183	ekin[j] = emin + Float_t(j)*eps;
184	de[j]=0;
185	deRad[j]=0;
186	}
187	// all losses
188	for (i=0; i< nl[1]; i++) {
189	((TGeant3*) gMC)->Gftmat(pABSO->GetMatId(mat[1][i]), 5, tmp, 100, ekin, dedx, pcut, ixst);
190	for (j=0; j< 100; j++) de[j]+=dedx[j]*(zmax[1][i]-zmin[1][i])/1000.;
191	}
192	// radative losses
193	for (i=0; i< nl[1]; i++) {
194	Float_t a, z, dens, radl, absl;
195	pABSO->AliGetMaterial(mat[1][i], name, a, z, dens, radl, absl);
196	for (Int_t j=0; j<100; j++) {
197	Float_t nor= kAvodens/a(zmax[1][i]-zmin[1][i]);
198	deRad[j] += (((TGeant3) gMC)->Gbrelm(z,ekin[j],1.e10))nor;
199	deRad[j] += (((TGeant3) gMC)->Gprelm(z,ekin[j],1.e10))nor;
200	}
201	}
202
203	TH2F *hr2 = new TH2F("hr3","Several graphs in the same pad",2,0,15,2,0,5);
204	hr2->SetXTitle("X title");
205	hr2->SetYTitle("Y title");
206	c1->cd(3);
207	hr2->Draw();
208
209	TGraph *eGraph2 = new TGraph(100, ekin, de);
210	TGraph *eGraphR2 = new TGraph(100, ekin, deRad);
211
212	eGraph2 ->Draw("LP");
213	eGraphR2->Draw("LP");
214
215	hr2 = new TH2F("hr4","Several graphs in the same pad",2,15,200,2,0,5);
216	c1->cd(4);
217	hr2->Draw();
218	eGraph2 ->Draw("LP");
219	eGraphR2->Draw("LP");
220
221
222	c->Update();
223
224	}
225
226
227
228