[u/mrichter/AliRoot.git] / STEER / AliLego.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
//                            
//  Utility class to evaluate the material budget from
//  a given radius to the surface of an arbitrary cylinder
//  along radial directions from the centre:
// 
//   - radiation length
//   - Interaction length
//   - g/cm2
// 
//  Geantinos are shot in the bins in the fNtheta bins in theta
//  and fNphi bins in phi with specified rectangular limits.
//  The statistics are accumulated per
//    fRadMin < r < fRadMax    and  <0 < z < fZMax
//
//  To activate this option, you can do:
//      Root > gAlice.RunLego();
//  or  Root > .x menu.C  then select menu item "RunLego"
//  Note that when calling gAlice->RunLego, an optional list
//  of arguments may be specified.
//
//Begin_Html
/*
<img src="picts/alilego.gif">
*/
//End_Html
//
//////////////////////////////////////////////////////////////

#include "TClonesArray.h"
#include "TH2.h"
#include "TMath.h"
#include "TString.h"

#include "AliConst.h"
#include "AliDebugVolume.h"
#include "AliLego.h"
#include "AliLegoGenerator.h"
#include "AliRun.h"

ClassImp(AliLego)

//_______________________________________________________________________
AliLego::AliLego():
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // Default constructor
  //
}

//_______________________________________________________________________
AliLego::AliLego(const AliLego &lego):
  TNamed(lego),
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // Copy constructor
  //
  lego.Copy(*this);
}


//_______________________________________________________________________
AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin, 
                 Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
                 Float_t rmin, Float_t rmax, Float_t zmax):
  TNamed("Lego Generator",title),
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // specify the angular limits and the size of the rectangular box
  //
   fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
		       nphi, phimin, phimax, rmin, rmax, zmax);
   fHistRadl = new TH2F("hradl","Radiation length map",    
			ntheta,thetamin,thetamax,nphi,phimin,phimax);
   fHistAbso = new TH2F("habso","Interaction length map",  
			ntheta,thetamin,thetamax,nphi,phimin,phimax);
   fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        
			ntheta,thetamin,thetamax,nphi,phimin,phimax);
//
   fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
   fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
   fDebug          = gAlice->GetDebug();
}

//_______________________________________________________________________
AliLego::AliLego(const char *title, AliLegoGenerator* generator):
  TNamed("Lego Generator",title),
  fGener(0),
  fTotRadl(0),
  fTotAbso(0),
  fTotGcm2(0),
  fHistRadl(0),
  fHistAbso(0),
  fHistGcm2(0),
  fHistReta(0),
  fVolumesFwd(0),
  fVolumesBwd(0),
  fStepBack(0),
  fStepsBackward(0),
  fStepsForward(0),
  fErrorCondition(0),
  fDebug(0)
{
  //
  // specify the angular limits and the size of the rectangular box
  //
  fGener = generator;
  Float_t c1min, c1max, c2min, c2max;
  Int_t n1 = fGener->NCoor1();
  Int_t n2 = fGener->NCoor2();
  
  fGener->Coor1Range(c1min, c1max);
  fGener->Coor2Range(c2min, c2max);   
  
  fHistRadl = new TH2F("hradl","Radiation length map",    
                       n2, c2min, c2max, n1, c1min, c1max);
  fHistAbso = new TH2F("habso","Interaction length map",  
                       n2, c2min, c2max, n1, c1min, c1max);
  fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",        
                       n2, c2min, c2max, n1, c1min, c1max);
  //
  //

  fVolumesFwd     = new TClonesArray("AliDebugVolume",1000);
  fVolumesBwd     = new TClonesArray("AliDebugVolume",1000);
  fDebug          = gAlice->GetDebug();
}

//_______________________________________________________________________
AliLego::~AliLego()
{
  //
  // Destructor
  //
  delete fHistRadl;
  delete fHistAbso;
  delete fHistGcm2;
  delete fGener;
  delete fVolumesFwd;
  delete fVolumesBwd;
}

//_______________________________________________________________________
void AliLego::BeginEvent()
{
  //
  // --- Set to 0 radiation length, absorption length and g/cm2 ---
  //
  fTotRadl = 0;
  fTotAbso = 0;
  fTotGcm2 = 0;

  if (fDebug) {
    if (fErrorCondition) DumpVolumes();
    fVolumesFwd->Delete();
    fVolumesBwd->Delete();
    fStepsForward    = 0;
    fStepsBackward   = 0;		  
    fErrorCondition  = 0;
    if (gAlice->GetCurrentTrackNumber() == 0) fStepBack = 0;
  }
}

