[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGenerator.cxx

#include <iostream>
#include "AliToyMCEventGenerator.h"
#include <AliTPCROC.h>
#include <AliTrackPointArray.h>
#include <TDatabasePDG.h>
#include <TRandom.h>
#include <TH2F.h>
#include <AliTrackerBase.h>
#include <AliCDBManager.h>
#include <AliTPCParam.h>
#include <AliGeomManager.h>
#include <AliTPCcalibDB.h>
#include <TGeoGlobalMagField.h>
#include <AliTPCclusterMI.h>
ClassImp(AliToyMCEventGenerator);


AliToyMCEventGenerator::AliToyMCEventGenerator()
  :TObject()
  ,fTPCParam(0x0)
{
  //TODO: Add method to set custom space charge file
  fSpaceCharge = new AliTPCSpaceCharge3D();
  fSpaceCharge->SetSCDataFileName("SC_NeCO2_eps5_50kHz.root");
  fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
  //fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
  fSpaceCharge->InitSpaceCharge3DDistortion();
  fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
  fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz"); 
  const char* ocdb="local://$ALICE_ROOT/OCDB/";
  AliCDBManager::Instance()->SetDefaultStorage(ocdb);
  AliCDBManager::Instance()->SetRun(0);   
  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  AliGeomManager::LoadGeometry("");
  fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  fTPCParam->ReadGeoMatrices();
}
//________________________________________________________________
AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
  :TObject(gen)
  ,fSpaceCharge(gen.fSpaceCharge)
{
  //
}
//________________________________________________________________
AliToyMCEventGenerator::~AliToyMCEventGenerator() 
{
  delete fSpaceCharge;
  delete fTPCParam;
}
//________________________________________________________________
Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0) {

  if(!fTPCParam) {
    fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
    fTPCParam->ReadGeoMatrices();
    std::cout << "init tpc param" << std::endl;
  }
  const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36);
  // const Double_t kRTPC0  =AliTPCROC::Instance()->GetPadRowRadii(0,0);
  // const Double_t kRTPC1  =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1);
  const Double_t kMaxSnp = 0.85;  
  const Double_t kSigmaY=0.1;
  const Double_t kSigmaZ=0.1;
  const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
  const Double_t kMaxZ0=fTPCParam->GetZLength();
  const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
    
  const Double_t iFCRadius=  83.5; //radius constants found in AliTPCCorrection.cxx 
  const Double_t oFCRadius= 254.5;
  
  AliToyMCTrack track(trackIn);
  const Int_t nMaxPoints = kNRows+50;
  AliTrackPointArray pointArray0(nMaxPoints);
  AliTrackPointArray pointArray1(nMaxPoints);
  Double_t xyz[3];

  if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) return 0;

  Int_t npoints=0;
  Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame


  //Simulate track from inner field cage radius to outer and save points along trajectory
  Int_t nMinPoints = 40;  
  for (Double_t radius=iFCRadius; radius<oFCRadius; radius++){

    if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
    track.GetXYZ(xyz);
    xyz[0]+=gRandom->Gaus(0,0.000005);
    xyz[1]+=gRandom->Gaus(0,0.000005);
    xyz[2]+=gRandom->Gaus(0,0.000005);

    if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
    if (TMath::Abs(track.GetX())<iFCRadius) continue;
    if (TMath::Abs(track.GetX())>oFCRadius) continue;

    AliTrackPoint pIn0;                               // space point          
    AliTrackPoint pIn1;
    Int_t sector= (xyz[2]>0)? 0:18;
    pointArray0.GetPoint(pIn0,npoints);
    pointArray1.GetPoint(pIn1,npoints);
    Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
    Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
    fSpaceCharge->DistortPoint(distPoint, sector);
    pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
    pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
    track.Rotate(alpha);
    AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
    AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
    prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
    prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
    pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
    pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
    pointArray0.AddPoint(npoints, &pIn0);      
    pointArray1.AddPoint(npoints, &pIn1);
    npoints++;
    if (npoints>=nMaxPoints) break;
  }

  if (npoints<nMinPoints) return 0;

  //save space points and make clusters of distorted points
  Int_t lastRow = 0;
  Int_t pntsInCurrentRow = 0;
  Double_t xt=0, yt=0, zt=0;
  for(Int_t iPoint = 0; iPoint<npoints; iPoint++){
    
    AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());

    const Float_t *xArr = pointArray0.GetX();
    const Float_t *yArr = pointArray0.GetY();
    const Float_t *zArr = pointArray0.GetZ();

    const Float_t *xArrDist = pointArray1.GetX();
    const Float_t *yArrDist = pointArray1.GetY();
    const Float_t *zArrDist = pointArray1.GetZ();

    tempCl->SetX(xArr[iPoint]);
    tempCl->SetY(yArr[iPoint]);
    tempCl->SetZ(zArr[iPoint]);
    
    Float_t xyz[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]};
    Int_t index[3] = {0};

    Float_t padRow = fTPCParam->GetPadRow(xyz,index);
    Int_t sector = index[1];

    if(padRow < 0 || (sector < 36 && padRow>62) || (sector > 35 &&padRow > 95) ) continue; //outside sensitive area

    
    if(lastRow!=padRow) {
      

      //make and save distorted cluster
      if(pntsInCurrentRow>0){
	xt/=pntsInCurrentRow;
	yt/=pntsInCurrentRow;
	zt/=pntsInCurrentRow;
	AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI());
	tempDistCl->SetX(xt);
	tempDistCl->SetY(yt);
	tempDistCl->SetZ(zt);
	Float_t clxyz[3] = {xt,yt,zt};
	Int_t clindex[3] = {0};
	Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex);
	Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow);
	Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx
	tempDistCl->SetPad(pad);
	tempDistCl->SetRow(clRow);
	tempDistCl->SetDetector(clindex[1]);
	tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0  + drift time from dist z coordinate to pad plane
      }
      xt = 0;
      yt = 0;
      zt = 0;
      pntsInCurrentRow = 0;
    }

    xt+=xArrDist[iPoint];
    yt+=yArrDist[iPoint];
    zt+=zArrDist[iPoint];
    pntsInCurrentRow++;
    lastRow = padRow;
  }
 
  return 1;

  
}
Commit	Line	Data
526ddf0e	1	#include <iostream>
de0014b7	2	#include "AliToyMCEventGenerator.h"
526ddf0e	3	#include <AliTPCROC.h>
	4	#include <AliTrackPointArray.h>
	5	#include <TDatabasePDG.h>
	6	#include <TRandom.h>
	7	#include <TH2F.h>
	8	#include <AliTrackerBase.h>
	9	#include <AliCDBManager.h>
	10	#include <AliTPCParam.h>
	11	#include <AliGeomManager.h>
	12	#include <AliTPCcalibDB.h>
	13	#include <TGeoGlobalMagField.h>
	14	#include <AliTPCclusterMI.h>
