[u/mrichter/AliRoot.git] / PWG1 / AliMCInfo.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.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////
/*

Origin: marian.ivanov@cern.ch
Container classes with MC infomation

The AliMCInfo contains the information about the particles properties 
during transportation throuch ALICE Detector

The base Information :
TParticle - fParticle             -  properties of the particle at creation point
AliTrackReference - fXXXRefernces -  TClonesArray of refernces in differnt detectors
fNXXXRef                          -  number of the track refernces in differnt detectors
AliTPCdigitRow    - fTPCRow       -  the map of the hitted rows - (will be repalced by TBits)
fRowsWith*                        - number of rows hitted by particle
fMCtracks         -               - number of turn over of the track inside of the TPC

++++
some additional information usable for tree draw - TO SPEED UP tree queries
IMPORTANT FOR PROOF FAST PROTOTYPING ANALYSIS 
                                      

*/

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdio.h>
//ROOT includes
#include "Rtypes.h"
#include "TClonesArray.h"
//ALIROOT includes
#include "AliTrackReference.h"
#include "AliMCInfo.h" 
#endif

//
// 

ClassImp(AliTPCdigitRow)
ClassImp(AliMCInfo)


////////////////////////////////////////////////////////////////////////
AliMCInfo::AliMCInfo():
  fTrackRef(),
  fTrackRefOut(),
  fTRdecay(),
  fPrimPart(0),
  fParticle(),
  fMass(0),
  fCharge(0),
  fLabel(0),
  fEventNr(),
  fMCtracks(),
  fPdg(0),
  fTPCdecay(0),
  fRowsWithDigitsInn(0),
  fRowsWithDigits(0),
  fRowsTrackLength(0),
  fPrim(0),
  fTPCRow(), 
  fNTPCRef(0),                    // tpc references counter
  fNITSRef(0),                    // ITS references counter
  fNTRDRef(0),                    // TRD references counter
  fNTOFRef(0),                    // TOF references counter
  fTPCReferences(0),
  fITSReferences(0),
  fTRDReferences(0),
  fTOFReferences(0)
{
  //
  // Default constructor
  //
  fTPCReferences  = new TClonesArray("AliTrackReference",10);
  fITSReferences  = new TClonesArray("AliTrackReference",10);
  fTRDReferences  = new TClonesArray("AliTrackReference",10);
  fTOFReferences  = new TClonesArray("AliTrackReference",10);
  fTRdecay.SetTrack(-1);
  fCharge = 0;
}

AliMCInfo::AliMCInfo(const AliMCInfo& info):
  TObject(),
  fTrackRef(info.fTrackRef),
  fTrackRefOut(info.fTrackRefOut),
  fTRdecay(info.GetTRdecay()),
  fPrimPart(info.fPrimPart),
  fParticle(info.fParticle),
  fMass(info.GetMass()),
  fCharge(info.fCharge),
  fLabel(info.fLabel),
  fEventNr(info.fEventNr),
  fMCtracks(info.fMCtracks),
  fPdg(info.fPdg),
  fTPCdecay(info.fTPCdecay),
  fRowsWithDigitsInn(info.fRowsWithDigitsInn),
  fRowsWithDigits(info.fRowsWithDigits),
  fRowsTrackLength(info.fRowsTrackLength),
  fPrim(info.fPrim),
  fTPCRow(info.fTPCRow), 
  fNTPCRef(info.fNTPCRef),                    // tpc references counter
  fNITSRef(info.fNITSRef),                    // ITS references counter
  fNTRDRef(info.fNTRDRef),                    // TRD references counter
  fNTOFRef(info.fNTOFRef),                    // TOF references counter
  fTPCReferences(0),
  fITSReferences(0),
  fTRDReferences(0),
  fTOFReferences(0)
{
  //
  // copy constructor
  //
  fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
  fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
  fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
  fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
}


AliMCInfo::~AliMCInfo()
{
  //
  // Destructor of the class
  //
  if (fTPCReferences) {
    delete fTPCReferences;
  }
  if (fITSReferences){
    delete fITSReferences;
  }
  if (fTRDReferences){    
    delete fTRDReferences;  
  }
  if (fTOFReferences){    
    delete fTOFReferences;  
  }

}