//_______________________________________________________________________
void AliLego::FinishEvent()
{
  //
  // Finish the event and update the histos
  //
  Double_t c1, c2;
  c1 = fGener->CurCoor1();
  c2 = fGener->CurCoor2();
  fHistRadl->Fill(c2,c1,fTotRadl);
  fHistAbso->Fill(c2,c1,fTotAbso);
  fHistGcm2->Fill(c2,c1,fTotGcm2);
}

//_______________________________________________________________________
void AliLego::FinishRun()
{
  //
  // Store histograms in current Root file
  //
  fHistRadl->Write();
  fHistAbso->Write();
  fHistGcm2->Write();

  // Delete histograms from memory
  fHistRadl->Delete(); fHistRadl=0;
  fHistAbso->Delete(); fHistAbso=0;
  fHistGcm2->Delete(); fHistGcm2=0;
  //
  if (fErrorCondition) DumpVolumes();
}

//_______________________________________________________________________
void AliLego::Copy(AliLego&) const
{
  //
  // Copy *this onto lego -- not implemented
  //
  Fatal("Copy","Not implemented!\n");
}

//_______________________________________________________________________
void AliLego::StepManager() 
{
  //
  // called from AliRun::Stepmanager from gustep.
  // Accumulate the 3 parameters step by step
  //
  static Float_t t;
  Float_t a,z,dens,radl,absl;
  Int_t i, id, copy;
  const char* vol;
  static Float_t vect[3], dir[3];
  
  TString tmp1, tmp2;
  copy = 1;
  id  = gMC->CurrentVolID(copy);
  vol = gMC->VolName(id);
  Float_t step  = gMC->TrackStep();
  
  TLorentzVector pos, mom; 
  gMC->TrackPosition(pos);  
  gMC->TrackMomentum(mom);
  
  Int_t status = 0;
  if (gMC->IsTrackEntering()) status = 1;
  if (gMC->IsTrackExiting())  status = 2; 
  
  if (! fStepBack) {
    //      printf("\n volume %s %d", vol, status);      
    //
    // Normal Forward stepping
    //
    if (fDebug) {
      //	  printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);	  
      
      //
      // store volume if different from previous
      //
	  
	  TClonesArray &lvols = *fVolumesFwd;
	  if (fStepsForward > 0) {
        AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsForward-1]);
        if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
		  fStepsForward -= 2;
		  fVolumesFwd->RemoveAt(fStepsForward);
		  fVolumesFwd->RemoveAt(fStepsForward+1);		  
        }
	  }
      
	  new(lvols[fStepsForward++]) 
        AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
      
    } // Debug
    //
    // Get current material properties
    
    gMC->CurrentMaterial(a,z,dens,radl,absl);
    
    if (z < 1) return;
    
    // --- See if we have to stop now
    if (TMath::Abs(pos[2]) > fGener->ZMax()  || 
        pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) {
      if (!gMC->IsNewTrack()) {
        // Not the first step, add past contribution
        fTotAbso += t/absl;
        fTotRadl += t/radl;
        fTotGcm2 += t*dens;
        if (fDebug) {
          //
          //  generate "mirror" particle flying back
          //
          fStepsBackward = fStepsForward;
          
          Float_t pmom[3], orig[3];
          Float_t polar[3] = {0.,0.,0.};
          Int_t ntr;
          pmom[0] = -dir[0];
          pmom[1] = -dir[1];	   
          pmom[2] = -dir[2];
          orig[0] =  vect[0];
          orig[1] =  vect[1];	   
          orig[2] =  vect[2];
          
          gAlice->PushTrack(1, gAlice->GetCurrentTrackNumber(), 
                           0, pmom, orig, polar, 0., kPNoProcess, ntr);
        } // debug
        
      } // not a new track !
      
      if (fDebug) fStepBack = 1;
      gMC->StopTrack();
      return;
    } // outside scoring region ?
    
    // --- See how long we have to go
    for(i=0;i<3;++i) {
      vect[i]=pos[i];
      dir[i]=mom[i];
    }
    
    t  = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
    
    if(step) {
      
      fTotAbso += step/absl;
      fTotRadl += step/radl;
      fTotGcm2 += step*dens;
    }
    
  } else {
      if (fDebug) {
        //
        // Geometry debugging
        // Fly back and compare volume sequence
        //
        TClonesArray &lvols = *fVolumesBwd;
        if (fStepsBackward < fStepsForward) {
	      AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[fStepsBackward]);
	      if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
            fStepsBackward += 2;
            fVolumesBwd->RemoveAt(fStepsBackward-1);
            fVolumesBwd->RemoveAt(fStepsBackward-2);		  
	      }
        } 
        
        fStepsBackward--;
        //	  printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);	  
        if (fStepsBackward < 0) {
	      gMC->StopTrack();
	      fStepBack = 0;
	      return;
        }
        
        new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
        
        AliDebugVolume* tmp = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[fStepsBackward]);
        if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition)) 
          {
            printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n", 
                   fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
            fErrorCondition = 1;
          } 
      } // Debug
  } // bwd/fwd
}

