[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 "AliConst.h"
#include "AliDebugVolume.h"
#include "AliLego.h"
#include "AliLegoGenerator.h"
#include "AliRun.h"
#include "AliMC.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->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) DumpVolumes();
}

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