35#include "GPUParam.inc"
41#include "../utils/qconfig.h"
50int32_t genEvents::GetSector(
double GlobalPhi)
52 double phi = GlobalPhi;
61 return (int32_t)(
phi / mSectorDAngle);
64int32_t genEvents::GetDSector(
double LocalPhi) {
return GetSector(LocalPhi + mSectorAngleOffset); }
66double genEvents::GetSectorAngle(int32_t iSector) {
return mSectorAngleOffset + iSector * mSectorDAngle; }
70 double phi = atan2(t.GetY(), t.GetX());
72 int32_t dSector = GetDSector(phi);
78 double dAlpha = dSector * mSectorDAngle;
90double genEvents::GetGaus(
double sigma)
94 x = gRandom->Gaus(0., sigma);
95 if (fabs(
x) <= 3.5 * sigma) {
104 const char*
rows[3] = {
"0-63",
"128-159",
"64-127"};
105 for (int32_t
i = 0;
i < 3;
i++) {
106 for (int32_t
j = 0;
j < 2;
j++) {
107 char name[1024], title[1024];
109 snprintf(
name, 1024,
"clError%s%d", (
j == 0 ?
"Y" :
"Z"),
i);
111 snprintf(title, 1024,
"Cluster %s Error for rows %s", (
j == 0 ?
"Y" :
"Z"),
rows[
i]);
113 mClusterError[
i][
j] =
new TH1F(
name, title, 1000, 0., .7);
114 mClusterError[
i][
j]->GetXaxis()->SetTitle(
"Cluster Error [cm]");
121 TFile* tout =
new TFile(
"generator.root",
"RECREATE");
122 TCanvas*
c =
new TCanvas(
"ClusterErrors",
"Cluste rErrors", 0, 0, 700, 700. * 2. / 3.);
125 for (int32_t
j = 0;
j < 2;
j++) {
126 for (int32_t
i = 0;
i < 3;
i++) {
136 mClusterError[k][
j]->Write();
139 mClusterError[k][
j]->Draw();
144 c->Print(
"plots/clusterErrors.pdf");
154 mRec->ClearIOPointers();
155 static int32_t iEvent = -1;
162 cout <<
"NTracks " << nTracks << endl;
163 std::vector<GPUTPCMCInfo> mcInfo(nTracks);
164 memset(mcInfo.data(), 0, nTracks *
sizeof(mcInfo[0]));
171 prop.SetToyMCEventsFlag(kTRUE);
173 prop.SetPolynomialField(&merger.Param().polynomialField);
178 std::vector<GenCluster> vClusters;
179 int32_t clusterId = 0;
183 for (int32_t itr = 0; itr < nTracks; itr++) {
187 mcInfo[itr].pid = -100;
189 mcInfo[itr].prim = 1;
190 mcInfo[itr].primDaughters = 0;
199 double dphi = mTwoPi / nTracks;
200 double phi = mSectorAngleOffset + dphi * itr;
201 double eta = gRandom->Uniform(-1.5, 1.5);
203 double theta = 2 * std::atan(1. /
exp(eta));
204 double lambda = theta - M_PI / 2;
206 double pt = .08 * std::pow(10, gRandom->Uniform(0, 2.2));
209 int32_t iSector = GetSector(phi);
210 phi =
phi - GetSectorAngle(iSector);
213 double x0 = cosf(phi);
214 double y0 = sinf(phi);
215 double z0 = tanf(lambda);
216 t.Set(
x0,
y0, z0, pt *
x0, pt *
y0, pt * z0, q);
218 if (RecalculateSector(t, iSector) != 0) {
219 std::cout <<
"Initial sector wrong!!!" << std::endl;
225 float xRow = GPUTPCGeometry::Row2X(iRow);
228 for (int32_t itry = 0; itry < 1; itry++) {
230 prop.GetBxByBz(GetSectorAngle(iSector), t.GetX(), t.GetY(), t.GetZ(),
B);
232 err = t.PropagateToXBxByBz(xRow,
B[0],
B[1],
B[2], dLp);
234 std::cout <<
"Can not propagate to x = " << xRow << std::endl;
238 if (fabsf(t.GetZ()) >= GPUTPCGeometry::TPCLength()) {
239 std::cout <<
"Can not propagate to x = " << xRow <<
": Z outside the volume" << std::endl;
245 int32_t isNewSector = RecalculateSector(t, iSector);
249 std::cout <<
"track " << itr <<
": new sector " << iSector <<
" at row " << iRow << std::endl;
262 tg.Rotate(-GetSectorAngle(iSector));
265 mcInfo[itr].charge = 3 * q;
266 mcInfo[itr].x = tg.GetX();
267 mcInfo[itr].y = tg.GetY();
268 mcInfo[itr].z = tg.GetZ();
269 mcInfo[itr].pX = tg.GetPx();
270 mcInfo[itr].pY = tg.GetPy();
271 mcInfo[itr].pZ = tg.GetPz();
279 const int32_t rowType = iRow < 64 ? 0 : iRow < 128 ? 2 : 1;
281 param.GetClusterErrors2(iSector, rowType, t.GetZ(), t.GetSinPhi(), t.GetDzDs(), -1.f, 0.f, 0.f, sigmaY, sigmaZ);
282 sigmaY = std::sqrt(sigmaY);
283 sigmaZ = std::sqrt(sigmaZ);
284 mClusterError[rowType][0]->Fill(sigmaY);
285 mClusterError[rowType][1]->Fill(sigmaZ);
289 c.sector = (t.GetZ() >= 0.) ? iSector : iSector + 18;
293 c.y = t.GetY() + GetGaus(sigmaY);
294 c.z = t.GetZ() + GetGaus(sigmaZ);
296 vClusters.push_back(
c);
300 std::vector<AliHLTTPCClusterMCLabel> labels;
306 int32_t nNumberOfHits = 0;
307 for (uint32_t
i = 0;
i < vClusters.size();
i++) {
308 if (vClusters[
i].sector == iSector) {
313 mRec->mIOPtrs.nClusterData[iSector] = nNumberOfHits;
318 for (uint32_t
i = 0;
i < vClusters.size();
i++) {
319 GenCluster&
c = vClusters[
i];
320 if (
c.sector == iSector) {
329 for (int32_t
j = 0;
j < 3;
j++) {
335 labels.push_back(clusterLabel);
338 mRec->mIOPtrs.clusterData[iSector] =
clusters;
343 mRec->mIOPtrs.nMCLabelsTPC = labels.size();
344 mRec->mIOPtrs.mcLabelsTPC = labels.data();
346 mRec->mIOPtrs.nMCInfosTPC = mcInfo.size();
347 mRec->mIOPtrs.mcInfosTPC = mcInfo.data();
349 mRec->mIOPtrs.mcInfosTPCCol = &mcColInfo;
350 mRec->mIOPtrs.nMCInfosTPCCol = 1;
363 mkdir(dirname, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
366 gen->InitEventGenerator();
373 gen->FinishEventGenerator();
default_random_engine gen(dev())
uint64_t exp(uint64_t base, uint8_t exp) noexcept
#define GPUCA_EVDUMP_FILE
static constexpr int32_t NSECTORS
const GPUParam & GetParam() const
void DumpSettings(const char *dir="")
float const GPUParam float float float float float int32_t int16_t int8_t float float float charge
static void RunEventGenerator(GPUChainTracking *rec)
void InitEventGenerator()
int32_t GenerateEvent(const GPUParam §orParam, char *filename)
void FinishEventGenerator()
GLuint const GLchar * name
GLuint GLfloat GLfloat y0
GPUSettingsStandalone configStandalone
Defining DataPointCompositeObject explicitly as copiable.
AliHLTTPCClusterMCWeight fClusterID[3]
std::vector< Cluster > clusters
std::vector< ReadoutWindowData > rows