[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->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)) 
	    {
		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
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"
5d12ce38	52	#include "TVirtualMC.h"
94de3818	53
21bf7095	54	#include "AliLog.h"
7eb618fa	55	#include "AliDebugVolume.h"
e2afb3b6	56	#include "AliLego.h"
8918e700	57	#include "AliLegoGenerator.h"
5d12ce38	58	#include "AliMC.h"
0742d588	59	#include "AliRun.h"
7eb618fa	60
fe4da5cc	61	ClassImp(AliLego)
fe4da5cc	62
e2afb3b6	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),
39e2942f	79	fDebug(0),
39e2942f	80	fStopped(0)
fe4da5cc	81	{
8918e700	82	//
	83	// Default constructor
	84	//
fe4da5cc	85	}
fe4da5cc	86
e2afb3b6	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),
39e2942f	104	fDebug(0),
39e2942f	105	fStopped(0)
e2afb3b6	106	{
	107	//
	108	// Copy constructor
	109	//
	110	lego.Copy(*this);
	111	}
	112
	113
	114	//_______________________________________________________________________
0b27b8c7	115	AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin,
e2afb3b6	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),
39e2942f	133	fDebug(0),
39e2942f	134	fStopped(0)
fe4da5cc	135	{
8918e700	136	//
	137	// specify the angular limits and the size of the rectangular box
	138	//
0b27b8c7	139	fGener = new AliLegoGenerator(ntheta, thetamin, thetamax,
b13db077	140	nphi, phimin, phimax, rmin, rmax, zmax);
b13db077	141	fHistRadl = new TH2F("hradl","Radiation length map",
0b27b8c7	142	ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077	143	fHistAbso = new TH2F("habso","Interaction length map",
0b27b8c7	144	ntheta,thetamin,thetamax,nphi,phimin,phimax);
b13db077	145	fHistGcm2 = new TH2F("hgcm2","g/cm2 length map",
0b27b8c7	146	ntheta,thetamin,thetamax,nphi,phimin,phimax);
7eb618fa	147	//
7eb618fa	148	fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
7eb618fa	149	fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
21bf7095	150	fDebug = AliDebugLevel();
0b27b8c7	151	}
0b27b8c7	152
e2afb3b6	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),
39e2942f	170	fDebug(0),
39e2942f	171	fStopped(0)
0b27b8c7	172	{
	173	//
	174	// specify the angular limits and the size of the rectangular box
	175	//
e2afb3b6	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	//
b13db077	192
e2afb3b6	193	fVolumesFwd = new TClonesArray("AliDebugVolume",1000);
e2afb3b6	194	fVolumesBwd = new TClonesArray("AliDebugVolume",1000);
21bf7095	195	fDebug = AliDebugLevel();
fe4da5cc	196	}
fe4da5cc	197
e2afb3b6	198	//_______________________________________________________________________
fe4da5cc	199	AliLego::~AliLego()
fe4da5cc	200	{
8918e700	201	//
	202	// Destructor
	203	//
	204	delete fHistRadl;
	205	delete fHistAbso;
	206	delete fHistGcm2;
8918e700	207	delete fGener;
e2afb3b6	208	delete fVolumesFwd;
e2afb3b6	209	delete fVolumesBwd;
fe4da5cc	210	}
fe4da5cc	211
e2afb3b6	212	//_______________________________________________________________________
dffd31ef	213	void AliLego::BeginEvent()
b13db077	214	{
8918e700	215	//
	216	// --- Set to 0 radiation length, absorption length and g/cm2 ---
	217	//
dffd31ef	218	fTotRadl = 0;
	219	fTotAbso = 0;
	220	fTotGcm2 = 0;
39e2942f	221	fStopped = 0;
39e2942f	222
a74d2699	223	if (fDebug) {
21bf7095	224	if (fErrorCondition) ToAliDebug(1, DumpVolumes());
e2afb3b6	225	fVolumesFwd->Delete();
	226	fVolumesBwd->Delete();
	227	fStepsForward = 0;
	228	fStepsBackward = 0;
	229	fErrorCondition = 0;
5d12ce38	230	if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0;
a74d2699	231	}
dffd31ef	232	}
dffd31ef	233
e2afb3b6	234	//_______________________________________________________________________
dffd31ef	235	void AliLego::FinishEvent()
dffd31ef	236	{
8918e700	237	//
	238	// Finish the event and update the histos
	239	//
0b27b8c7	240	Double_t c1, c2;
	241	c1 = fGener->CurCoor1();
	242	c2 = fGener->CurCoor2();
10bd7535	243	fHistRadl->Fill(c2,c1,fTotRadl);
	244	fHistAbso->Fill(c2,c1,fTotAbso);
	245	fHistGcm2->Fill(c2,c1,fTotGcm2);
b13db077	246	}
b13db077	247
e2afb3b6	248	//_______________________________________________________________________
dffd31ef	249	void AliLego::FinishRun()
dffd31ef	250	{
8918e700	251	//
	252	// Store histograms in current Root file
	253	//
dffd31ef	254	fHistRadl->Write();
	255	fHistAbso->Write();
	256	fHistGcm2->Write();
dffd31ef	257
	258	// Delete histograms from memory
	259	fHistRadl->Delete(); fHistRadl=0;
	260	fHistAbso->Delete(); fHistAbso=0;
	261	fHistGcm2->Delete(); fHistGcm2=0;
a74d2699	262	//
21bf7095	263	if (fErrorCondition) ToAliError(DumpVolumes());
dffd31ef	264	}
dffd31ef	265
e2afb3b6	266	//_______________________________________________________________________
6c4904c2	267	void AliLego::Copy(TObject&) const
8918e700	268	{
	269	//
	270	// Copy *this onto lego -- not implemented
	271	//
21bf7095	272	AliFatal("Not implemented!");
8918e700	273	}
dffd31ef	274
e2afb3b6	275	//_______________________________________________________________________
	276	void AliLego::StepManager()
	277	{
	278	//
	279	// called from AliRun::Stepmanager from gustep.
	280	// Accumulate the 3 parameters step by step
	281	//
b75a7a9d	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;
e2afb3b6	301
b75a7a9d	302	if (! fStepBack) {
e2afb3b6	303	// printf("\n volume %s %d", vol, status);
	304	//
	305	// Normal Forward stepping
	306	//
b75a7a9d	307	if (fDebug) {
e2afb3b6	308	// printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step);
	309
	310	//
	311	// store volume if different from previous
	312	//
7eb618fa	313
b75a7a9d	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);
39e2942f	332
b75a7a9d	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
39e2942f	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);
b75a7a9d	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;
39e2942f	373	fStopped = kTRUE;
b75a7a9d	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;
39e2942f	391	// printf("%13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %s %13.3f\n",
39e2942f	392	// pos[2], TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1]), step, a, z, radl, absl, gMC->CurrentVolName(), fTotRadl);
b75a7a9d	393	}
39e2942f	394
b75a7a9d	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->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->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
b13db077	430	}
b13db077	431
e2afb3b6	432	//_______________________________________________________________________
a74d2699	433	void AliLego::DumpVolumes()
a74d2699	434	{
e2afb3b6	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--)
a74d2699	442	{
e2afb3b6	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");
a74d2699	457	}
e2afb3b6	458	printf("\n ==============================================\n");
a74d2699	459	}