TPC/Upgrade/AliToyMCEventGenerator.cxx

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