o add setters for space charge files
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGenerator.cxx
CommitLineData
526ddf0e 1#include <iostream>
a1a695e5 2
526ddf0e 3#include <TDatabasePDG.h>
4#include <TRandom.h>
5#include <TH2F.h>
a1a695e5 6#include <TGeoGlobalMagField.h>
32438f4e 7#include <TSpline.h>
cd8ed0ac 8#include <TObjString.h>
a1a695e5 9
b0c26b55 10#include <AliLog.h>
a1a695e5 11#include <AliTPCROC.h>
12#include <AliTrackPointArray.h>
526ddf0e 13#include <AliTrackerBase.h>
14#include <AliCDBManager.h>
15#include <AliTPCParam.h>
16#include <AliGeomManager.h>
17#include <AliTPCcalibDB.h>
526ddf0e 18#include <AliTPCclusterMI.h>
a1a695e5 19#include <AliTPCSpaceCharge3D.h>
32438f4e 20#include <AliTPCROC.h>
a1a695e5 21
22#include "AliToyMCEvent.h"
23#include "AliToyMCTrack.h"
24
25#include "AliToyMCEventGenerator.h"
26
de0014b7 27ClassImp(AliToyMCEventGenerator);
526ddf0e 28
29
de0014b7 30AliToyMCEventGenerator::AliToyMCEventGenerator()
526ddf0e 31 :TObject()
32 ,fTPCParam(0x0)
a1a695e5 33 ,fEvent(0x0)
34 ,fSpaceCharge(0x0)
cd8ed0ac 35 ,fSpaceChargeFile("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root")
a1a695e5 36 ,fOutputFileName("toyMC.root")
37 ,fOutFile(0x0)
38 ,fOutTree(0x0)
526ddf0e 39{
526ddf0e 40 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
41 fTPCParam->ReadGeoMatrices();
42}
43//________________________________________________________________
de0014b7 44AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
526ddf0e 45 :TObject(gen)
a1a695e5 46 ,fTPCParam(gen.fTPCParam)
47 ,fEvent(0x0)
526ddf0e 48 ,fSpaceCharge(gen.fSpaceCharge)
cd8ed0ac 49 ,fSpaceChargeFile(gen.fSpaceChargeFile)
a1a695e5 50 ,fOutputFileName(gen.fOutputFileName)
51 ,fOutFile(0x0)
52 ,fOutTree(0x0)
526ddf0e 53{
54 //
55}
56//________________________________________________________________
de0014b7 57AliToyMCEventGenerator::~AliToyMCEventGenerator()
526ddf0e 58{
59 delete fSpaceCharge;
526ddf0e 60}
32438f4e 61
526ddf0e 62//________________________________________________________________
a1a695e5 63Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
64{
65 //
66 //
67 //
68
526ddf0e 69 if(!fTPCParam) {
70 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
71 fTPCParam->ReadGeoMatrices();
526ddf0e 72 }
32438f4e 73
74 // make it big enough to hold all points
75 // store real number of generated points in the unique id
76 const Int_t nMaxPoints=3000;
77 AliTrackPointArray pointArray0(nMaxPoints); //undistorted
78 AliTrackPointArray pointArray1(nMaxPoints); //distorted
79
80 //Create space point of undistorted and distorted clusters along the propagated track trajectory
81 CreateSpacePoints(trackIn,pointArray0,pointArray1);
82 //Convert the space points into clusters in the local frame
83 //for undistorted and distorted clusters using the same function
84 ConvertTrackPointsToLocalClusters(pointArray0,trackIn,t0,0);
85 ConvertTrackPointsToLocalClusters(pointArray1,trackIn,t0,1);
526ddf0e 86
32438f4e 87 return 1;
88}
526ddf0e 89
32438f4e 90//________________________________________________________________
91void AliToyMCEventGenerator::CreateSpacePoints(AliToyMCTrack &trackIn,
92 AliTrackPointArray &arrUdist,
93 AliTrackPointArray &arrDist)
94{
95 //
96 // sample the track from the inner to the outer wall of the TPC
97 // a graph is filled in local coordinates for later
98 //
99
100 const Double_t kMaxSnp = 0.85;
101 const Double_t kMaxZ0 = fTPCParam->GetZLength();
102 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
103
104 const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
105 const Double_t oFCRadius = 254.5;
106
107 AliToyMCTrack track(trackIn);
108 //!!! TODO: make this adjustable perhaps
109 const Double_t stepSize=0.1;
110 Double_t xyz[3] = {0.,0.,0.};
111 Float_t xyzf[3] = {0.,0.,0.};
112
a1a695e5 113 //!!! when does the propagation not work, how often does it happen?
b0c26b55 114 if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) {
115 AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
32438f4e 116 return;
b0c26b55 117 }
32438f4e 118
526ddf0e 119 Int_t npoints=0;
b0c26b55 120
32438f4e 121 for (Double_t radius=iFCRadius; radius<oFCRadius; radius+=stepSize){
122
b0c26b55 123 //!!! changed from return 0 to continue -> Please check
32438f4e 124 if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,1,kTRUE,kMaxSnp)) {
125 AliError(Form("Propagation to %.2f failed\n",radius));
b0c26b55 126 continue;
127 }
526ddf0e 128 track.GetXYZ(xyz);
32438f4e 129
b0c26b55 130 //!!! Why is this smeared
526ddf0e 131 xyz[0]+=gRandom->Gaus(0,0.000005);
132 xyz[1]+=gRandom->Gaus(0,0.000005);
133 xyz[2]+=gRandom->Gaus(0,0.000005);
32438f4e 134
135 xyzf[0]=Float_t(xyz[0]);
136 xyzf[1]=Float_t(xyz[1]);
137 xyzf[2]=Float_t(xyz[2]);
138
526ddf0e 139 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
32438f4e 140 // if (TMath::Abs(track.GetX())<iFCRadius) continue;
141 // if (TMath::Abs(track.GetX())>oFCRadius) continue;
142
143
144 AliTrackPoint pUdist; // undistorted space point
145 AliTrackPoint pDist; // distorted space point
146 // Set undistorted point
147 SetPoint(xyzf,pUdist);
148 arrUdist.AddPoint(npoints, &pUdist);
149 Int_t sector=pUdist.GetVolumeID();
150 // abuse volume ID for the sector number
151 pUdist.SetVolumeID(sector);
152
153 // set distorted point
526ddf0e 154 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
155 fSpaceCharge->DistortPoint(distPoint, sector);
32438f4e 156 SetPoint(distPoint, pDist);
157 arrDist.AddPoint(npoints, &pDist);
158 // abuse volume ID for the sector number
159 sector=pDist.GetVolumeID();
160 pDist.SetVolumeID(sector);
161
162 ++npoints;
526ddf0e 163 }
32438f4e 164
165 arrUdist.SetUniqueID(npoints);
166 arrDist.SetUniqueID(npoints);
167}
526ddf0e 168
32438f4e 169//________________________________________________________________
170void AliToyMCEventGenerator::SetPoint(Float_t xyz[3], AliTrackPoint &point)
171{
172 //
173 // make AliTrackPoint out of AliTPCclusterMI
174 //
175
176 //covariance at the local frame
177 //assume 1mm distortion in y and z
178 Int_t i[3]={0,0,0};
179 const Double_t kSigmaY=0.1;
180 const Double_t kSigmaZ=0.1;
181 Float_t cov[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};
182
183 const Float_t alpha = -TMath::ATan2(xyz[1],xyz[0]);
184 const Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
185
186 Float_t newcov[6];
187 newcov[0] = cov[0]*cos*cos + 2*cov[1]*sin*cos + cov[3]*sin*sin;
188 newcov[1] = cov[1]*(cos*cos-sin*sin) + (cov[3]-cov[0])*sin*cos;
189 newcov[2] = cov[2]*cos + cov[4]*sin;
190 newcov[3] = cov[0]*sin*sin - 2*cov[1]*sin*cos + cov[3]*cos*cos;
191 newcov[4] = cov[4]*cos - cov[2]*sin;
192 newcov[5] = cov[5];
193
194 // voluem ID to add later ....
195 point.SetXYZ(xyz);
196 point.SetCov(newcov);
197 point.SetVolumeID(fTPCParam->Transform0to1(xyz,i));
b0c26b55 198
32438f4e 199 // TODO: Add sampled dE/dx (use SetCharge)
200}
b0c26b55 201
32438f4e 202//________________________________________________________________
203void AliToyMCEventGenerator::ConvertTrackPointsToLocalClusters(AliTrackPointArray &arrPoints,
204 AliToyMCTrack &tr, Double_t t0, Int_t type)
205{
206 //
207 //
208 //
526ddf0e 209
526ddf0e 210
32438f4e 211 const Int_t npoints=Int_t(arrPoints.GetUniqueID());
212 Int_t secOld=-1;
213
214 // create an array for graphs which are used for local interpolation
215 // we need a new graph if the sector changed since then the local frame will also change
216 TObjArray arrGraphsXY(72);
217 arrGraphsXY.SetOwner();
218 TObjArray arrGraphsXZ(72);
219 arrGraphsXZ.SetOwner();
220 AliTrackPoint p;
221 //create initial graph
222 TGraph *grXY=0x0;
223 TGraph *grXZ=0x0;
224 //row -> sector mapping
225 Int_t rowMap[159];
226 for (Int_t irow=0; irow<159; ++irow) rowMap[irow]=-1;
227
228 // 1. Step
229 // Make from the list of global space points graphs in the local frame
230 // one graph per sector is needed
231 for (Int_t ipoint=0; ipoint<npoints; ++ipoint){
232 arrPoints.GetPoint(p,ipoint);
233 Float_t xyz[3] = {p.GetX(),p.GetY(),p.GetZ()};
234 Int_t index[3] = {0,0,0};
235 Int_t row = fTPCParam->GetPadRow(xyz,index);
236 // rotate space point to local frame
237 // the angle is given by the VolumeID which was set in CrateSpacePoints
238 const Int_t sec=p.GetVolumeID();
239 if (row<0 || (sec<36 && row>62) || row>95 ) continue;
240 Double_t angle=((sec%18)*20.+10.)/TMath::RadToDeg();
241 AliTrackPoint pRot=p.Rotate(angle);
242 const Int_t secrow=row+(sec>35)*63;
243 if (rowMap[secrow]==-1) rowMap[secrow]=sec;
244 // check if we need a new graph (sector change)
245 if (secOld!=sec){
246 grXY=new TGraph;
247 grXZ=new TGraph;
248 arrGraphsXY.AddAt(grXY,sec);
249 arrGraphsXZ.AddAt(grXZ,sec);
250 }
b0c26b55 251
32438f4e 252 //add coordinates in local frame for later interpolation
253 grXY->SetPoint(grXY->GetN(), pRot.GetX(), pRot.GetY());
254 grXZ->SetPoint(grXZ->GetN(), pRot.GetX(), pRot.GetZ());
255 secOld=sec;
256 }
b0c26b55 257
32438f4e 258 // 2. Step
259 // create in the center of each row a space point by using the graph to interpolate
260 // the the center of the row. This is done in xy and xz
261 TSpline3 *splXY=0x0;
262 TSpline3 *splXZ=0x0;
263 AliTPCclusterMI tempCl;
264 secOld=-1;
265 for (Int_t irow=0; irow<159; ++irow ){
266 const Int_t sec = rowMap[irow];
267 if (sec==-1) continue;
268 const Int_t secrow = irow<63?irow:irow-63;
269 Double_t localX = fTPCParam->GetPadRowRadii(sec,secrow);
270 // get graph for the current row
271 if (sec!=secOld){
272 delete splXY;
273 splXY=0x0;
274 delete splXZ;
275 splXZ=0x0;
276
277 grXY=(TGraph*)arrGraphsXY.At(sec);
278 grXZ=(TGraph*)arrGraphsXZ.At(sec);
279 if (!grXY) continue;
b0c26b55 280
32438f4e 281 splXY=new TSpline3("splXY",grXY);
282 splXZ=new TSpline3("splXZ",grXZ);
283 }
284 secOld=sec;
b0c26b55 285
32438f4e 286 // check we are in an active area
287 if (splXY->FindX(localX)<1 || splXZ->FindX(localX)<1) continue;
b0c26b55 288
32438f4e 289 //get interpolated value at the center for the pad row
290 // using splines
291 const Double_t localY=splXY->Eval(localX/*,0x0,"S"*/);
292 const Double_t localZ=splXZ->Eval(localX/*,0x0,"S"*/);
293 Float_t xyz[3]={localX,localY,localZ};
526ddf0e 294
32438f4e 295 if (!SetupCluster(tempCl,xyz,sec,t0)) continue;
296
297 if (type==0) tr.AddSpacePoint(tempCl);
298 else tr.AddDistortedSpacePoint(tempCl);
299// printf("SetupCluster %3d: (%.2f, %.2f, %.2f), %d, %.2f\n",irow,xyz[0],xyz[1],xyz[2],sec,t0);
300 }
526ddf0e 301
32438f4e 302 delete splXY;
303 splXY=0x0;
304 delete splXZ;
305 splXZ=0x0;
306
307}
526ddf0e 308
32438f4e 309//________________________________________________________________
310Bool_t AliToyMCEventGenerator::SetupCluster(AliTPCclusterMI &tempCl, Float_t xyz[3], Int_t sec, Double_t t0)
311{
312 //
313 //
314 //
526ddf0e 315
32438f4e 316 const Double_t kSigmaY = 0.1;
317 const Double_t kSigmaZ = 0.1;
318 const Double_t kMaxZ0 = fTPCParam->GetZLength();
319 //TODO: Get this from the OCDB at some point?
320 const Double_t kDriftVel = fTPCParam->GetDriftV();
526ddf0e 321
32438f4e 322 tempCl.SetX(xyz[0]);
323 tempCl.SetY(xyz[1]);
324 tempCl.SetZ(xyz[2]);
526ddf0e 325
32438f4e 326 tempCl.SetSigmaY2(kSigmaY*kSigmaY);
327 tempCl.SetSigmaZ2(kSigmaZ*kSigmaZ);
328
329 // transform from the local coordinates to the coordinates expressed in pad coordinates
330 Int_t index[3] = {0,sec,0};
331 fTPCParam->Transform2to3(xyz,index);
332 fTPCParam->Transform3to4(xyz,index);
333 fTPCParam->Transform4to8(xyz,index);
334
335 const Int_t row = index[2];
336 const Int_t nPads = fTPCParam->GetNPads(sec, row);
337 // pad is fractional, but it needs to be shifted from the center
338 // to the edge of the row
339 const Float_t pad = xyz[1] + nPads/2;
340
341 tempCl.SetRow(row);
342 tempCl.SetPad(pad);
343 Float_t timeBin=Float_t(t0 + (kMaxZ0-TMath::Abs(tempCl.GetZ()))/kDriftVel);
344 tempCl.SetTimeBin(timeBin); // set time as t0 + drift time from dist z
345 tempCl.SetDetector(sec);
346
347 //check if we are in the active area
348 if (pad<0 || pad>=nPads) return kFALSE;
349
350 return kTRUE;
526ddf0e 351}
32438f4e 352
a1a695e5 353//________________________________________________________________
354Bool_t AliToyMCEventGenerator::ConnectOutputFile()
355{
356 //
357 // Create the output file name and tree and connect the event
358 //
359
360 fOutFile = new TFile(fOutputFileName.Data(),"recreate");
361
362 if (!fOutFile || !fOutFile->IsOpen()){
363 delete fOutFile;
364 fOutFile=0x0;
365 return kFALSE;
366 }
367
368 fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
369 fOutTree->Branch("event","AliToyMCEvent",&fEvent);
370
371 return kTRUE;
372}
373
374//________________________________________________________________
375Bool_t AliToyMCEventGenerator::CloseOutputFile()
376{
377 //
378 // close the output file
379 //
380 if (!fOutFile) return kFALSE;
381 fOutFile->Write();
382 fOutFile->Close();
383 delete fOutFile;
384 fOutFile=0x0;
385
386 return kTRUE;
387}
388
389//________________________________________________________________
390void AliToyMCEventGenerator::FillTree()
391{
392 // fill the tree
393 if (fOutTree) fOutTree->Fill();
394}
395
cd8ed0ac 396//________________________________________________________________
397void AliToyMCEventGenerator::SetSpaceCharge(EEpsilon epsilon, EGasType gasType/*=kNeCO2_9010*/, ECollRate collRate/*=k50kHz*/)
398{
399 //
400 // Set the space charge conditions
401 //
402 fSpaceChargeFile="$ALICE_ROOT/TPC/Calib/maps/SC";
403 switch (gasType) {
404 case kNeCO2_9010:
405 fSpaceChargeFile.Append("_NeCO2");
406 break;
407 }
408 switch (epsilon) {
409 case kEps5:
410 fSpaceChargeFile.Append("_eps5");
411 break;
412 case kEps10:
413 fSpaceChargeFile.Append("_eps10");
414 break;
415 case kEps20:
416 fSpaceChargeFile.Append("_eps20");
417 break;
418 }
419 switch (collRate) {
420 case k50kHz:
421 fSpaceChargeFile.Append("_50kHz");
422 break;
423 }
424 fSpaceChargeFile.Append("_precal.root");
425}
426
427//________________________________________________________________
428void AliToyMCEventGenerator::InitSpaceCharge()
429{
430 //
431 // init the space charge conditions
432 // this should be called after the tree was connected
433 //
434
435 AliInfo(Form("Using space charge map file: '%s'",fSpaceChargeFile.Data()));
436
437 TFile f(fSpaceChargeFile.Data());
438 fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
439
440 if (fOutTree){
441 AliInfo("Attaching space charge map file name to the tree");
442 fOutTree->GetUserInfo()->Add(new TObjString(fSpaceChargeFile.Data()));
443 }
444}