35#include "GPUParam.inc"
42#include "../utils/qconfig.h"
51int32_t genEvents::GetSector(
double GlobalPhi)
53 double phi = GlobalPhi;
62 return (int32_t)(
phi / mSectorDAngle);
65int32_t genEvents::GetDSector(
double LocalPhi) {
return GetSector(LocalPhi + mSectorAngleOffset); }
67double genEvents::GetSectorAngle(int32_t iSector) {
return mSectorAngleOffset + iSector * mSectorDAngle; }
71 double phi = atan2(t.GetY(), t.GetX());
73 int32_t dSector = GetDSector(phi);
79 double dAlpha = dSector * mSectorDAngle;
91double genEvents::GetGaus(
double sigma)
95 x = gRandom->Gaus(0., sigma);
96 if (fabs(
x) <= 3.5 * sigma) {
105 const char*
rows[3] = {
"0-63",
"128-159",
"64-127"};
106 for (int32_t
i = 0;
i < 3;
i++) {
107 for (int32_t
j = 0;
j < 2;
j++) {
108 char name[1024], title[1024];
110 snprintf(
name, 1024,
"clError%s%d", (
j == 0 ?
"Y" :
"Z"),
i);
112 snprintf(title, 1024,
"Cluster %s Error for rows %s", (
j == 0 ?
"Y" :
"Z"),
rows[
i]);
114 mClusterError[
i][
j] =
new TH1F(
name, title, 1000, 0., .7);
115 mClusterError[
i][
j]->GetXaxis()->SetTitle(
"Cluster Error [cm]");
122 TFile* tout =
new TFile(
"generator.root",
"RECREATE");
123 TCanvas*
c =
new TCanvas(
"ClusterErrors",
"Cluste rErrors", 0, 0, 700, 700. * 2. / 3.);
126 for (int32_t
j = 0;
j < 2;
j++) {
127 for (int32_t
i = 0;
i < 3;
i++) {
137 mClusterError[k][
j]->Write();
140 mClusterError[k][
j]->Draw();
145 c->Print(
"plots/clusterErrors.pdf");
155 mRec->ClearIOPointers();
156 static int32_t iEvent = -1;
163 cout <<
"NTracks " << nTracks << endl;
164 std::vector<GPUTPCMCInfo> mcInfo(nTracks);
165 memset(mcInfo.data(), 0, nTracks *
sizeof(mcInfo[0]));
172 prop.SetToyMCEventsFlag(kTRUE);
174 prop.SetPolynomialField(&merger.Param().polynomialField);
179 std::vector<GenCluster> vClusters;
180 int32_t clusterId = 0;
184 for (int32_t itr = 0; itr < nTracks; itr++) {
188 mcInfo[itr].pid = -100;
190 mcInfo[itr].prim = 1;
191 mcInfo[itr].primDaughters = 0;
200 double dphi = mTwoPi / nTracks;
201 double phi = mSectorAngleOffset + dphi * itr;
202 double eta = gRandom->Uniform(-1.5, 1.5);
204 double theta = 2 * std::atan(1. /
exp(eta));
205 double lambda = theta - M_PI / 2;
207 double pt = .08 * std::pow(10, gRandom->Uniform(0, 2.2));
210 int32_t iSector = GetSector(phi);
211 phi =
phi - GetSectorAngle(iSector);
214 double x0 = cosf(phi);
215 double y0 = sinf(phi);
216 double z0 = tanf(lambda);
217 t.Set(
x0,
y0, z0, pt *
x0, pt *
y0, pt * z0, q);
219 if (RecalculateSector(t, iSector) != 0) {
220 std::cout <<
"Initial sector wrong!!!" << std::endl;
226 float xRow = GPUTPCGeometry::Row2X(iRow);
229 for (int32_t itry = 0; itry < 1; itry++) {
231 prop.GetBxByBz(GetSectorAngle(iSector), t.GetX(), t.GetY(), t.GetZ(),
B);
233 err = t.PropagateToXBxByBz(xRow,
B[0],
B[1],
B[2], dLp);
235 std::cout <<
"Can not propagate to x = " << xRow << std::endl;
239 if (fabsf(t.GetZ()) >= GPUTPCGeometry::TPCLength()) {
240 std::cout <<
"Can not propagate to x = " << xRow <<
": Z outside the volume" << std::endl;
246 int32_t isNewSector = RecalculateSector(t, iSector);
250 std::cout <<
"track " << itr <<
": new sector " << iSector <<
" at row " << iRow << std::endl;
263 tg.Rotate(-GetSectorAngle(iSector));
266 mcInfo[itr].charge = 3 * q;
267 mcInfo[itr].x = tg.GetX();
268 mcInfo[itr].y = tg.GetY();
269 mcInfo[itr].z = tg.GetZ();
270 mcInfo[itr].pX = tg.GetPx();
271 mcInfo[itr].pY = tg.GetPy();
272 mcInfo[itr].pZ = tg.GetPz();
280 const int32_t rowType = iRow < 64 ? 0 : iRow < 128 ? 2 : 1;
282 param.GetClusterErrors2(iSector, rowType, t.GetZ(), t.GetSinPhi(), t.GetDzDs(), -1.f, 0.f, 0.f, sigmaY, sigmaZ);
283 sigmaY = std::sqrt(sigmaY);
284 sigmaZ = std::sqrt(sigmaZ);
285 mClusterError[rowType][0]->Fill(sigmaY);
286 mClusterError[rowType][1]->Fill(sigmaZ);
290 c.sector = (t.GetZ() >= 0.) ? iSector : iSector + 18;
294 c.y = t.GetY() + GetGaus(sigmaY);
295 c.z = t.GetZ() + GetGaus(sigmaZ);
297 vClusters.push_back(
c);
301 std::vector<AliHLTTPCClusterMCLabel> labels;
307 int32_t nNumberOfHits = 0;
308 for (uint32_t
i = 0;
i < vClusters.size();
i++) {
309 if (vClusters[
i].sector == iSector) {
314 mRec->mIOPtrs.nClusterData[iSector] = nNumberOfHits;
319 for (uint32_t
i = 0;
i < vClusters.size();
i++) {
320 GenCluster&
c = vClusters[
i];
321 if (
c.sector == iSector) {
330 for (int32_t
j = 0;
j < 3;
j++) {
336 labels.push_back(clusterLabel);
339 mRec->mIOPtrs.clusterData[iSector] =
clusters;
344 mRec->mIOPtrs.nMCLabelsTPC = labels.size();
345 mRec->mIOPtrs.mcLabelsTPC = labels.data();
347 mRec->mIOPtrs.nMCInfosTPC = mcInfo.size();
348 mRec->mIOPtrs.mcInfosTPC = mcInfo.data();
350 mRec->mIOPtrs.mcInfosTPCCol = &mcColInfo;
351 mRec->mIOPtrs.nMCInfosTPCCol = 1;
364 mkdir(dirname, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
367 gen->InitEventGenerator();
374 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