void AliMCInfo::Update()
{
  //
  // Update MC info
  // Calculates some derived variables
  //
  //
  fMCtracks =1;
  if (!fTPCReferences) {
    fNTPCRef =0;
    return;
  }
  Float_t direction=1;
  //Float_t rlast=0;
  fNTPCRef = fTPCReferences->GetEntriesFast();
  fNITSRef = fITSReferences->GetEntriesFast();
  fNTRDRef = fTRDReferences->GetEntriesFast();
  fNTOFRef = fTOFReferences->GetEntriesFast();
  
  for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
    AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
    //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
    Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
    if (iref==0) direction = newdirection;
    if ( newdirection*direction<0){
      //changed direction
      direction = newdirection;
      fMCtracks+=1;
    }
    //rlast=r;			    
  }
  //
  // decay info
  fTPCdecay=kFALSE;
  if (fTRdecay.GetTrack()>0){
    fDecayCoord[0] = fTRdecay.X();
    fDecayCoord[1] = fTRdecay.Y();
    fDecayCoord[2] = fTRdecay.Z();
    if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
      fTPCdecay=kTRUE;     
    }
    else{
      fDecayCoord[0] = 0;
      fDecayCoord[1] = 0;
      fDecayCoord[2] = 0;
    }
  }
  //
  //
  //digits information update
  fRowsWithDigits    = fTPCRow.RowsOn();    
  fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
  fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
  //
  //
  // calculate primary ionization per cm  
  if (fParticle.GetPDG()){
    fMass = fParticle.GetMass();  
    fCharge = fParticle.GetPDG()->Charge();
    if (fTPCReferences->GetEntriesFast()>0){
      fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
    }
    if (fMass>0){
      Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
			      fTrackRef.Py()*fTrackRef.Py()+
			      fTrackRef.Pz()*fTrackRef.Pz());    
      if (p>0.001){
	Float_t betagama = p /fMass;
	fPrim = TPCBetheBloch(betagama);
      }else fPrim=0;
    }
  }else{
    fMass =0;
    fPrim =0;
  }  
}


////////////////////////////////////////////////////////////////////////
AliTPCdigitRow::AliTPCdigitRow()
{
  Reset();
}
////////////////////////////////////////////////////////////////////////
AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
{
  for (Int_t i = 0; i<32; i++) fDig[i] = digOld.fDig[i];
  return (*this);
}
////////////////////////////////////////////////////////////////////////
void AliTPCdigitRow::SetRow(Int_t row) 
{
  //
  // set bit mask for given row
  //
  if (row >= 8*32) {
    //    cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
    return;
  }
  Int_t iC = row/8;
  Int_t iB = row%8;
  SETBIT(fDig[iC],iB);
}

////////////////////////////////////////////////////////////////////////
Bool_t AliTPCdigitRow::TestRow(Int_t row) const
{
//
// return kTRUE if row is on
//
  Int_t iC = row/8;
  Int_t iB = row%8;
  return TESTBIT(fDig[iC],iB);
}
////////////////////////////////////////////////////////////////////////
Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
{
//
// returns number of rows with a digit  
// count only rows less equal row number upto
//
  Int_t total = 0;
  for (Int_t i = 0; i<32; i++) {
    for (Int_t j = 0; j < 8; j++) {
      if (i*8+j > upto) return total;
      if (TESTBIT(fDig[i],j))  total++;
    }
  }
  return total;
}
////////////////////////////////////////////////////////////////////////
void AliTPCdigitRow::Reset()
{
//
// resets all rows to zero
//
  for (Int_t i = 0; i<32; i++) {
    fDig[i] <<= 8;
  }
}
////////////////////////////////////////////////////////////////////////
Int_t AliTPCdigitRow::Last() const
{
//
// returns the last row number with a digit
// returns -1 if now digits 
//
  for (Int_t i = 32-1; i>=0; i--) {
    for (Int_t j = 7; j >= 0; j--) {
      if TESTBIT(fDig[i],j) return i*8+j;
    }
  }
  return -1;
}
////////////////////////////////////////////////////////////////////////
Int_t AliTPCdigitRow::First() const
{
//
// returns the first row number with a digit
// returns -1 if now digits 
//
  for (Int_t i = 0; i<32; i++) {
    for (Int_t j = 0; j < 8; j++) {
      if (TESTBIT(fDig[i],j)) return i*8+j;
    }
  }
  return -1;
}