de0014b7	15	ClassImp(AliToyMCEventGenerator);
526ddf0e	16
526ddf0e	17
de0014b7	18	AliToyMCEventGenerator::AliToyMCEventGenerator()
526ddf0e	19	:TObject()
	20	,fTPCParam(0x0)
	21	{
	22	//TODO: Add method to set custom space charge file
	23	fSpaceCharge = new AliTPCSpaceCharge3D();
	24	fSpaceCharge->SetSCDataFileName("SC_NeCO2_eps5_50kHz.root");
	25	fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
	26	//fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
	27	fSpaceCharge->InitSpaceCharge3DDistortion();
	28	fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
	29	fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
	30	const char* ocdb="local://$ALICE_ROOT/OCDB/";
	31	AliCDBManager::Instance()->SetDefaultStorage(ocdb);
	32	AliCDBManager::Instance()->SetRun(0);
	33	TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
	34	AliGeomManager::LoadGeometry("");
	35	fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
	36	fTPCParam->ReadGeoMatrices();
	37	}
	38	//________________________________________________________________
de0014b7	39	AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
526ddf0e	40	:TObject(gen)
	41	,fSpaceCharge(gen.fSpaceCharge)
	42	{
	43	//
	44	}
	45	//________________________________________________________________
de0014b7	46	AliToyMCEventGenerator::~AliToyMCEventGenerator()
526ddf0e	47	{
	48	delete fSpaceCharge;
	49	delete fTPCParam;
	50	}
	51	//________________________________________________________________
de0014b7	52	Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0) {
526ddf0e	53
	54	if(!fTPCParam) {
	55	fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
	56	fTPCParam->ReadGeoMatrices();
	57	std::cout << "init tpc param" << std::endl;
	58	}
	59	const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36);
	60	// const Double_t kRTPC0 =AliTPCROC::Instance()->GetPadRowRadii(0,0);
	61	// const Double_t kRTPC1 =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1);
	62	const Double_t kMaxSnp = 0.85;
	63	const Double_t kSigmaY=0.1;
	64	const Double_t kSigmaZ=0.1;
	65	const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
	66	const Double_t kMaxZ0=fTPCParam->GetZLength();
	67	const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
	68
	69	const Double_t iFCRadius= 83.5; //radius constants found in AliTPCCorrection.cxx
	70	const Double_t oFCRadius= 254.5;
	71
