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