//_____________________________________________________________________________
Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
{
  //
  // Bethe-Bloch energy loss formula
  //
  const Double_t kp1=0.76176e-1;
  const Double_t kp2=10.632;
  const Double_t kp3=0.13279e-4;
  const Double_t kp4=1.8631;
  const Double_t kp5=1.9479;

  Double_t dbg = (Double_t) bg;

  Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);

  Double_t aa = TMath::Power(beta,kp4);
  Double_t bb = TMath::Power(1./dbg,kp5);

  bb=TMath::Log(kp3+bb);
  
  return ((Float_t)((kp2-aa-bb)*kp1/aa));
}
Commit	Line	Data
5c7ef659	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
	17	///////////////////////////////////////////////////////////////////////////
	18	/*
	19
	20	Origin: marian.ivanov@cern.ch
	21	Container classes with MC infomation
	22
7ff5e4cf	23	The AliMCInfo contains the information about the particles properties
	24	during transportation throuch ALICE Detector
	25
	26	The base Information :
	27	TParticle - fParticle - properties of the particle at creation point
	28	AliTrackReference - fXXXRefernces - TClonesArray of refernces in differnt detectors
	29	fNXXXRef - number of the track refernces in differnt detectors
	30	AliTPCdigitRow - fTPCRow - the map of the hitted rows - (will be repalced by TBits)
	31	fRowsWith* - number of rows hitted by particle
	32	fMCtracks - - number of turn over of the track inside of the TPC
	33
	34	++++
	35	some additional information usable for tree draw - TO SPEED UP tree queries
	36	IMPORTANT FOR PROOF FAST PROTOTYPING ANALYSIS
	37
	38
	39
	40
5c7ef659	41	*/
	42
	43	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	44	#include <stdio.h>
