[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

*/

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdio.h>
#include <string.h>
//ROOT includes
#include "TROOT.h"
#include "Rtypes.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TCut.h"
#include "TString.h"
#include "TStopwatch.h"
#include "TParticle.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TGeometry.h"
#include "TPolyLine3D.h"

//ALIROOT includes
#include "AliRun.h"
#include "AliStack.h"
#include "AliSimDigits.h"
#include "AliTPCParam.h"
#include "AliTPC.h"
#include "AliTPCLoader.h"
#include "AliDetector.h"
#include "AliTrackReference.h"
#include "AliTPCParamSR.h"
#include "AliTracker.h"
#include "AliMagF.h"
#include "AliHelix.h"
#include "AliTrackPointArray.h"

#endif
#include "AliMCInfo.h" 
//
// 

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)
{
  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)
{
  fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
  fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
  fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
  fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
}


AliMCInfo::~AliMCInfo()
{
  if (fTPCReferences) {
    delete fTPCReferences;
  }
  if (fITSReferences){
    delete fITSReferences;
  }
  if (fTRDReferences){    
    delete fTRDReferences;  
  }
  if (fTOFReferences){    
    delete fTOFReferences;  
  }

}


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