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 #include <TLinearFitter.h>
  14
  15 #include <AliLog.h>
  16 #include <AliTPCROC.h>
  17 #include <AliTrackPointArray.h>
  18 #include <AliTrackerBase.h>
  19 #include <AliCDBManager.h>
  20 #include <AliTPCParam.h>
  21 #include <AliGeomManager.h>
  22 #include <AliTPCcalibDB.h>
  23 #include <AliTPCclusterMI.h>
  24 #include <AliTPCSpaceCharge3D.h>
  25 #include <AliTPCROC.h>
  26 #include <AliExternalTrackParam.h>
  27
  28 #include "AliToyMCEvent.h"
  29 #include "AliToyMCTrack.h"
  30
  31 #include "AliToyMCEventGenerator.h"
  32
  33 ClassImp(AliToyMCEventGenerator);
  34
  35
  36 AliToyMCEventGenerator::AliToyMCEventGenerator()
  37   :TObject()
  38   ,fTPCParam(0x0)
  39   ,fEvent(0x0)
  40   ,fCurrentTrack(0)
  41   ,fTPCCorrection(0x0)
  42   ,fTPCCorrectionAv(0x0)
  43   ,fSCList(0x0)
  44   ,fSCListFile()
  45   ,fCorrectionFile("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.lookup.root")
  46   ,fOutputFileName("toyMC.root")
  47   ,fOutFile(0x0)
  48   ,fOutTree(0x0)
  49   ,fUseStepCorrection(kFALSE)
  50   ,fUseMaterialBudget(kFALSE)
  51   ,fIsLaser(kTRUE)
  52   ,fPrereadSCList(kFALSE)
  53 {
  54   fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  55   fTPCParam->ReadGeoMatrices();
  56   gRandom->SetSeed();
  57 }
  58 //________________________________________________________________
  59 AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
  60   :TObject(gen)
  61   ,fTPCParam(gen.fTPCParam)
  62   ,fEvent(0x0)
  63   ,fCurrentTrack(0)
  64   ,fTPCCorrection(gen.fTPCCorrection)
  65   ,fTPCCorrectionAv(0x0)
  66   ,fSCList(gen.fSCList)
  67   ,fSCListFile(gen.fSCListFile)
  68   ,fCorrectionFile(gen.fCorrectionFile)
  69   ,fOutputFileName(gen.fOutputFileName)
  70   ,fOutFile(0x0)
  71   ,fOutTree(0x0)
  72   ,fUseStepCorrection(gen.fUseStepCorrection)
  73   ,fUseMaterialBudget(gen.fUseMaterialBudget)
  74   ,fIsLaser(gen.fIsLaser)
  75   ,fPrereadSCList(gen.fPrereadSCList)
  76 {
  77   //
  78   gRandom->SetSeed();
  79 }
  80 //________________________________________________________________
  81 AliToyMCEventGenerator::~AliToyMCEventGenerator()
  82 {
  83   if (HasSCList() &&!fPrereadSCList) delete fTPCCorrection;
  84   delete fSCList;
  85 }
  86
  87 //________________________________________________________________
  88 Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
  89 {
  90   //
  91   //
  92   //
  93
  94   if(!fTPCParam) {
  95     fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  96     fTPCParam->ReadGeoMatrices();
  97   }
  98
  99   MakeITSClusters(trackIn/*,t0*/);
 100   MakeTPCClusters(trackIn, t0);
 101   MakeTRDClusters(trackIn/*,t0*/);
 102
 103   return kTRUE;
 104 }
 105 //________________________________________________________________
 106 void 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
 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
 119   const Double_t kMaxSnp = 0.85;
 120 //   const Double_t kMaxZ0  = fTPCParam->GetZLength();
 121   const Double_t kMass   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 122
 123   AliExternalTrackParam track(trackIn);
 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     }
 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
 141     track.GetXYZ(xyz);
 142
 143 //     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
 144     if (TMath::Abs(track.GetZ())>lengthITS[iLayer]/2) continue;
 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);
 154
 155     trackIn.AddITSPoint(point)->SetUniqueID(trackIn.GetUniqueID());
 156   }
 157
 158 }
 159 //________________________________________________________________
 160 void 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
 169   const Float_t sigmaY = 0.06;
 170   const Float_t sigmaZ = 0.2;
 171
 172   AliTrackPoint point;
 173
 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
 178   AliExternalTrackParam track(trackIn);
 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
 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());
 207   }
 208
 209 }
 210 //________________________________________________________________
 211 void AliToyMCEventGenerator::MakeTPCClusters(AliToyMCTrack &trackIn, Double_t t0)
 212 {
 213
 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;
 217   static AliTrackPointArray pointArray0(nMaxPoints);  //undistorted
 218   static AliTrackPointArray pointArray1(nMaxPoints);  //distorted
 219
 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);
 226
 227
 228 }
 229 //________________________________________________________________
 230 void 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
 239   Double_t kMaxSnp = 0.85;
 240   if (fIsLaser) kMaxSnp=0.99;
 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;
 246
 247   const Float_t kSigmaY=0.1;
 248   const Float_t kSigmaZ=0.1;
 249
 250   AliExternalTrackParam track(trackIn);
 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
 256   //!!! when does the propagation not work, how often does it happen?
 257   if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget) && !fIsLaser) {
 258     AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
 259     return;
 260   }
 261
 262   Int_t npoints=0;
 263
 264   for (Double_t radius=iFCRadius; radius<oFCRadius; radius+=stepSize){
 265     //!!! changed from return 0 to continue -> Please check
 266     if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
 267       AliError(Form("Propagation to r=%.2f (snp=%.2f) failed\n",radius,track.GetSnp()));
 268       continue;
 269     }
 270     track.GetXYZ(xyz);
 271
 272     //!!! Why is this smeared
 273
 274     xyzf[0]=Float_t(xyz[0]);
 275     xyzf[1]=Float_t(xyz[1]);
 276     xyzf[2]=Float_t(xyz[2]);
 277
 278     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
 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
 286     SetPoint(xyzf,kSigmaY,kSigmaZ,pUdist);
 287     arrUdist.AddPoint(npoints, &pUdist);
 288     Int_t sector=pUdist.GetVolumeID();
 289
 290     // set distorted point
 291     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 292     Float_t dxyz[3]={0.,0.,0.};
 293     if (!fUseStepCorrection){
 294       fTPCCorrection->DistortPoint(distPoint, sector);
 295     } else {
 296       fTPCCorrection->GetCorrectionIntegralDz(distPoint,sector,dxyz,5);
 297       distPoint[0]-=dxyz[0];
 298       distPoint[1]-=dxyz[1];
 299       distPoint[2]-=dxyz[2];
 300     }
 301     SetPoint(distPoint,kSigmaY,kSigmaZ,pDist);
 302     arrDist.AddPoint(npoints, &pDist);
 303
 304     ++npoints;
 305   }
 306
 307   arrUdist.SetUniqueID(npoints);
 308   arrDist.SetUniqueID(npoints);
 309 }
 310
 311 //________________________________________________________________
 312 void AliToyMCEventGenerator::SetPoint(Float_t xyz[3], Float_t sigmaY, Float_t sigmaZ, AliTrackPoint &point)
 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};
 321   Float_t cov[6]={0,0,0, sigmaY*sigmaY,0,sigmaZ*sigmaZ};
 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);
 337   // abuse volume ID for the sector number
 338   point.SetVolumeID(fTPCParam->Transform0to1(xyz,i));
 339
 340   // TODO: Add sampled dE/dx  (use SetCharge)
 341 }
 342
 343 //________________________________________________________________
 344 void AliToyMCEventGenerator::ConvertTrackPointsToLocalClusters(AliTrackPointArray &arrPoints,
 345                                                                AliToyMCTrack &tr, Double_t t0, Int_t type)
 346 {
 347   //
 348   //
 349   //
 350
 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     }
 391
 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   }
 397
 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;
 420
 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 {
 426         //TODO: make a cluster also in the sector w only one space point?
 427 //      continue;
 428         // Double_t tempX=0., tempY = 0., tempZ = 0.;
 429
 430         // grXY->GetPoint(0,tempX,localY);
 431         // grXZ->GetPoint(0,tempX,localZ);
 432 //       }
 433
 434     }
 435     secOld=sec;
 436
 437     // check we are in an active area
 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;
 440
 441     //get interpolated value at the center for the pad row
 442     //  using splines
 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"*/);
 447     Float_t xyz[3]={localX,localY,localZ};
 448
 449     if (!SetupCluster(tempCl,xyz,sec,t0)) continue;
 450     tempCl.SetLabel(tr.GetUniqueID(), 0);
 451
 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   }
 456
 457   delete splXY;
 458   splXY=0x0;
 459   delete splXZ;
 460   splXZ=0x0;
 461
 462 }
 463
 464 //________________________________________________________________
 465 Bool_t AliToyMCEventGenerator::SetupCluster(AliTPCclusterMI &tempCl, Float_t xyz[3], Int_t sec, Double_t t0)
 466 {
 467   //
 468   //
 469   //
 470
 471   // intrinsic cluster resolution is 1mm
 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();
 477
 478   // smear the ideal positions with the cluster resolution
 479   xyz[1]+=gRandom->Gaus(0,kSigmaY);
 480   xyz[2]+=gRandom->Gaus(0,kSigmaZ);
 481
 482   tempCl.SetX(xyz[0]);
 483   tempCl.SetY(xyz[1]);
 484   tempCl.SetZ(xyz[2]);
 485
 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;
 511 }
 512
 513 //________________________________________________________________
 514 Bool_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
 531   gROOT->cd();
 532
 533   return kTRUE;
 534 }
 535
 536 //________________________________________________________________
 537 Bool_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 //________________________________________________________________
 552 void AliToyMCEventGenerator::FillTree()
 553 {
 554   // fill the tree
 555   if (fOutTree&&fEvent) fOutTree->Fill();
 556 }
 557
 558 //________________________________________________________________
 559 void AliToyMCEventGenerator::SetSpaceCharge(EEpsilon epsilon, EGasType gasType/*=kNeCO2_9010*/,
 560                                             ECollRate collRate/*=k50kHz*/, ECorrection corrType/*=kLookup*/)
 561 {
 562   //
 563   // Set the space charge conditions
 564   //
 565   fCorrectionFile="$ALICE_ROOT/TPC/Calib/maps/SC";
 566   switch (gasType) {
 567     case kNeCO2_9010:
 568       fCorrectionFile.Append("_NeCO2");
 569       break;
 570     case kNeCO2N2_90105:
 571       fCorrectionFile.Append("_NeCO2N2");
 572       break;
 573   }
 574   switch (epsilon) {
 575     case kEps5:
 576       fCorrectionFile.Append("_eps5");
 577       break;
 578     case kEps10:
 579       fCorrectionFile.Append("_eps10");
 580       break;
 581     case kEps20:
 582       fCorrectionFile.Append("_eps20");
 583       break;
 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;
 596   }
 597   switch (collRate) {
 598     case k50kHz:
 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");
 608       break;
 609   }
 610 }
 611
 612 //________________________________________________________________
 613 void AliToyMCEventGenerator::InitSpaceCharge()
 614 {
 615   //
 616   // init the space charge conditions
 617   // this should be called after the tree was connected
 618   //
 619
 620   if (HasSCList()) {
 621     InitSpaceChargeList();
 622     return;
 623   }
 624
 625   AliInfo(Form("Using space charge map file: '%s'",fCorrectionFile.Data()));
 626
 627   SetCorrectionFromFile(fCorrectionFile, fTPCCorrection);
 628
 629   if (fOutTree){
 630     AliInfo("Attaching space charge map file name to the tree");
 631     fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
 632   }
 633 }
 634
 635 //________________________________________________________________
 636 void 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
 666   // check for an average map
 667   if (!fCorrectionFile.IsNull()) {
 668     SetCorrectionFromFile(fCorrectionFile, fTPCCorrectionAv);
 669   }
 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) {
 686     TString scname=arr->At(ifile)->GetName();
 687     AliTPCCorrection *sc=0x0;
 688     SetCorrectionFromFile(scname, sc);
 689     if (!sc) continue;
 690     fSCList->Add(sc);
 691   }
 692 }
 693
 694 //________________________________________________________________
 695 void 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
 721   SetCorrectionFromFile(file, fTPCCorrection);
 722
 723   if (!fTPCCorrection) {
 724     AliFatal(Form("Could not read SC map from SC file '%s'",file.Data()));
 725     return;
 726   }
 727
 728 }
 729
 730 //________________________________________________________________
 731 Float_t AliToyMCEventGenerator::GetSCScalingFactor(AliTPCCorrection *corr, AliTPCCorrection *averageCorr, Float_t &chi2)
 732 {
 733   //
 734   // calculate the scaling factor
 735   // between the fluctuation map and the average map
 736   //
 737
 738   TLinearFitter fitter(2,"pol1");
 739
 740   for (Float_t iz=-245; iz<=245; iz+=10) {
 741     Short_t roc=(iz>=0)?0:18;
 742     for (Float_t ir=86; ir<250; ir+=10) {
 743       for (Float_t iphi=0; iphi<TMath::TwoPi(); iphi+=10*TMath::DegToRad()){
 744         Float_t x=ir*(Float_t)TMath::Cos(iphi);
 745         Float_t y=ir*(Float_t)TMath::Sin(iphi);
 746         Float_t x3[3]    = {x,y,iz};
 747         Float_t dx3[3]   = {0.,0.,0.};
 748         Float_t dx3av[3] = {0.,0.,0.};
 749
 750         corr->GetDistortion(x3,roc,dx3);
 751         averageCorr->GetDistortion(x3,roc,dx3av);
 752
 753         Double_t dr   = TMath::Sqrt(dx3[0]*dx3[0]     +  dx3[1]*dx3[1]);
 754         Double_t drav = TMath::Sqrt(dx3av[0]*dx3av[0] +  dx3av[1]*dx3av[1]);
 755
 756         fitter.AddPoint(&drav,dr);
 757       }
 758     }
 759   }
 760
 761   fitter.Eval();
 762   chi2 = fitter.GetChisquare()/fitter.GetNpoints();
 763   return fitter.GetParameter(1);
 764 }
 765
 766 //________________________________________________________________
 767 void AliToyMCEventGenerator::SetSCScalingFactor()
 768 {
 769   //
 770   // if running with a SC list calculate the scaling factor
 771   // between the fluctuation map and the average map
 772   //
 773
 774   if ( !(HasSCList() && fTPCCorrection && fTPCCorrectionAv && fEvent) ) return;
 775
 776   // loop over several z, r and phi bins and find a factor that minimises
 777   // the distortions between the current map and the average map
 778
 779   Float_t chi2   = 0.;
 780   Float_t factor = GetSCScalingFactor(fTPCCorrection, fTPCCorrectionAv, chi2);
 781
 782   fEvent->SetSCscale(factor);
 783   fEvent->SetSCscaleChi2(chi2);
 784 }
 785
 786 //________________________________________________________________
 787 void AliToyMCEventGenerator::SetCorrectionFromFile(TString file, AliTPCCorrection* &corr)
 788 {
 789   //
 790   // read the correction from file and set it to corr
 791   //
 792
 793   corr=0x0;
 794   TString corrName("map");
 795
 796   // allow for specifying an object name for the AliTPCCorrection in the file name
 797   // separated by a ':'
 798   TObjArray *arr=file.Tokenize(":");
 799   if (arr->GetEntriesFast()>1) {
 800     file=arr->At(0)->GetName();
 801     corrName=arr->At(1)->GetName();
 802   }
 803   delete arr;
 804
 805
 806   TFile f(file.Data());
 807   if (!f.IsOpen()||f.IsZombie()) {
 808     printf("E-AliToyMCEventGenerator::Could not open SC file '%s'",file.Data());
 809     return;
 810   }
 811
 812   corr=(AliTPCSpaceCharge3D*)f.Get(corrName.Data());
 813   if (corr) corr->SetName(f.GetName());
 814 }
 815
 816