STRUCT/EvalAbso.C

   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++) {
  39         printf("\n                        A            Z            dens        radL        absL");
  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.*dz*zmin[j][i]+dz*dz);
  48             printf("\n %3d %14s %12.3f %12.3f %12.3f %12.3f %12.3f",
  49                    i+1, name, a, z, dens, radl, absl);
  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.*dz*zmin0[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.*dz*zmin[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= kAvo*dens/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= kAvo*dens/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