[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 "TVirtualMC.h"

#include "AliLog.h"
#include "AliDebugVolume.h"
#include "AliLego.h"
#include "AliLegoGenerator.h"
#include "AliMC.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),
  fStopped(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),
  fStopped(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),
  fStopped(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          = AliDebugLevel();
}

//_______________________________________________________________________
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),
  fStopped(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          = AliDebugLevel();
}

//_______________________________________________________________________
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;
  fStopped = 0;
  
  if (fDebug) {
    if (fErrorCondition) ToAliDebug(1, DumpVolumes());
    fVolumesFwd->Delete();
    fVolumesBwd->Delete();
    fStepsForward    = 0;
    fStepsBackward   = 0;		  
    fErrorCondition  = 0;
    if (gAlice->GetMCApp()->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) ToAliError(DumpVolumes());
}

//_______________________________________________________________________
void AliLego::Copy(TObject&) const
{
  //
  // Copy *this onto lego -- not implemented
  //
  AliFatal("Not implemented!");
}

//_______________________________________________________________________
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
		if (!fStopped) {
		    if (absl) fTotAbso += t/absl;
		    if (radl) fTotRadl += t/radl;
		    fTotGcm2 += t*dens;
		}
		
//		printf("We will stop now %5d %13.3f !\n", fStopped, t);
//		printf("%13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %s %13.3f\n",
//		       pos[2], TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1]), step, a, z, radl, absl, gMC->CurrentVolName(), fTotRadl);
		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->GetMCApp()->PushTrack(1, gAlice->GetMCApp()->GetCurrentTrackNumber(), 
						  0, pmom, orig, polar, 0., kPNoProcess, ntr);
		} // debug
		
	    } // not a new track !
	    
	    if (fDebug) fStepBack = 1;
	    fStopped = kTRUE;
	    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) {
	    
	    if (absl) fTotAbso += step/absl;
	    if (radl) fTotRadl += step/radl;
	    fTotGcm2 += step*dens;
//	     printf("%13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %s %13.3f\n",
//	     pos[2],  TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1]), step, a, z, radl, absl, gMC->CurrentVolName(), fTotRadl);
	}

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