Front absorber analysis routine. (Geant3 dependent !)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Jan 2001 13:18:03 +0000 (13:18 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Jan 2001 13:18:03 +0000 (13:18 +0000)
STRUCT/EvalAbso.C [new file with mode: 0644]

diff --git a/STRUCT/EvalAbso.C b/STRUCT/EvalAbso.C
new file mode 100644 (file)
index 0000000..9ac26d6
--- /dev/null
@@ -0,0 +1,228 @@
+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            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",
+                  i+1, name, a, z, 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();
+
+}
+
+
+
+