de0014b7	72	AliToyMCTrack track(trackIn);
526ddf0e	73	const Int_t nMaxPoints = kNRows+50;
	74	AliTrackPointArray pointArray0(nMaxPoints);
	75	AliTrackPointArray pointArray1(nMaxPoints);
	76	Double_t xyz[3];
	77
	78	if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) return 0;
	79
	80	Int_t npoints=0;
	81	Float_t covPoint[6]={0,0,0, kSigmaYkSigmaY,0,kSigmaZkSigmaZ}; //covariance at the local frame
	82
	83
	84	//Simulate track from inner field cage radius to outer and save points along trajectory
	85	Int_t nMinPoints = 40;
	86	for (Double_t radius=iFCRadius; radius<oFCRadius; radius++){
	87
	88	if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
	89	track.GetXYZ(xyz);
	90	xyz[0]+=gRandom->Gaus(0,0.000005);
	91	xyz[1]+=gRandom->Gaus(0,0.000005);
	92	xyz[2]+=gRandom->Gaus(0,0.000005);
	93
	94	if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
	95	if (TMath::Abs(track.GetX())<iFCRadius) continue;
	96	if (TMath::Abs(track.GetX())>oFCRadius) continue;
	97
	98	AliTrackPoint pIn0; // space point
	99	AliTrackPoint pIn1;
	100	Int_t sector= (xyz[2]>0)? 0:18;
	101	pointArray0.GetPoint(pIn0,npoints);
	102	pointArray1.GetPoint(pIn1,npoints);
	103	Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
	104	Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
	105	fSpaceCharge->DistortPoint(distPoint, sector);
	106	pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
	107	pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
	108	track.Rotate(alpha);
	109	AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
	110	AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
	111	prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
	112	prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
	113	pIn0=prot0.Rotate(-alpha); // rotate back to global frame
	114	pIn1=prot1.Rotate(-alpha); // rotate back to global frame
	115	pointArray0.AddPoint(npoints, &pIn0);
	116	pointArray1.AddPoint(npoints, &pIn1);
	117	npoints++;
	118	if (npoints>=nMaxPoints) break;
	119	}
	120
	121	if (npoints<nMinPoints) return 0;
	122
	123	//save space points and make clusters of distorted points
	124	Int_t lastRow = 0;
	125	Int_t pntsInCurrentRow = 0;
	126	Double_t xt=0, yt=0, zt=0;
	127	for(Int_t iPoint = 0; iPoint<npoints; iPoint++){
	128
	129	AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());
	130
	131	const Float_t *xArr = pointArray0.GetX();
	132	const Float_t *yArr = pointArray0.GetY();
	133	const Float_t *zArr = pointArray0.GetZ();
	134
	135	const Float_t *xArrDist = pointArray1.GetX();
	136	const Float_t *yArrDist = pointArray1.GetY();
137	const Float_t *zArrDist = pointArray1.GetZ();
138
139	tempCl->SetX(xArr[iPoint]);
140	tempCl->SetY(yArr[iPoint]);
141	tempCl->SetZ(zArr[iPoint]);
142
143	Float_t xyz[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]};
144	Int_t index[3] = {0};
145
146	Float_t padRow = fTPCParam->GetPadRow(xyz,index);
147	Int_t sector = index[1];
148
149	if(padRow < 0 \|\| (sector < 36 && padRow>62) \|\| (sector > 35 &&padRow > 95) ) continue; //outside sensitive area
150
151
152	if(lastRow!=padRow) {
153
154
155	//make and save distorted cluster
156	if(pntsInCurrentRow>0){
157	xt/=pntsInCurrentRow;
158	yt/=pntsInCurrentRow;
159	zt/=pntsInCurrentRow;
160	AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI());
161	tempDistCl->SetX(xt);
162	tempDistCl->SetY(yt);
163	tempDistCl->SetZ(zt);
164	Float_t clxyz[3] = {xt,yt,zt};
165	Int_t clindex[3] = {0};
166	Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex);
167	Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow);
168	Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx
169	tempDistCl->SetPad(pad);
170	tempDistCl->SetRow(clRow);
171	tempDistCl->SetDetector(clindex[1]);
172	tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0 + drift time from dist z coordinate to pad plane
173	}
174	xt = 0;
175	yt = 0;
176	zt = 0;
177	pntsInCurrentRow = 0;
178	}
179
180	xt+=xArrDist[iPoint];
181	yt+=yArrDist[iPoint];
182	zt+=zArrDist[iPoint];
183	pntsInCurrentRow++;
184	lastRow = padRow;
185	}
186
187	return 1;
188
189
190
191	}