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