//_______________________________________________________________________
void AliLego::DumpVolumes()
{
  //
  // Dump volume sequence in case of error
  //
  printf("\n Dumping Volume Sequence:");
  printf("\n ==============================================");
  
  for (Int_t i = fStepsForward-1; i>=0; i--)
    {
      AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume*>((*fVolumesFwd)[i]);
      AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume*>((*fVolumesBwd)[i]);
      if (tmp1)
        printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
               , i, 
               tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(), 
               tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
      if (tmp2 && i>= fStepsBackward)
        printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
               , i, 
               tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(), 
               tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
      
      printf("\n ............................................................................\n");
    }
  printf("\n ==============================================\n");
}
Commit	Line	Data
4c039060	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	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	**************************************************************************/
	15
acd84897	16	/* $Id$ */
4c039060	17
fe4da5cc	18	//////////////////////////////////////////////////////////////
	19	//////////////////////////////////////////////////////////////
	20	//
	21	// Utility class to evaluate the material budget from
	22	// a given radius to the surface of an arbitrary cylinder
	23	// along radial directions from the centre:
	24	//
	25	// - radiation length
	26	// - Interaction length
	27	// - g/cm2
	28	//
	29	// Geantinos are shot in the bins in the fNtheta bins in theta
	30	// and fNphi bins in phi with specified rectangular limits.
	31	// The statistics are accumulated per
	32	// fRadMin < r < fRadMax and <0 < z < fZMax
	33	//
	34	// To activate this option, you can do:
	35	// Root > gAlice.RunLego();
	36	// or Root > .x menu.C then select menu item "RunLego"
	37	// Note that when calling gAlice->RunLego, an optional list
	38	// of arguments may be specified.
	39	//
	40	//Begin_Html
	41	/*
1439f98e	42	<img src="picts/alilego.gif">
fe4da5cc	43	*/
	44	//End_Html
	45	//
	46	//////////////////////////////////////////////////////////////
	47
e2afb3b6	48	#include "TClonesArray.h"
e2afb3b6	49	#include "TH2.h"
fe4da5cc	50	#include "TMath.h"
e2afb3b6	51	#include "TString.h"
94de3818	52
e2afb3b6	53	#include "AliConst.h"
7eb618fa	54	#include "AliDebugVolume.h"
e2afb3b6	55	#include "AliLego.h"
8918e700	56	#include "AliLegoGenerator.h"
e2afb3b6	57	#include "AliRun.h"
7eb618fa	58
fe4da5cc	59	ClassImp(AliLego)
fe4da5cc	60
e2afb3b6	61	//_______________________________________________________________________
	62	AliLego::AliLego():
	63	fGener(0),
	64	fTotRadl(0),
	65	fTotAbso(0),
	66	fTotGcm2(0),
	67	fHistRadl(0),
	68	fHistAbso(0),
	69	fHistGcm2(0),
	70	fHistReta(0),
	71	fVolumesFwd(0),
	72	fVolumesBwd(0),
	73	fStepBack(0),
	74	fStepsBackward(0),
	75	fStepsForward(0),
	76	fErrorCondition(0),
	77	fDebug(0)
