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