5c7ef659	45	//ROOT includes
5c7ef659	46	#include "Rtypes.h"
7ff5e4cf	47	#include "TClonesArray.h"
5c7ef659	48	//ALIROOT includes
5c7ef659	49	#include "AliTrackReference.h"
5c7ef659	50	#include "AliMCInfo.h"
7ff5e4cf	51	#endif
7ff5e4cf	52
5c7ef659	53	//
	54	//
	55
	56	ClassImp(AliTPCdigitRow)
	57	ClassImp(AliMCInfo)
	58
	59
	60
	61	////////////////////////////////////////////////////////////////////////
	62	AliMCInfo::AliMCInfo():
	63	fTrackRef(),
	64	fTrackRefOut(),
	65	fTRdecay(),
	66	fPrimPart(0),
	67	fParticle(),
	68	fMass(0),
	69	fCharge(0),
	70	fLabel(0),
	71	fEventNr(),
	72	fMCtracks(),
	73	fPdg(0),
	74	fTPCdecay(0),
	75	fRowsWithDigitsInn(0),
	76	fRowsWithDigits(0),
	77	fRowsTrackLength(0),
	78	fPrim(0),
	79	fTPCRow(),
	80	fNTPCRef(0), // tpc references counter
	81	fNITSRef(0), // ITS references counter
	82	fNTRDRef(0), // TRD references counter
	83	fNTOFRef(0), // TOF references counter
	84	fTPCReferences(0),
	85	fITSReferences(0),
	86	fTRDReferences(0),
	87	fTOFReferences(0)
	88	{
7ff5e4cf	89	//
	90	// Default constructor
	91	//
5c7ef659	92	fTPCReferences = new TClonesArray("AliTrackReference",10);
	93	fITSReferences = new TClonesArray("AliTrackReference",10);
	94	fTRDReferences = new TClonesArray("AliTrackReference",10);
	95	fTOFReferences = new TClonesArray("AliTrackReference",10);
	96	fTRdecay.SetTrack(-1);
	97	fCharge = 0;
	98	}
	99
	100	AliMCInfo::AliMCInfo(const AliMCInfo& info):
	101	TObject(),
	102	fTrackRef(info.fTrackRef),
	103	fTrackRefOut(info.fTrackRefOut),
	104	fTRdecay(info.GetTRdecay()),
	105	fPrimPart(info.fPrimPart),
	106	fParticle(info.fParticle),
	107	fMass(info.GetMass()),
	108	fCharge(info.fCharge),
	109	fLabel(info.fLabel),
	110	fEventNr(info.fEventNr),
	111	fMCtracks(info.fMCtracks),
	112	fPdg(info.fPdg),
	113	fTPCdecay(info.fTPCdecay),
	114	fRowsWithDigitsInn(info.fRowsWithDigitsInn),
	115	fRowsWithDigits(info.fRowsWithDigits),
	116	fRowsTrackLength(info.fRowsTrackLength),
	117	fPrim(info.fPrim),
	118	fTPCRow(info.fTPCRow),
	119	fNTPCRef(info.fNTPCRef), // tpc references counter
	120	fNITSRef(info.fNITSRef), // ITS references counter
	121	fNTRDRef(info.fNTRDRef), // TRD references counter
	122	fNTOFRef(info.fNTOFRef), // TOF references counter
	123	fTPCReferences(0),
	124	fITSReferences(0),
	125	fTRDReferences(0),
	126	fTOFReferences(0)
	127	{
7ff5e4cf	128	//
	129	// copy constructor
	130	//
5c7ef659	131	fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
	132	fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
	133	fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
	134	fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
	135	}
	136
	137
	138	AliMCInfo::~AliMCInfo()
	139	{
7ff5e4cf	140	//
	141	// Destructor of the class
	142	//
5c7ef659	143	if (fTPCReferences) {
	144	delete fTPCReferences;
	145	}
	146	if (fITSReferences){
	147	delete fITSReferences;
	148	}
	149	if (fTRDReferences){
	150	delete fTRDReferences;
	151	}
	152	if (fTOFReferences){
	153	delete fTOFReferences;
	154	}
	155
	156	}
	157
	158
	159
	160	void AliMCInfo::Update()
	161	{
7ff5e4cf	162	//
	163	// Update MC info
	164	// Calculates some derived variables
5c7ef659	165	//
	166	//
	167	fMCtracks =1;
	168	if (!fTPCReferences) {
	169	fNTPCRef =0;
	170	return;
	171	}
	172	Float_t direction=1;
	173	//Float_t rlast=0;
	174	fNTPCRef = fTPCReferences->GetEntriesFast();
	175	fNITSRef = fITSReferences->GetEntriesFast();
	176	fNTRDRef = fTRDReferences->GetEntriesFast();
	177	fNTOFRef = fTOFReferences->GetEntriesFast();
	178
	179	for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
	180	AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
	181	//Float_t r = (ref->X()ref->X()+ref->Y()ref->Y());
	182	Float_t newdirection = ref->X()ref->Px()+ref->Y()ref->Py(); //inside or outside
	183	if (iref==0) direction = newdirection;
	184	if ( newdirection*direction<0){
	185	//changed direction
	186	direction = newdirection;
	187	fMCtracks+=1;
	188	}
	189	//rlast=r;
	190	}
	191	//
	192	// decay info
	193	fTPCdecay=kFALSE;
	194	if (fTRdecay.GetTrack()>0){
	195	fDecayCoord[0] = fTRdecay.X();
	196	fDecayCoord[1] = fTRdecay.Y();
	197	fDecayCoord[2] = fTRdecay.Z();
	198	if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
	199	fTPCdecay=kTRUE;
	200	}
	201	else{
	202	fDecayCoord[0] = 0;
	203	fDecayCoord[1] = 0;
	204	fDecayCoord[2] = 0;
	205	}
	206	}
	207	//
	208	//
	209	//digits information update
	210	fRowsWithDigits = fTPCRow.RowsOn();
	211	fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
	212	fRowsTrackLength = fTPCRow.Last() - fTPCRow.First();
	213	//
	214	//
	215	// calculate primary ionization per cm
	216	if (fParticle.GetPDG()){
	217	fMass = fParticle.GetMass();
	218	fCharge = fParticle.GetPDG()->Charge();
	219	if (fTPCReferences->GetEntriesFast()>0){
	220	fTrackRef = ((AliTrackReference)fTPCReferences->At(0));
	221	}
	222	if (fMass>0){
	223	Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
	224	fTrackRef.Py()*fTrackRef.Py()+
	225	fTrackRef.Pz()*fTrackRef.Pz());
	226	if (p>0.001){
	227	Float_t betagama = p /fMass;
	228	fPrim = TPCBetheBloch(betagama);
229	}else fPrim=0;
230	}
231	}else{
232	fMass =0;
233	fPrim =0;
234	}
235	}
236
237
238
239
240	////////////////////////////////////////////////////////////////////////
241	AliTPCdigitRow::AliTPCdigitRow()
242	{
243	Reset();
244	}
245	////////////////////////////////////////////////////////////////////////
246	AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
247	{
7ff5e4cf	248	for (Int_t i = 0; i<32; i++) fDig[i] = digOld.fDig[i];
5c7ef659	249	return (*this);
	250	}
	251	////////////////////////////////////////////////////////////////////////
	252	void AliTPCdigitRow::SetRow(Int_t row)
	253	{
7ff5e4cf	254	//
	255	// set bit mask for given row
	256	//
	257	if (row >= 8*32) {
	258	// cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
5c7ef659	259	return;
	260	}
	261	Int_t iC = row/8;
	262	Int_t iB = row%8;
	263	SETBIT(fDig[iC],iB);
	264	}
	265
	266	////////////////////////////////////////////////////////////////////////
	267	Bool_t AliTPCdigitRow::TestRow(Int_t row) const
	268	{
	269	//
	270	// return kTRUE if row is on
	271	//
	272	Int_t iC = row/8;
	273	Int_t iB = row%8;
	274	return TESTBIT(fDig[iC],iB);
	275	}
	276	////////////////////////////////////////////////////////////////////////
	277	Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
	278	{
	279	//
	280	// returns number of rows with a digit
	281	// count only rows less equal row number upto
	282	//
	283	Int_t total = 0;
7ff5e4cf	284	for (Int_t i = 0; i<32; i++) {
5c7ef659	285	for (Int_t j = 0; j < 8; j++) {
	286	if (i*8+j > upto) return total;
	287	if (TESTBIT(fDig[i],j)) total++;
	288	}
	289	}
	290	return total;
	291	}
	292	////////////////////////////////////////////////////////////////////////
	293	void AliTPCdigitRow::Reset()
	294	{
	295	//
	296	// resets all rows to zero
	297	//
7ff5e4cf	298	for (Int_t i = 0; i<32; i++) {
5c7ef659	299	fDig[i] <<= 8;
	300	}
	301	}
	302	////////////////////////////////////////////////////////////////////////
	303	Int_t AliTPCdigitRow::Last() const
	304	{
	305	//
	306	// returns the last row number with a digit
	307	// returns -1 if now digits
	308	//
7ff5e4cf	309	for (Int_t i = 32-1; i>=0; i--) {
5c7ef659	310	for (Int_t j = 7; j >= 0; j--) {
	311	if TESTBIT(fDig[i],j) return i*8+j;
	312	}
	313	}
	314	return -1;
	315	}
	316	////////////////////////////////////////////////////////////////////////
	317	Int_t AliTPCdigitRow::First() const
	318	{
	319	//
	320	// returns the first row number with a digit
	321	// returns -1 if now digits
	322	//
7ff5e4cf	323	for (Int_t i = 0; i<32; i++) {
5c7ef659	324	for (Int_t j = 0; j < 8; j++) {
	325	if (TESTBIT(fDig[i],j)) return i*8+j;
	326	}
	327	}
	328	return -1;
	329	}
	330
	331
	332	//_____________________________________________________________________________
	333	Float_t AliMCInfo::TPCBetheBloch(Float_t bg)
	334	{
	335	//
	336	// Bethe-Bloch energy loss formula
	337	//
	338	const Double_t kp1=0.76176e-1;
	339	const Double_t kp2=10.632;
	340	const Double_t kp3=0.13279e-4;
	341	const Double_t kp4=1.8631;
	342	const Double_t kp5=1.9479;
	343
	344	Double_t dbg = (Double_t) bg;
	345
	346	Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
	347
	348	Double_t aa = TMath::Power(beta,kp4);
	349	Double_t bb = TMath::Power(1./dbg,kp5);
	350
	351	bb=TMath::Log(kp3+bb);
	352
	353	return ((Float_t)((kp2-aa-bb)*kp1/aa));
	354	}
	355