fe4da5cc	78	{
8918e700	79	//
	80	// Default constructor
	81	//
fe4da5cc	82	}
fe4da5cc	83
e2afb3b6	84	//_______________________________________________________________________
	85	AliLego::AliLego(const AliLego &lego):
	86	TNamed(lego),
	87	fGener(0),
	88	fTotRadl(0),
	89	fTotAbso(0),
	90	fTotGcm2(0),
	91	fHistRadl(0),
	92	fHistAbso(0),
	93	fHistGcm2(0),
	94	fHistReta(0),
	95	fVolumesFwd(0),
	96	fVolumesBwd(0),
	97	fStepBack(0),
	98	fStepsBackward(0),
	99	fStepsForward(0),
	100	fErrorCondition(0),
	101	fDebug(0)
	102	{
	103	//
	104	// Copy constructor
	105	//
	106	lego.Copy(*this);
	107	}
	108
	109
	110	//_______________________________________________________________________
0b27b8c7	111	AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
e2afb3b6	112	Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax,
	113	Float_t rmin, Float_t rmax, Float_t zmax):
	114	TNamed("Lego Generator",title),
	115	fGener(0),
	116	fTotRadl(0),
	117	fTotAbso(0),
	118	fTotGcm2(0),
	119	fHistRadl(0),
	120	fHistAbso(0),
	121	fHistGcm2(0),
	122	fHistReta(0),
	123	fVolumesFwd(0),
	124	fVolumesBwd(0),
	125	fStepBack(0),
	126	fStepsBackward(0),
	127	fStepsForward(0),
	128	fErrorCondition(0),
	129	fDebug(0)
fe4da5cc	130	{
8918e700	131	//
	132	// specify the angular limits and the size of the rectangular box
	133	//
0b27b8c7	134	fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
b13db077	135	nphi, phimin, phimax, rmin, rmax, zmax);
b13db077	136	fHistRadl = new TH2F("hradl","Radiation length map",
0b27b8c7	137	ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077	138	fHistAbso = new TH2F("habso","Interaction length map",
0b27b8c7	139	ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077	140	fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
0b27b8c7	141	ntheta,thetamin,thetamax,nphi,phimin,phimax);
7eb618fa	142	//
7eb618fa	143	fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
7eb618fa	144	fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
19237f5f	145	fDebug = gAlice->GetDebug();
0b27b8c7	146	}
0b27b8c7	147
e2afb3b6	148	//_______________________________________________________________________
	149	AliLego::AliLego(const char title, AliLegoGenerator generator):
	150	TNamed("Lego Generator",title),
	151	fGener(0),
	152	fTotRadl(0),
	153	fTotAbso(0),
	154	fTotGcm2(0),
	155	fHistRadl(0),
	156	fHistAbso(0),
	157	fHistGcm2(0),
	158	fHistReta(0),
	159	fVolumesFwd(0),
	160	fVolumesBwd(0),
	161	fStepBack(0),
	162	fStepsBackward(0),
	163	fStepsForward(0),
	164	fErrorCondition(0),
	165	fDebug(0)
0b27b8c7	166	{
	167	//
	168	// specify the angular limits and the size of the rectangular box
	169	//
e2afb3b6	170	fGener = generator;
	171	Float_t c1min, c1max, c2min, c2max;
	172	Int_t n1 = fGener->NCoor1();
	173	Int_t n2 = fGener->NCoor2();
	174
	175	fGener->Coor1Range(c1min, c1max);
	176	fGener->Coor2Range(c2min, c2max);
	177
	178	fHistRadl = new TH2F("hradl","Radiation length map",
	179	n2, c2min, c2max, n1, c1min, c1max);
	180	fHistAbso = new TH2F("habso","Interaction length map",
	181	n2, c2min, c2max, n1, c1min, c1max);
	182	fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
	183	n2, c2min, c2max, n1, c1min, c1max);
	184	//
	185	//
b13db077	186
e2afb3b6	187	fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
	188	fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
	189	fDebug = gAlice->GetDebug();
