{
  TCanvas *fMainCanvas;
  fMainCanvas = new TCanvas("MainCanvas","Main Canvas of Litrani",100,5,740,600);
  fMainCanvas->SetFillColor(30);
  fMainCanvas->SetGrid();
  const Float_t xlowPad1 = 0.02;   //xlow of fMainPad1 in fMainCanvas
  const Float_t ylowPad1 = 0.02;   //ylow of fMainPad1 in fMainCanvas
  const Float_t xupPad1  = 0.98;   //xup of fMainPad1 in fMainCanvas
  const Float_t yupPad1  = 0.98;   //yup of fMainPad1 in fMainCanvas
  const Color_t colPad1  = 31;     //color of fMainPad1
  const Short_t bszPad   = 3;      //border size of both pads
  fMainPad1 = new TPad("MainPad1","Pad1 in MainCanvas",xlowPad1,ylowPad1,xupPad1,yupPad1,
        colPad1,bszPad);
  fMainPad1->SetGridx();
  fMainPad1->SetGridy();
  fMainPad1->Draw();
  fMainPad1->cd();
  char c;
  cout << endl;
  TMaterial *pbwo4;
  pbwo4 = new TMaterial("PbWO4","PbWO4 crystal",1.0,1.0,1.0);
  TMaterial *glue;
  glue = new TMaterial("Glue","Glue",1.0,1.0,1.0);
  TMaterial *Silicium;
  Silicium = new TMaterial("Silicium","Silicium of APD",1.0,1.0,1.0);
  TMaterial *grease;
  grease = new TMaterial("Grease","Grease",1.0,1.0,1.0);
  TMaterial *plexi;
  plexi = new TMaterial("Plexi","Plexi",1.0,1.0,1.0);
//
//    GEOMETRY
//
  const Double_t crystal_dx1 = 1.5;    //half x length at -dz in cm
  const Double_t crystal_dx2 = 1.0;    //half x length at +dz in cm
  const Double_t crystal_dy1 = 2.5;    //half y length at -dz in cm
  const Double_t crystal_dy2 = 2.0;    //half y length at +dz in cm
  const Double_t crystal_dz  = 6.0;    //half-length along the z-axis in cm
  TTRD2 *crystal;
  crystal = new TTRD2("crystal","CMS crystal","PbWO4",crystal_dx1,
                       crystal_dx2,crystal_dy1,crystal_dy2,crystal_dz);
  crystal->SetLineColor(28);
  TNode *node1 = new TNode("node1","node1",crystal);
  node1->cd();
  const Double_t zero = 0.0;
  const Double_t z05 = 0.5;
  const Double_t pideg = 180.0;
  const Double_t Dmin = 1.0;         //false width of the bellows, if possible.
  const Double_t deltaDstart = -20.0;//starting value for delta in degree
  const Double_t deltaDstep = 5.0;   //step for delta in degree
  const Double_t bow_dz = 0.8;       //false half-length of the TSTRAP along its z axis in cm
  const Double_t bow_h1 = 0.8;       //false half-length along y of the face at -DZ of the TSTRAP
  const Double_t fib_r = 0.6;        //false radius of the fibre
  const Double_t fib_dz = 1.0;       //length of the fibre (irrelevant)
  Double_t fib_x;                    //x of the center of the fibre in WCS
  Double_t fib_y;                    //y of the center of the fibre in WCS
  Double_t fib_z;                    //z of the center of the fibre in WCS
  Double_t bow_alf1;                 //alpha1 = alpha2 of the TSTRAP
  Double_t deltastart;               //starting value for delta in radian
  Double_t deltastep;                //step for delta in radian
  Double_t torad = TMath::Pi()/pideg;
  Double_t D;                        //width of the bellows (=Dmin if possible)
  Double_t bow_center;               //shift in z to put bowel in place
  Double_t deltaD;                   //angle delta in degree
  Double_t delta;                    //angle delta in radian
  Double_t tgd;                      //tangent of angle delta
  Double_t htgd;                     //bow_h1*tgd
  Double_t bl1;                      //half-length along x of the side at -H1 in y of the
                                     // face at -DZ in z (of the TSTRAP)
  Double_t tl1;                      //half-length along x of the side at +H1 in y of the
                                     // face at -DZ in z (of the TSTRAP)
  deltastart = torad*deltaDstart;
  deltastep  = torad*deltaDstep;
  delta      = deltastart;
  deltaD     = deltaDstart;
  TTRAP *bowels;
  TRotMatrix *matrixbow;
  matrixbow = new TRotMatrix("MatrixBow", "MatrixBow",180,0,90,90,90,0);
  TNode *node4;
  TTUBE *fibre;
  TRotMatrix *matrixfib;
  TNode *node5;
  Double_t yy;
  for (Int_t run=96;run<105;run++) {
    cout << endl;
    cout << "begin loop    " << "run = " << run << endl;
    cout << endl;
    tgd  = TMath::Tan(delta);
    htgd = bow_h1*tgd;
    D    = TMath::Max(Dmin,TMath::Abs(htgd));
    bl1  = z05*(D + htgd);
    tl1  = z05*(D - htgd);
    bow_center = crystal_dz + z05*D;
    bow_alf1 = TMath::ATan((z05*(bl1-tl1))/bow_h1)/torad;
    bowels = new TTRAP("bowels","bowels","Grease",bow_dz,zero,zero,
                         bow_h1,bl1,tl1,bow_alf1,bow_h1,bl1,tl1,bow_alf1);
    bowels->SetLineColor(28);
    node4 = new TNode("node4","node4",bowels,zero,zero,bow_center,matrixbow);
    fibre = new TTUBE("fibre","fibre","Plexi",0.0,fib_r,fib_dz);
    fibre->SetLineColor(28);
    fib_x = zero;
    fib_y = fib_dz*TMath::Sin(delta);
    fib_z = crystal_dz + D + fib_dz*TMath::Cos(delta);
    cout << "deltaD : " << deltaD << endl;
    cout << "new matrixfib" << endl;
//
//!!!!!  Here the program dies for run == 103 !!!!!!  Why ???
//
    matrixfib = new TRotMatrix("Matrixfib","Matrixfib",90.0,0.0,deltaD+90.0,90.0,deltaD,90.0);
    cout << "matrixfib : " << matrixfib << endl;
    cout << "new node5" << endl;
    node5 = new TNode("node5","node5",fibre,fib_x,fib_y,fib_z,matrixfib);
    if (run==96) node1->Draw();
//
//end of loop
//
    fMainPad1->Modified();
    fMainPad1->Update();
    cout << "Type any character to continue ";
    cin >> c;
    delete node5;
    cout << "node5 : " << node5 << endl;
    node5 = 0;
    yy = matrixfib->GetTheta();
    cout << "theta matrixfib : " << yy << endl;
    delete matrixfib;
    cout << "matrixfib : " << matrixfib << endl;
    matrixfib=0;
    delete fibre;
    cout << "fibre : " << fibre << endl;
    fibre = 0;
    delete node4;
    cout << "node4 : " << node4 << endl;
    node4 = 0;
    delete bowels;
    cout << "bowels : " << bowels << endl;
    bowels = 0;
    delta  += deltastep;
    deltaD += deltaDstep;
    cout << "end loop" << endl;
  }
}