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