fe4da5cc	190	}
fe4da5cc	191
e2afb3b6	192	//_______________________________________________________________________
fe4da5cc	193	AliLego::~AliLego()
fe4da5cc	194	{
8918e700	195	//
	196	// Destructor
	197	//
	198	delete fHistRadl;
	199	delete fHistAbso;
	200	delete fHistGcm2;
8918e700	201	delete fGener;
e2afb3b6	202	delete fVolumesFwd;
e2afb3b6	203	delete fVolumesBwd;
fe4da5cc	204	}
fe4da5cc	205
e2afb3b6	206	//_______________________________________________________________________
dffd31ef	207	void AliLego::BeginEvent()
b13db077	208	{
8918e700	209	//
	210	// --- Set to 0 radiation length, absorption length and g/cm2 ---
	211	//
dffd31ef	212	fTotRadl = 0;
	213	fTotAbso = 0;
	214	fTotGcm2 = 0;
19237f5f	215
a74d2699	216	if (fDebug) {
e2afb3b6	217	if (fErrorCondition) DumpVolumes();
	218	fVolumesFwd->Delete();
	219	fVolumesBwd->Delete();
	220	fStepsForward = 0;
	221	fStepsBackward = 0;
	222	fErrorCondition = 0;
642f15cf	223	if (gAlice->GetCurrentTrackNumber() == 0) fStepBack = 0;
a74d2699	224	}
dffd31ef	225	}
dffd31ef	226
e2afb3b6	227	//_______________________________________________________________________
dffd31ef	228	void AliLego::FinishEvent()
dffd31ef	229	{
8918e700	230	//
	231	// Finish the event and update the histos
	232	//
0b27b8c7	233	Double_t c1, c2;
	234	c1 = fGener->CurCoor1();
	235	c2 = fGener->CurCoor2();
10bd7535	236	fHistRadl->Fill(c2,c1,fTotRadl);
	237	fHistAbso->Fill(c2,c1,fTotAbso);
	238	fHistGcm2->Fill(c2,c1,fTotGcm2);
b13db077	239	}
b13db077	240
e2afb3b6	241	//_______________________________________________________________________
dffd31ef	242	void AliLego::FinishRun()
dffd31ef	243	{
8918e700	244	//
	245	// Store histograms in current Root file
	246	//
dffd31ef	247	fHistRadl->Write();
	248	fHistAbso->Write();
	249	fHistGcm2->Write();
dffd31ef	250
	251	// Delete histograms from memory
	252	fHistRadl->Delete(); fHistRadl=0;
	253	fHistAbso->Delete(); fHistAbso=0;
	254	fHistGcm2->Delete(); fHistGcm2=0;
a74d2699	255	//
a74d2699	256	if (fErrorCondition) DumpVolumes();
dffd31ef	257	}
dffd31ef	258
e2afb3b6	259	//_______________________________________________________________________
e2afb3b6	260	void AliLego::Copy(AliLego&) const
8918e700	261	{
	262	//
	263	// Copy *this onto lego -- not implemented
	264	//
	265	Fatal("Copy","Not implemented!\n");
	266	}
dffd31ef	267
e2afb3b6	268	//_______________________________________________________________________
	269	void AliLego::StepManager()
	270	{
	271	//
	272	// called from AliRun::Stepmanager from gustep.
	273	// Accumulate the 3 parameters step by step
	274	//
	275	static Float_t t;
	276	Float_t a,z,dens,radl,absl;
	277	Int_t i, id, copy;
	278	const char* vol;
	279	static Float_t vect[3], dir[3];
	280
	281	TString tmp1, tmp2;
	282	copy = 1;
	283	id = gMC->CurrentVolID(copy);
	284	vol = gMC->VolName(id);
	285	Float_t step = gMC->TrackStep();
	286
	287	TLorentzVector pos, mom;
	288	gMC->TrackPosition(pos);
	289	gMC->TrackMomentum(mom);
	290
	291	Int_t status = 0;
	292	if (gMC->IsTrackEntering()) status = 1;
	293	if (gMC->IsTrackExiting()) status = 2;
	294
7eb618fa	295	if (! fStepBack) {
e2afb3b6	296	// printf("\n volume %s %d", vol, status);
	297	//
	298	// Normal Forward stepping
	299	//
	300	if (fDebug) {
	301	// printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
	302
	303	//
	304	// store volume if different from previous
	305	//
7eb618fa	306
a74d2699	307	TClonesArray &lvols = *fVolumesFwd;
7eb618fa	308	if (fStepsForward > 0) {
e2afb3b6	309	AliDebugVolume* tmp = dynamic_cast<AliDebugVolume>((fVolumesFwd)[fStepsForward-1]);
e2afb3b6	310	if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
a74d2699	311	fStepsForward -= 2;
5f3dc1c9	312	fVolumesFwd->RemoveAt(fStepsForward);
5f3dc1c9	313	fVolumesFwd->RemoveAt(fStepsForward+1);
e2afb3b6	314	}
7eb618fa	315	}
e2afb3b6	316
a74d2699	317	new(lvols[fStepsForward++])
e2afb3b6	318	AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
7eb618fa	319
e2afb3b6	320	} // Debug
	321	//
	322	// Get current material properties
	323
	324	gMC->CurrentMaterial(a,z,dens,radl,absl);
	325
	326	if (z < 1) return;
	327
	328	// --- See if we have to stop now
	329	if (TMath::Abs(pos[2]) > fGener->ZMax() \|\|
	330	pos[0]pos[0] +pos[1]pos[1] > fGener->RadMax()*fGener->RadMax()) {
	331	if (!gMC->IsNewTrack()) {
	332	// Not the first step, add past contribution
	333	fTotAbso += t/absl;
	334	fTotRadl += t/radl;
	335	fTotGcm2 += t*dens;
	336	if (fDebug) {
	337	//
	338	// generate "mirror" particle flying back
	339	//
	340	fStepsBackward = fStepsForward;
	341
	342	Float_t pmom[3], orig[3];
	343	Float_t polar[3] = {0.,0.,0.};
	344	Int_t ntr;
	345	pmom[0] = -dir[0];
	346	pmom[1] = -dir[1];
	347	pmom[2] = -dir[2];
	348	orig[0] = vect[0];
	349	orig[1] = vect[1];
	350	orig[2] = vect[2];
	351
642f15cf	352	gAlice->PushTrack(1, gAlice->GetCurrentTrackNumber(),
e2afb3b6	353	0, pmom, orig, polar, 0., kPNoProcess, ntr);
	354	} // debug
	355
	356	} // not a new track !
7eb618fa	357
e2afb3b6	358	if (fDebug) fStepBack = 1;
	359	gMC->StopTrack();
	360	return;
	361	} // outside scoring region ?
	362
	363	// --- See how long we have to go
	364	for(i=0;i<3;++i) {
	365	vect[i]=pos[i];
	366	dir[i]=mom[i];
	367	}
	368
	369	t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax());
	370
	371	if(step) {
	372
	373	fTotAbso += step/absl;
	374	fTotRadl += step/radl;
	375	fTotGcm2 += step*dens;
	376	}
	377
