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 }