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