7eb618fa	378	} else {
7eb618fa	379	if (fDebug) {
e2afb3b6	380	//
	381	// Geometry debugging
	382	// Fly back and compare volume sequence
	383	//
	384	TClonesArray &lvols = *fVolumesBwd;
	385	if (fStepsBackward < fStepsForward) {
	386	AliDebugVolume* tmp = dynamic_cast<AliDebugVolume>((fVolumesBwd)[fStepsBackward]);
	387	if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) {
	388	fStepsBackward += 2;
	389	fVolumesBwd->RemoveAt(fStepsBackward-1);
	390	fVolumesBwd->RemoveAt(fStepsBackward-2);
7eb618fa	391	}
e2afb3b6	392	}
	393
	394	fStepsBackward--;
	395	// printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition);
	396	if (fStepsBackward < 0) {
7eb618fa	397	gMC->StopTrack();
a74d2699	398	fStepBack = 0;
7eb618fa	399	return;
e2afb3b6	400	}
	401
	402	new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status);
	403
	404	AliDebugVolume* tmp = dynamic_cast<AliDebugVolume>((fVolumesFwd)[fStepsBackward]);
	405	if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition))
	406	{
	407	printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n",
	408	fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step);
	409	fErrorCondition = 1;
	410	}
7eb618fa	411	} // Debug
7eb618fa	412	} // bwd/fwd
b13db077	413	}
b13db077	414
e2afb3b6	415	//_______________________________________________________________________
a74d2699	416	void AliLego::DumpVolumes()
a74d2699	417	{
e2afb3b6	418	//
	419	// Dump volume sequence in case of error
	420	//
	421	printf("\n Dumping Volume Sequence:");
	422	printf("\n ==============================================");
	423
	424	for (Int_t i = fStepsForward-1; i>=0; i--)
a74d2699	425	{
e2afb3b6	426	AliDebugVolume* tmp1 = dynamic_cast<AliDebugVolume>((fVolumesFwd)[i]);
	427	AliDebugVolume* tmp2 = dynamic_cast<AliDebugVolume>((fVolumesBwd)[i]);
	428	if (tmp1)
	429	printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
	430	, i,
	431	tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(),
	432	tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status());
	433	if (tmp2 && i>= fStepsBackward)
	434	printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n"
	435	, i,
	436	tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(),
	437	tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status());
	438
	439	printf("\n ............................................................................\n");
a74d2699	440	}
e2afb3b6	441	printf("\n ==============================================\n");
a74d2699	442	}