]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/Upgrade/AliToyMCEventGenerator.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
o add usage of SC lists
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGenerator.cxx
CommitLineData
526ddf0e 1#include <iostream>
2dd02504 2#include <fstream>
a1a695e5 3
526ddf0e 4#include <TDatabasePDG.h>
5#include <TRandom.h>
6#include <TH2F.h>
a1a695e5 7#include <TGeoGlobalMagField.h>
32438f4e 8#include <TSpline.h>
cd8ed0ac 9#include <TObjString.h>
0403120d 10#include <TROOT.h>
2dd02504 11#include <TSystem.h>
12#include <TObjArray.h>
a1a695e5 13
b0c26b55 14#include <AliLog.h>
a1a695e5 15#include <AliTPCROC.h>
16#include <AliTrackPointArray.h>
526ddf0e 17#include <AliTrackerBase.h>
18#include <AliCDBManager.h>
19#include <AliTPCParam.h>
20#include <AliGeomManager.h>
21#include <AliTPCcalibDB.h>
526ddf0e 22#include <AliTPCclusterMI.h>
a1a695e5 23#include <AliTPCSpaceCharge3D.h>
32438f4e 24#include <AliTPCROC.h>
0403120d 25#include <AliExternalTrackParam.h>
a1a695e5 26
27#include "AliToyMCEvent.h"
28#include "AliToyMCTrack.h"
29
30#include "AliToyMCEventGenerator.h"
31
de0014b7 32ClassImp(AliToyMCEventGenerator);
526ddf0e 33
34
de0014b7 35AliToyMCEventGenerator::AliToyMCEventGenerator()
526ddf0e 36 :TObject()
37 ,fTPCParam(0x0)
a1a695e5 38 ,fEvent(0x0)
223d9e38 39 ,fCurrentTrack(0)
40 ,fTPCCorrection(0x0)
2dd02504 41 ,fSCList(0x0)
42 ,fSCListFile()
223d9e38 43 ,fCorrectionFile("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.lookup.root")
a1a695e5 44 ,fOutputFileName("toyMC.root")
45 ,fOutFile(0x0)
46 ,fOutTree(0x0)
d92b8630 47 ,fUseStepCorrection(kFALSE)
d1cf83f5 48 ,fUseMaterialBudget(kFALSE)
0403120d 49 ,fIsLaser(kTRUE)
2dd02504 50 ,fPrereadSCList(kFALSE)
526ddf0e 51{
526ddf0e 52 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
53 fTPCParam->ReadGeoMatrices();
223d9e38 54 gRandom->SetSeed();
526ddf0e 55}
56//________________________________________________________________
de0014b7 57AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
526ddf0e 58 :TObject(gen)
a1a695e5 59 ,fTPCParam(gen.fTPCParam)
60 ,fEvent(0x0)
223d9e38 61 ,fCurrentTrack(0)
62 ,fTPCCorrection(gen.fTPCCorrection)
2dd02504 63 ,fSCList(gen.fSCList)
64 ,fSCListFile(gen.fSCListFile)
223d9e38 65 ,fCorrectionFile(gen.fCorrectionFile)
a1a695e5 66 ,fOutputFileName(gen.fOutputFileName)
67 ,fOutFile(0x0)
68 ,fOutTree(0x0)
d92b8630 69 ,fUseStepCorrection(gen.fUseStepCorrection)
d1cf83f5 70 ,fUseMaterialBudget(gen.fUseMaterialBudget)
0403120d 71 ,fIsLaser(gen.fIsLaser)
2dd02504 72 ,fPrereadSCList(gen.fPrereadSCList)
526ddf0e 73{
74 //
223d9e38 75 gRandom->SetSeed();
526ddf0e 76}
77//________________________________________________________________
de0014b7 78AliToyMCEventGenerator::~AliToyMCEventGenerator()
526ddf0e 79{
2dd02504 80 if (HasSCList() &&!fPrereadSCList) delete fTPCCorrection;
81 delete fSCList;
526ddf0e 82}
32438f4e 83
526ddf0e 84//________________________________________________________________
a1a695e5 85Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
86{
87 //
88 //
89 //
90
526ddf0e 91 if(!fTPCParam) {
92 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
93 fTPCParam->ReadGeoMatrices();
526ddf0e 94 }
32438f4e 95
1e62e876 96 MakeITSClusters(trackIn/*,t0*/);
97 MakeTPCClusters(trackIn, t0);
98 MakeTRDClusters(trackIn/*,t0*/);
99
0403120d 100 return kTRUE;
1e62e876 101}
102//________________________________________________________________
103void AliToyMCEventGenerator::MakeITSClusters(AliToyMCTrack &trackIn/*, Double_t t0*/)
104{
105 //Upgrade ITS parameters
106 const Int_t nITSLayers = 7;
107 const Double_t ITSRadii[nITSLayers] = {2.2, 2.8, 3.6, 20.0, 22.0, 41.0, 43.0};
108 const Double_t lengthITS[nITSLayers] = {22.4, 24.2, 26.8, 78.0, 83.6, 142.4, 148.6};
109
e83fd282 110 //resolution of the point is 10um
111 const Float_t sigmaY = 0.001;
112 const Float_t sigmaZ = 0.001;
113
114 AliTrackPoint point;
115
1e62e876 116 const Double_t kMaxSnp = 0.85;
e83fd282 117// const Double_t kMaxZ0 = fTPCParam->GetZLength();
1e62e876 118 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
119
0403120d 120 AliExternalTrackParam track(trackIn);
1e62e876 121 Double_t xyz[3] = {0.,0.,0.};
122
123 if (!AliTrackerBase::PropagateTrackTo(&track,ITSRadii[0],kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
124 AliError(Form("Propagation to %.2f failed\n",ITSRadii[0]));
125 return;
126 }
127
128 for(Int_t iLayer = 0; iLayer<nITSLayers; iLayer++){
129
130 if (!AliTrackerBase::PropagateTrackTo(&track,ITSRadii[iLayer],kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
131 AliError(Form("Propagation to %.2f failed\n",ITSRadii[iLayer]));
132 continue;
133 }
e83fd282 134 // since I don't know how to extract the volumeid of the ITS layers I use the following strategy:
135 // - rotate the track to the next integer angle
136 // - use coordinates there and set as volumeid the integer angle
137 // - in the reco one can then rotate the cluster to the global frame using the angle stored in the volume id
1e62e876 138 track.GetXYZ(xyz);
139
e83fd282 140// if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
1e62e876 141 if (TMath::Abs(track.GetZ())>lengthITS[iLayer]/2) continue;
e83fd282 142
143
144 // smear the ideal positions with the cluster resolution
145 xyz[1]+=gRandom->Gaus(0,sigmaY);
146 xyz[2]+=gRandom->Gaus(0,sigmaZ);
147
148 Float_t xyzf[3]={xyz[0],xyz[1],xyz[2]};
149
150 SetPoint(xyzf,sigmaY,sigmaZ,point);
1e62e876 151
e83fd282 152 trackIn.AddITSPoint(point)->SetUniqueID(trackIn.GetUniqueID());
1e62e876 153 }
154
155}
156//________________________________________________________________
157void AliToyMCEventGenerator::MakeTRDClusters(AliToyMCTrack &trackIn/*, Double_t t0*/)
158
159{
160 //Uses current TRD parameters
161 const Int_t nTRDLayers = 6;
162 const Double_t distToMid = 3.2 + 30./2; //dist to middle of drift region (radiator + half drift region)
163 const Double_t TRDRadii[nTRDLayers] = {294.5 + distToMid, 307.1 + distToMid, 319.7 + distToMid, 332.3 + distToMid, 344.9 + distToMid, 357.5 + distToMid};
164 const Double_t lengthTRD[nTRDLayers] = {604.0, 634.0, 656.0, 686.0, 700.0, 700.0};
165
e83fd282 166 const Float_t sigmaY = 0.06;
167 const Float_t sigmaZ = 0.2;
1e62e876 168
e83fd282 169 AliTrackPoint point;
170
1e62e876 171 const Double_t kMaxSnp = 0.85;
172 const Double_t kMaxZ0 = fTPCParam->GetZLength();
173 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
174
0403120d 175 AliExternalTrackParam track(trackIn);
1e62e876 176 Double_t xyz[3] = {0.,0.,0.};
177
178 if (!AliTrackerBase::PropagateTrackTo(&track,TRDRadii[0],kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
179 AliError(Form("Propagation to %.2f failed\n",TRDRadii[0]));
180 return;
181 }
182
183 for(Int_t iLayer = 0; iLayer<nTRDLayers; iLayer++){
184
185 if (!AliTrackerBase::PropagateTrackTo(&track,TRDRadii[iLayer],kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
186 AliError(Form("Propagation to %.2f failed\n",TRDRadii[iLayer]));
187 continue;
188 }
189 track.GetXYZ(xyz);
190
191 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
192 if (TMath::Abs(track.GetZ())>lengthTRD[iLayer]/2) continue;
193
e83fd282 194
195 // smear the ideal positions with the cluster resolution
196 xyz[1]+=gRandom->Gaus(0,sigmaY);
197 xyz[2]+=gRandom->Gaus(0,sigmaZ);
198
199 Float_t xyzf[3]={xyz[0],xyz[1],xyz[2]};
200
201 SetPoint(xyzf,sigmaY,sigmaZ,point);
202
203 trackIn.AddTRDPoint(point)->SetUniqueID(trackIn.GetUniqueID());
1e62e876 204 }
205
206}
207//________________________________________________________________
208void AliToyMCEventGenerator::MakeTPCClusters(AliToyMCTrack &trackIn, Double_t t0)
209{
210
32438f4e 211 // make it big enough to hold all points
212 // store real number of generated points in the unique id
213 const Int_t nMaxPoints=3000;
0403120d 214 static AliTrackPointArray pointArray0(nMaxPoints); //undistorted
215 static AliTrackPointArray pointArray1(nMaxPoints); //distorted
1e62e876 216
32438f4e 217 //Create space point of undistorted and distorted clusters along the propagated track trajectory
218 CreateSpacePoints(trackIn,pointArray0,pointArray1);
219 //Convert the space points into clusters in the local frame
220 //for undistorted and distorted clusters using the same function
221 ConvertTrackPointsToLocalClusters(pointArray0,trackIn,t0,0);
222 ConvertTrackPointsToLocalClusters(pointArray1,trackIn,t0,1);
526ddf0e 223
526ddf0e 224
1e62e876 225}
32438f4e 226//________________________________________________________________
227void AliToyMCEventGenerator::CreateSpacePoints(AliToyMCTrack &trackIn,
228 AliTrackPointArray &arrUdist,
229 AliTrackPointArray &arrDist)
230{
231 //
232 // sample the track from the inner to the outer wall of the TPC
233 // a graph is filled in local coordinates for later
234 //
235
0403120d 236 Double_t kMaxSnp = 0.85;
237 if (fIsLaser) kMaxSnp=0.99;
32438f4e 238 const Double_t kMaxZ0 = fTPCParam->GetZLength();
239 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
240
241 const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
242 const Double_t oFCRadius = 254.5;
e83fd282 243
244 const Float_t kSigmaY=0.1;
245 const Float_t kSigmaZ=0.1;
32438f4e 246
0403120d 247 AliExternalTrackParam track(trackIn);
32438f4e 248 //!!! TODO: make this adjustable perhaps
249 const Double_t stepSize=0.1;
250 Double_t xyz[3] = {0.,0.,0.};
251 Float_t xyzf[3] = {0.,0.,0.};
252
a1a695e5 253 //!!! when does the propagation not work, how often does it happen?
0403120d 254 if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget) && !fIsLaser) {
b0c26b55 255 AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
32438f4e 256 return;
b0c26b55 257 }
32438f4e 258
526ddf0e 259 Int_t npoints=0;
b0c26b55 260
32438f4e 261 for (Double_t radius=iFCRadius; radius<oFCRadius; radius+=stepSize){
b0c26b55 262 //!!! changed from return 0 to continue -> Please check
d1cf83f5 263 if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
0403120d 264 AliError(Form("Propagation to r=%.2f (snp=%.2f) failed\n",radius,track.GetSnp()));
b0c26b55 265 continue;
266 }
526ddf0e 267 track.GetXYZ(xyz);
32438f4e 268
b0c26b55 269 //!!! Why is this smeared
32438f4e 270
271 xyzf[0]=Float_t(xyz[0]);
272 xyzf[1]=Float_t(xyz[1]);
273 xyzf[2]=Float_t(xyz[2]);
274
526ddf0e 275 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
32438f4e 276 // if (TMath::Abs(track.GetX())<iFCRadius) continue;
277 // if (TMath::Abs(track.GetX())>oFCRadius) continue;
278
279
280 AliTrackPoint pUdist; // undistorted space point
281 AliTrackPoint pDist; // distorted space point
282 // Set undistorted point
e83fd282 283 SetPoint(xyzf,kSigmaY,kSigmaZ,pUdist);
32438f4e 284 arrUdist.AddPoint(npoints, &pUdist);
285 Int_t sector=pUdist.GetVolumeID();
32438f4e 286
287 // set distorted point
526ddf0e 288 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
d92b8630 289 Float_t dxyz[3]={0.,0.,0.};
290 if (!fUseStepCorrection){
223d9e38 291 fTPCCorrection->DistortPoint(distPoint, sector);
d92b8630 292 } else {
223d9e38 293 fTPCCorrection->GetCorrectionIntegralDz(distPoint,sector,dxyz,5);
d92b8630 294 distPoint[0]-=dxyz[0];
295 distPoint[1]-=dxyz[1];
296 distPoint[2]-=dxyz[2];
297 }
e83fd282 298 SetPoint(distPoint,kSigmaY,kSigmaZ,pDist);
32438f4e 299 arrDist.AddPoint(npoints, &pDist);
32438f4e 300
301 ++npoints;
526ddf0e 302 }
32438f4e 303
304 arrUdist.SetUniqueID(npoints);
305 arrDist.SetUniqueID(npoints);
306}
526ddf0e 307
32438f4e 308//________________________________________________________________
e83fd282 309void AliToyMCEventGenerator::SetPoint(Float_t xyz[3], Float_t sigmaY, Float_t sigmaZ, AliTrackPoint &point)
32438f4e 310{
311 //
312 // make AliTrackPoint out of AliTPCclusterMI
313 //
314
315 //covariance at the local frame
316 //assume 1mm distortion in y and z
317 Int_t i[3]={0,0,0};
e83fd282 318 Float_t cov[6]={0,0,0, sigmaY*sigmaY,0,sigmaZ*sigmaZ};
32438f4e 319
320 const Float_t alpha = -TMath::ATan2(xyz[1],xyz[0]);
321 const Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
322
323 Float_t newcov[6];
324 newcov[0] = cov[0]*cos*cos + 2*cov[1]*sin*cos + cov[3]*sin*sin;
325 newcov[1] = cov[1]*(cos*cos-sin*sin) + (cov[3]-cov[0])*sin*cos;
326 newcov[2] = cov[2]*cos + cov[4]*sin;
327 newcov[3] = cov[0]*sin*sin - 2*cov[1]*sin*cos + cov[3]*cos*cos;
328 newcov[4] = cov[4]*cos - cov[2]*sin;
329 newcov[5] = cov[5];
330
331 // voluem ID to add later ....
332 point.SetXYZ(xyz);
333 point.SetCov(newcov);
d92b8630 334 // abuse volume ID for the sector number
32438f4e 335 point.SetVolumeID(fTPCParam->Transform0to1(xyz,i));
b0c26b55 336
32438f4e 337 // TODO: Add sampled dE/dx (use SetCharge)
338}
b0c26b55 339
32438f4e 340//________________________________________________________________
341void AliToyMCEventGenerator::ConvertTrackPointsToLocalClusters(AliTrackPointArray &arrPoints,
342 AliToyMCTrack &tr, Double_t t0, Int_t type)
343{
344 //
345 //
346 //
526ddf0e 347
32438f4e 348 const Int_t npoints=Int_t(arrPoints.GetUniqueID());
349 Int_t secOld=-1;
350
351 // create an array for graphs which are used for local interpolation
352 // we need a new graph if the sector changed since then the local frame will also change
353 TObjArray arrGraphsXY(72);
354 arrGraphsXY.SetOwner();
355 TObjArray arrGraphsXZ(72);
356 arrGraphsXZ.SetOwner();
357 AliTrackPoint p;
358 //create initial graph
359 TGraph *grXY=0x0;
360 TGraph *grXZ=0x0;
361 //row -> sector mapping
362 Int_t rowMap[159];
363 for (Int_t irow=0; irow<159; ++irow) rowMap[irow]=-1;
364
365 // 1. Step
366 // Make from the list of global space points graphs in the local frame
367 // one graph per sector is needed
368 for (Int_t ipoint=0; ipoint<npoints; ++ipoint){
369 arrPoints.GetPoint(p,ipoint);
370 Float_t xyz[3] = {p.GetX(),p.GetY(),p.GetZ()};
371 Int_t index[3] = {0,0,0};
372 Int_t row = fTPCParam->GetPadRow(xyz,index);
373 // rotate space point to local frame
374 // the angle is given by the VolumeID which was set in CrateSpacePoints
375 const Int_t sec=p.GetVolumeID();
376 if (row<0 || (sec<36 && row>62) || row>95 ) continue;
377 Double_t angle=((sec%18)*20.+10.)/TMath::RadToDeg();
378 AliTrackPoint pRot=p.Rotate(angle);
379 const Int_t secrow=row+(sec>35)*63;
380 if (rowMap[secrow]==-1) rowMap[secrow]=sec;
381 // check if we need a new graph (sector change)
382 if (secOld!=sec){
383 grXY=new TGraph;
384 grXZ=new TGraph;
385 arrGraphsXY.AddAt(grXY,sec);
386 arrGraphsXZ.AddAt(grXZ,sec);
387 }
b0c26b55 388
32438f4e 389 //add coordinates in local frame for later interpolation
390 grXY->SetPoint(grXY->GetN(), pRot.GetX(), pRot.GetY());
391 grXZ->SetPoint(grXZ->GetN(), pRot.GetX(), pRot.GetZ());
392 secOld=sec;
393 }
b0c26b55 394
32438f4e 395 // 2. Step
396 // create in the center of each row a space point by using the graph to interpolate
397 // the the center of the row. This is done in xy and xz
398 TSpline3 *splXY=0x0;
399 TSpline3 *splXZ=0x0;
400 AliTPCclusterMI tempCl;
401 secOld=-1;
402 for (Int_t irow=0; irow<159; ++irow ){
403 const Int_t sec = rowMap[irow];
404 if (sec==-1) continue;
405 const Int_t secrow = irow<63?irow:irow-63;
406 Double_t localX = fTPCParam->GetPadRowRadii(sec,secrow);
407 // get graph for the current row
408 if (sec!=secOld){
409 delete splXY;
410 splXY=0x0;
411 delete splXZ;
412 splXZ=0x0;
413
414 grXY=(TGraph*)arrGraphsXY.At(sec);
415 grXZ=(TGraph*)arrGraphsXZ.At(sec);
416 if (!grXY) continue;
b0c26b55 417
05da1b4e 418// if(grXY->GetN()>1 && grXZ->GetN()>1) { //causes segmentation violation if N==1
419// splXY=new TSpline3("splXY",grXY);
420// splXZ=new TSpline3("splXZ",grXZ);
421// }
422// else {
544076ea 423 //TODO: make a cluster also in the sector w only one space point?
05da1b4e 424// continue;
544076ea 425 // Double_t tempX=0., tempY = 0., tempZ = 0.;
426
427 // grXY->GetPoint(0,tempX,localY);
428 // grXZ->GetPoint(0,tempX,localZ);
05da1b4e 429// }
544076ea 430
32438f4e 431 }
432 secOld=sec;
b0c26b55 433
32438f4e 434 // check we are in an active area
05da1b4e 435// if (splXY->FindX(localX)<1 || splXZ->FindX(localX)<1) continue;
436 if ( localX<grXY->GetX()[0] || localX>grXY->GetX()[grXY->GetN()-1] || localX<grXZ->GetX()[0] || localX>grXZ->GetX()[grXZ->GetN()-1]) continue;
b0c26b55 437
32438f4e 438 //get interpolated value at the center for the pad row
439 // using splines
05da1b4e 440// const Double_t localY=splXY->Eval(localX/*,0x0,"S"*/);
441// const Double_t localZ=splXZ->Eval(localX/*,0x0,"S"*/);
442 const Double_t localY=grXY->Eval(localX/*,0x0,"S"*/);
443 const Double_t localZ=grXZ->Eval(localX/*,0x0,"S"*/);
32438f4e 444 Float_t xyz[3]={localX,localY,localZ};
526ddf0e 445
32438f4e 446 if (!SetupCluster(tempCl,xyz,sec,t0)) continue;
d1cf83f5 447 tempCl.SetLabel(tr.GetUniqueID(), 0);
0403120d 448
32438f4e 449 if (type==0) tr.AddSpacePoint(tempCl);
450 else tr.AddDistortedSpacePoint(tempCl);
451// printf("SetupCluster %3d: (%.2f, %.2f, %.2f), %d, %.2f\n",irow,xyz[0],xyz[1],xyz[2],sec,t0);
452 }
526ddf0e 453
32438f4e 454 delete splXY;
455 splXY=0x0;
456 delete splXZ;
457 splXZ=0x0;
458
459}
526ddf0e 460
32438f4e 461//________________________________________________________________
462Bool_t AliToyMCEventGenerator::SetupCluster(AliTPCclusterMI &tempCl, Float_t xyz[3], Int_t sec, Double_t t0)
463{
464 //
465 //
466 //
526ddf0e 467
e83fd282 468 // intrinsic cluster resolution is 1mm
32438f4e 469 const Double_t kSigmaY = 0.1;
470 const Double_t kSigmaZ = 0.1;
471 const Double_t kMaxZ0 = fTPCParam->GetZLength();
472 //TODO: Get this from the OCDB at some point?
473 const Double_t kDriftVel = fTPCParam->GetDriftV();
e83fd282 474
475 // smear the ideal positions with the cluster resolution
476 xyz[1]+=gRandom->Gaus(0,kSigmaY);
477 xyz[2]+=gRandom->Gaus(0,kSigmaZ);
526ddf0e 478
32438f4e 479 tempCl.SetX(xyz[0]);
480 tempCl.SetY(xyz[1]);
481 tempCl.SetZ(xyz[2]);
526ddf0e 482
32438f4e 483 tempCl.SetSigmaY2(kSigmaY*kSigmaY);
484 tempCl.SetSigmaZ2(kSigmaZ*kSigmaZ);
485
486 // transform from the local coordinates to the coordinates expressed in pad coordinates
487 Int_t index[3] = {0,sec,0};
488 fTPCParam->Transform2to3(xyz,index);
489 fTPCParam->Transform3to4(xyz,index);
490 fTPCParam->Transform4to8(xyz,index);
491
492 const Int_t row = index[2];
493 const Int_t nPads = fTPCParam->GetNPads(sec, row);
494 // pad is fractional, but it needs to be shifted from the center
495 // to the edge of the row
496 const Float_t pad = xyz[1] + nPads/2;
497
498 tempCl.SetRow(row);
499 tempCl.SetPad(pad);
500 Float_t timeBin=Float_t(t0 + (kMaxZ0-TMath::Abs(tempCl.GetZ()))/kDriftVel);
501 tempCl.SetTimeBin(timeBin); // set time as t0 + drift time from dist z
502 tempCl.SetDetector(sec);
503
504 //check if we are in the active area
505 if (pad<0 || pad>=nPads) return kFALSE;
506
507 return kTRUE;
526ddf0e 508}
32438f4e 509
a1a695e5 510//________________________________________________________________
511Bool_t AliToyMCEventGenerator::ConnectOutputFile()
512{
513 //
514 // Create the output file name and tree and connect the event
515 //
516
517 fOutFile = new TFile(fOutputFileName.Data(),"recreate");
518
519 if (!fOutFile || !fOutFile->IsOpen()){
520 delete fOutFile;
521 fOutFile=0x0;
522 return kFALSE;
523 }
524
525 fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
526 fOutTree->Branch("event","AliToyMCEvent",&fEvent);
527
0403120d 528 gROOT->cd();
529
a1a695e5 530 return kTRUE;
531}
532
533//________________________________________________________________
534Bool_t AliToyMCEventGenerator::CloseOutputFile()
535{
536 //
537 // close the output file
538 //
539 if (!fOutFile) return kFALSE;
540 fOutFile->Write();
541 fOutFile->Close();
542 delete fOutFile;
543 fOutFile=0x0;
544
545 return kTRUE;
546}
547
548//________________________________________________________________
549void AliToyMCEventGenerator::FillTree()
550{
551 // fill the tree
0403120d 552 if (fOutTree&&fEvent) fOutTree->Fill();
a1a695e5 553}
554
cd8ed0ac 555//________________________________________________________________
223d9e38 556void AliToyMCEventGenerator::SetSpaceCharge(EEpsilon epsilon, EGasType gasType/*=kNeCO2_9010*/,
557 ECollRate collRate/*=k50kHz*/, ECorrection corrType/*=kLookup*/)
cd8ed0ac 558{
559 //
560 // Set the space charge conditions
561 //
223d9e38 562 fCorrectionFile="$ALICE_ROOT/TPC/Calib/maps/SC";
cd8ed0ac 563 switch (gasType) {
564 case kNeCO2_9010:
223d9e38 565 fCorrectionFile.Append("_NeCO2");
cd8ed0ac 566 break;
d5715c32 567 case kNeCO2N2_90105:
568 fCorrectionFile.Append("_NeCO2N2");
569 break;
cd8ed0ac 570 }
571 switch (epsilon) {
572 case kEps5:
223d9e38 573 fCorrectionFile.Append("_eps5");
cd8ed0ac 574 break;
575 case kEps10:
223d9e38 576 fCorrectionFile.Append("_eps10");
cd8ed0ac 577 break;
578 case kEps20:
223d9e38 579 fCorrectionFile.Append("_eps20");
cd8ed0ac 580 break;
d5715c32 581 case kEps25:
582 fCorrectionFile.Append("_eps25");
583 break;
584 case kEps30:
585 fCorrectionFile.Append("_eps30");
586 break;
587 case kEps35:
588 fCorrectionFile.Append("_eps35");
589 break;
590 case kEps40:
591 fCorrectionFile.Append("_eps40");
592 break;
cd8ed0ac 593 }
594 switch (collRate) {
595 case k50kHz:
223d9e38 596 fCorrectionFile.Append("_50kHz");
597 break;
598 }
599 switch (corrType) {
600 case kLookup:
601 fCorrectionFile.Append("_precal.lookup.root");
602 break;
603 case kSpaceChargeFile:
604 fCorrectionFile.Append("_precal.root");
cd8ed0ac 605 break;
606 }
cd8ed0ac 607}
608
609//________________________________________________________________
610void AliToyMCEventGenerator::InitSpaceCharge()
611{
612 //
613 // init the space charge conditions
614 // this should be called after the tree was connected
615 //
616
2dd02504 617 if (HasSCList()) {
618 InitSpaceChargeList();
619 return;
620 }
621
223d9e38 622 AliInfo(Form("Using space charge map file: '%s'",fCorrectionFile.Data()));
623
624 TString corrName("map");
625
626 // allow for specifying an object name for the AliTPCCorrection in the file name
627 // separated by a ':'
628 TObjArray *arr=fCorrectionFile.Tokenize(":");
629 if (arr->GetEntriesFast()>1) {
630 fCorrectionFile=arr->At(0)->GetName();
631 corrName=arr->At(1)->GetName();
632 }
633 delete arr;
634
cd8ed0ac 635
223d9e38 636 TFile f(fCorrectionFile.Data());
637 fTPCCorrection=(AliTPCSpaceCharge3D*)f.Get("map");
cd8ed0ac 638
639 if (fOutTree){
640 AliInfo("Attaching space charge map file name to the tree");
223d9e38 641 fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
cd8ed0ac 642 }
643}
2dd02504 644
645//________________________________________________________________
646void AliToyMCEventGenerator::InitSpaceChargeList()
647{
648 //
649 // init space charge conditions from a list of files
650 // this should be called after the tree was connected
651 //
652
653 std::ifstream file(fSCListFile.Data());
654 TString list;
655 list.ReadFile(file);
656
657 TObjArray *arr=list.Tokenize("\n");
658 if (!arr->GetEntriesFast()) {
659 delete arr;
660 AliFatal(Form("No SC file initialised. SC list '%s' seems empty",fSCListFile.Data()));
661 return;
662 }
663
664 // it is assumed that in case of an input list
665 // fCorrectionFile contains the name of the average correction
666 // to be then used in the reconstruction
667 if (fOutTree){
668 if (fCorrectionFile.IsNull()) {
669 AliFatal("List of SC files set, but no average map is specified. Use 'SetSpaceChargeFile' to do so");
670 return;
671 }
672 AliInfo("Attaching average space charge map file name to the tree");
673 fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
674 }
675
676
677 // case of non preread
678 // store the names of the files
679 if (!fPrereadSCList) {
680 if (fSCList) delete fSCList;
681 fSCList=arr;
682 return;
683 }
684
685 // case of preread
686 // load all SC files and set them to the list
687 if (!fSCList) fSCList=new TObjArray;
688 fSCList->SetOwner();
689 fSCList->Delete();
690
691 for (Int_t ifile=0; ifile<arr->GetEntriesFast(); ++ifile) {
692 TFile f(arr->At(ifile)->GetName());
693 if (!f.IsOpen()||f.IsZombie()) continue;
694 AliTPCSpaceCharge3D *sc=(AliTPCSpaceCharge3D*)f.Get("map");
695 if (!sc) continue;
696 sc->SetName(f.GetName());
697 fSCList->Add(sc);
698 }
699}
700
701//________________________________________________________________
702void AliToyMCEventGenerator::IterateSC(Int_t ipos)
703{
704 //
705 // In case a list of SC files is set iterate over them
706 //
707
708 if (!HasSCList()) return;
709
710 if (ipos<0) {
711 if (fOutTree) ipos=fOutTree->GetEntries();
712 else ipos=0;
713 }
714
715 TObject *sc=fSCList->At(ipos%fSCList->GetEntriesFast());
716 AliInfo(Form("Event: %d - SC: %s",ipos,sc->GetName()));
717 // case SC files have been preread
718 if (fPrereadSCList) {
719 fTPCCorrection=(AliTPCCorrection*)sc;
720 return;
721 }
722
723 // case no preread was done
724 TString &file=((TObjString*)sc)->String();
725 delete fTPCCorrection;
726 fTPCCorrection=0x0;
727
728 TFile f(file);
729 if (!f.IsOpen()||f.IsZombie()) {
730 AliFatal(Form("Could not open SC file '%s'",file.Data()));
731 return;
732 }
733
734 fTPCCorrection=(AliTPCSpaceCharge3D*)f.Get("map");
735 if (!fTPCCorrection) {
736 AliFatal(Form("Could not read SC map from SC file '%s'",file.Data()));
737 return;
738 }
739
740}