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
  10 #include <AliLog.h>
  11 #include <AliTPCROC.h>
  12 #include <AliTrackPointArray.h>
  13 #include <AliTrackerBase.h>
  14 #include <AliCDBManager.h>
  15 #include <AliTPCParam.h>
  16 #include <AliGeomManager.h>
  17 #include <AliTPCcalibDB.h>
  18 #include <AliTPCclusterMI.h>
  19 #include <AliTPCSpaceCharge3D.h>
  20 #include <AliTPCROC.h>
  21
  22 #include "AliToyMCEvent.h"
  23 #include "AliToyMCTrack.h"
  24
  25 #include "AliToyMCEventGenerator.h"
  26
  27 ClassImp(AliToyMCEventGenerator);
  28
  29
  30 AliToyMCEventGenerator::AliToyMCEventGenerator()
  31   :TObject()
  32   ,fTPCParam(0x0)
  33   ,fEvent(0x0)
  34   ,fSpaceCharge(0x0)
  35   ,fSpaceChargeFile("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root")
  36   ,fOutputFileName("toyMC.root")
  37   ,fOutFile(0x0)
  38   ,fOutTree(0x0)
  39   ,fUseStepCorrection(kFALSE)
  40   ,fUseMaterialBudget(kFALSE)
  41 {
  42   fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  43   fTPCParam->ReadGeoMatrices();
  44 }
  45 //________________________________________________________________
  46 AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
  47   :TObject(gen)
  48   ,fTPCParam(gen.fTPCParam)
  49   ,fEvent(0x0)
  50   ,fSpaceCharge(gen.fSpaceCharge)
  51   ,fSpaceChargeFile(gen.fSpaceChargeFile)
  52   ,fOutputFileName(gen.fOutputFileName)
  53   ,fOutFile(0x0)
  54   ,fOutTree(0x0)
  55   ,fUseStepCorrection(gen.fUseStepCorrection)
  56   ,fUseMaterialBudget(gen.fUseMaterialBudget)
  57 {
  58   //
  59 }
  60 //________________________________________________________________
  61 AliToyMCEventGenerator::~AliToyMCEventGenerator()
  62 {
  63   delete fSpaceCharge;
  64 }
  65
  66 //________________________________________________________________
  67 Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
  68 {
  69   //
  70   //
  71   //
  72
  73   if(!fTPCParam) {
  74     fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  75     fTPCParam->ReadGeoMatrices();
  76   }
  77
  78   // make it big enough to hold all points
  79   // store real number of generated points in the unique id
  80   const Int_t nMaxPoints=3000;
  81   AliTrackPointArray pointArray0(nMaxPoints);  //undistorted
  82   AliTrackPointArray pointArray1(nMaxPoints);  //distorted
  83
  84   //Create space point of undistorted and distorted clusters along the propagated track trajectory
  85   CreateSpacePoints(trackIn,pointArray0,pointArray1);
  86   //Convert the space points into clusters in the local frame
  87   //for undistorted and distorted clusters using the same function
  88   ConvertTrackPointsToLocalClusters(pointArray0,trackIn,t0,0);
  89   ConvertTrackPointsToLocalClusters(pointArray1,trackIn,t0,1);
  90
  91   return 1;
  92 }
  93
  94 //________________________________________________________________
  95 void AliToyMCEventGenerator::CreateSpacePoints(AliToyMCTrack &trackIn,
  96                                                AliTrackPointArray &arrUdist,
  97                                                AliTrackPointArray &arrDist)
  98 {
  99   //
 100   // sample the track from the inner to the outer wall of the TPC
 101   // a graph is filled in local coordinates for later
 102   //
 103
 104   const Double_t kMaxSnp = 0.85;
 105   const Double_t kMaxZ0  = fTPCParam->GetZLength();
 106   const Double_t kMass   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 107
 108   const Double_t iFCRadius =  83.5; //radius constants found in AliTPCCorrection.cxx
 109   const Double_t oFCRadius = 254.5;
 110
 111   AliToyMCTrack track(trackIn);
 112   //!!! TODO: make this adjustable perhaps
 113   const Double_t stepSize=0.1;
 114   Double_t xyz[3]  = {0.,0.,0.};
 115   Float_t  xyzf[3] = {0.,0.,0.};
 116
 117   //!!! when does the propagation not work, how often does it happen?
 118   if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
 119     AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
 120     return;
 121   }
 122
 123   Int_t npoints=0;
 124
 125   for (Double_t radius=iFCRadius; radius<oFCRadius; radius+=stepSize){
 126
 127     //!!! changed from return 0 to continue -> Please check
 128     if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
 129       AliError(Form("Propagation to %.2f failed\n",radius));
 130       continue;
 131     }
 132     track.GetXYZ(xyz);
 133
 134     //!!! Why is this smeared
 135 //     xyz[0]+=gRandom->Gaus(0,0.000005);
 136 //     xyz[1]+=gRandom->Gaus(0,0.000005);
 137 //     xyz[2]+=gRandom->Gaus(0,0.000005);
 138
 139     xyzf[0]=Float_t(xyz[0]);
 140     xyzf[1]=Float_t(xyz[1]);
 141     xyzf[2]=Float_t(xyz[2]);
 142
 143     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
 144     //       if (TMath::Abs(track.GetX())<iFCRadius) continue;
 145     //       if (TMath::Abs(track.GetX())>oFCRadius) continue;
 146
 147
 148     AliTrackPoint pUdist;                               // undistorted space point
 149     AliTrackPoint pDist;                                // distorted space point
 150     // Set undistorted point
 151     SetPoint(xyzf,pUdist);
 152     arrUdist.AddPoint(npoints, &pUdist);
 153     Int_t sector=pUdist.GetVolumeID();
 154
 155     // set distorted point
 156     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 157     Float_t dxyz[3]={0.,0.,0.};
 158     if (!fUseStepCorrection){
 159       fSpaceCharge->DistortPoint(distPoint, sector);
 160     } else {
 161       fSpaceCharge->GetCorrectionIntegralDz(distPoint,sector,dxyz,5);
 162       distPoint[0]-=dxyz[0];
 163       distPoint[1]-=dxyz[1];
 164       distPoint[2]-=dxyz[2];
 165     }
 166     SetPoint(distPoint, pDist);
 167     arrDist.AddPoint(npoints, &pDist);
 168
 169     ++npoints;
 170   }
 171
 172   arrUdist.SetUniqueID(npoints);
 173   arrDist.SetUniqueID(npoints);
 174 }
 175
 176 //________________________________________________________________
 177 void AliToyMCEventGenerator::SetPoint(Float_t xyz[3], AliTrackPoint &point)
 178 {
 179   //
 180   // make AliTrackPoint out of AliTPCclusterMI
 181   //
 182
 183   //covariance at the local frame
 184   //assume 1mm distortion in y and z
 185   Int_t i[3]={0,0,0};
 186   const Double_t kSigmaY=0.1;
 187   const Double_t kSigmaZ=0.1;
 188   Float_t cov[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};
 189
 190   const Float_t alpha = -TMath::ATan2(xyz[1],xyz[0]);
 191   const Float_t sin   = TMath::Sin(alpha), cos = TMath::Cos(alpha);
 192
 193   Float_t newcov[6];
 194   newcov[0] = cov[0]*cos*cos           + 2*cov[1]*sin*cos         + cov[3]*sin*sin;
 195   newcov[1] = cov[1]*(cos*cos-sin*sin) + (cov[3]-cov[0])*sin*cos;
 196   newcov[2] = cov[2]*cos               + cov[4]*sin;
 197   newcov[3] = cov[0]*sin*sin           - 2*cov[1]*sin*cos          +  cov[3]*cos*cos;
 198   newcov[4] = cov[4]*cos               - cov[2]*sin;
 199   newcov[5] = cov[5];
 200
 201   // voluem ID to add later ....
 202   point.SetXYZ(xyz);
 203   point.SetCov(newcov);
 204   // abuse volume ID for the sector number
 205   point.SetVolumeID(fTPCParam->Transform0to1(xyz,i));
 206
 207   // TODO: Add sampled dE/dx  (use SetCharge)
 208 }
 209
 210 //________________________________________________________________
 211 void AliToyMCEventGenerator::ConvertTrackPointsToLocalClusters(AliTrackPointArray &arrPoints,
 212                                                                AliToyMCTrack &tr, Double_t t0, Int_t type)
 213 {
 214   //
 215   //
 216   //
 217
 218
 219   const Int_t npoints=Int_t(arrPoints.GetUniqueID());
 220   Int_t secOld=-1;
 221
 222   // create an array for graphs which are used for local interpolation
 223   // we need a new graph if the sector changed since then the local frame will also change
 224   TObjArray arrGraphsXY(72);
 225   arrGraphsXY.SetOwner();
 226   TObjArray arrGraphsXZ(72);
 227   arrGraphsXZ.SetOwner();
 228   AliTrackPoint p;
 229   //create initial graph
 230   TGraph *grXY=0x0;
 231   TGraph *grXZ=0x0;
 232   //row -> sector mapping
 233   Int_t rowMap[159];
 234   for (Int_t irow=0; irow<159; ++irow) rowMap[irow]=-1;
 235
 236   // 1. Step
 237   // Make from the list of global space points graphs in the local frame
 238   // one graph per sector is needed
 239   for (Int_t ipoint=0; ipoint<npoints; ++ipoint){
 240     arrPoints.GetPoint(p,ipoint);
 241     Float_t xyz[3] = {p.GetX(),p.GetY(),p.GetZ()};
 242     Int_t index[3] = {0,0,0};
 243     Int_t row = fTPCParam->GetPadRow(xyz,index);
 244     // rotate space point to local frame
 245     // the angle is given by the VolumeID which was set in CrateSpacePoints
 246     const Int_t sec=p.GetVolumeID();
 247     if (row<0 || (sec<36 && row>62) || row>95 ) continue;
 248     Double_t angle=((sec%18)*20.+10.)/TMath::RadToDeg();
 249     AliTrackPoint pRot=p.Rotate(angle);
 250     const Int_t secrow=row+(sec>35)*63;
 251     if (rowMap[secrow]==-1) rowMap[secrow]=sec;
 252     // check if we need a new graph (sector change)
 253     if (secOld!=sec){
 254       grXY=new TGraph;
 255       grXZ=new TGraph;
 256       arrGraphsXY.AddAt(grXY,sec);
 257       arrGraphsXZ.AddAt(grXZ,sec);
 258     }
 259
 260     //add coordinates in local frame for later interpolation
 261     grXY->SetPoint(grXY->GetN(), pRot.GetX(), pRot.GetY());
 262     grXZ->SetPoint(grXZ->GetN(), pRot.GetX(), pRot.GetZ());
 263     secOld=sec;
 264   }
 265
 266   // 2. Step
 267   // create in the center of each row a space point by using the graph to interpolate
 268   // the the center of the row. This is done in xy and xz
 269   TSpline3 *splXY=0x0;
 270   TSpline3 *splXZ=0x0;
 271   AliTPCclusterMI tempCl;
 272   secOld=-1;
 273   for (Int_t irow=0; irow<159; ++irow ){
 274     const Int_t sec     = rowMap[irow];
 275     if (sec==-1) continue;
 276     const Int_t secrow = irow<63?irow:irow-63;
 277     Double_t localX = fTPCParam->GetPadRowRadii(sec,secrow);
 278     // get graph for the current row
 279     if (sec!=secOld){
 280       delete splXY;
 281       splXY=0x0;
 282       delete splXZ;
 283       splXZ=0x0;
 284
 285       grXY=(TGraph*)arrGraphsXY.At(sec);
 286       grXZ=(TGraph*)arrGraphsXZ.At(sec);
 287       if (!grXY) continue;
 288
 289 //       if(grXY->GetN()>1 && grXZ->GetN()>1) { //causes segmentation violation if N==1
 290 //              splXY=new TSpline3("splXY",grXY);
 291 //      splXZ=new TSpline3("splXZ",grXZ);
 292 //       }
 293 //       else {
 294         //TODO: make a cluster also in the sector w only one space point?
 295 //      continue;
 296         // Double_t tempX=0., tempY = 0., tempZ = 0.;
 297
 298         // grXY->GetPoint(0,tempX,localY);
 299         // grXZ->GetPoint(0,tempX,localZ);
 300 //       }
 301
 302     }
 303     secOld=sec;
 304
 305     // check we are in an active area
 306 //     if (splXY->FindX(localX)<1 || splXZ->FindX(localX)<1) continue;
 307     if ( localX<grXY->GetX()[0] || localX>grXY->GetX()[grXY->GetN()-1] || localX<grXZ->GetX()[0] || localX>grXZ->GetX()[grXZ->GetN()-1]) continue;
 308
 309     //get interpolated value at the center for the pad row
 310     //  using splines
 311 //     const Double_t localY=splXY->Eval(localX/*,0x0,"S"*/);
 312 //     const Double_t localZ=splXZ->Eval(localX/*,0x0,"S"*/);
 313     const Double_t localY=grXY->Eval(localX/*,0x0,"S"*/);
 314     const Double_t localZ=grXZ->Eval(localX/*,0x0,"S"*/);
 315     Float_t xyz[3]={localX,localY,localZ};
 316
 317     if (!SetupCluster(tempCl,xyz,sec,t0)) continue;
 318     tempCl.SetLabel(tr.GetUniqueID(), 0);
 319
 320     if (type==0) tr.AddSpacePoint(tempCl);
 321     else tr.AddDistortedSpacePoint(tempCl);
 322 //     printf("SetupCluster %3d: (%.2f, %.2f, %.2f), %d, %.2f\n",irow,xyz[0],xyz[1],xyz[2],sec,t0);
 323   }
 324
 325   delete splXY;
 326   splXY=0x0;
 327   delete splXZ;
 328   splXZ=0x0;
 329
 330 }
 331
 332 //________________________________________________________________
 333 Bool_t AliToyMCEventGenerator::SetupCluster(AliTPCclusterMI &tempCl, Float_t xyz[3], Int_t sec, Double_t t0)
 334 {
 335   //
 336   //
 337   //
 338
 339   const Double_t kSigmaY   = 0.1;
 340   const Double_t kSigmaZ   = 0.1;
 341   const Double_t kMaxZ0    = fTPCParam->GetZLength();
 342   //TODO: Get this from the OCDB at some point?
 343   const Double_t kDriftVel = fTPCParam->GetDriftV();
 344
 345   tempCl.SetX(xyz[0]);
 346   tempCl.SetY(xyz[1]);
 347   tempCl.SetZ(xyz[2]);
 348
 349   tempCl.SetSigmaY2(kSigmaY*kSigmaY);
 350   tempCl.SetSigmaZ2(kSigmaZ*kSigmaZ);
 351
 352   // transform from the local coordinates to the coordinates expressed in pad coordinates
 353   Int_t index[3] = {0,sec,0};
 354   fTPCParam->Transform2to3(xyz,index);
 355   fTPCParam->Transform3to4(xyz,index);
 356   fTPCParam->Transform4to8(xyz,index);
 357
 358   const Int_t   row   = index[2];
 359   const Int_t   nPads = fTPCParam->GetNPads(sec, row);
 360   // pad is fractional, but it needs to be shifted from the center
 361   // to the edge of the row
 362   const Float_t pad   = xyz[1] + nPads/2;
 363
 364   tempCl.SetRow(row);
 365   tempCl.SetPad(pad);
 366   Float_t timeBin=Float_t(t0 + (kMaxZ0-TMath::Abs(tempCl.GetZ()))/kDriftVel);
 367   tempCl.SetTimeBin(timeBin); // set time as t0  + drift time from dist z
 368   tempCl.SetDetector(sec);
 369
 370   //check if we are in the active area
 371   if (pad<0 || pad>=nPads) return kFALSE;
 372
 373   return kTRUE;
 374 }
 375
 376 //________________________________________________________________
 377 Bool_t AliToyMCEventGenerator::ConnectOutputFile()
 378 {
 379   //
 380   // Create the output file name and tree and connect the event
 381   //
 382
 383   fOutFile = new TFile(fOutputFileName.Data(),"recreate");
 384
 385   if (!fOutFile || !fOutFile->IsOpen()){
 386     delete fOutFile;
 387     fOutFile=0x0;
 388     return kFALSE;
 389   }
 390
 391   fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
 392   fOutTree->Branch("event","AliToyMCEvent",&fEvent);
 393
 394   return kTRUE;
 395 }
 396
 397 //________________________________________________________________
 398 Bool_t AliToyMCEventGenerator::CloseOutputFile()
 399 {
 400   //
 401   // close the output file
 402   //
 403   if (!fOutFile) return kFALSE;
 404   fOutFile->Write();
 405   fOutFile->Close();
 406   delete fOutFile;
 407   fOutFile=0x0;
 408
 409   return kTRUE;
 410 }
 411
 412 //________________________________________________________________
 413 void AliToyMCEventGenerator::FillTree()
 414 {
 415   // fill the tree
 416   if (fOutTree) fOutTree->Fill();
 417 }
 418
 419 //________________________________________________________________
 420 void AliToyMCEventGenerator::SetSpaceCharge(EEpsilon epsilon, EGasType gasType/*=kNeCO2_9010*/, ECollRate collRate/*=k50kHz*/)
 421 {
 422   //
 423   // Set the space charge conditions
 424   //
 425   fSpaceChargeFile="$ALICE_ROOT/TPC/Calib/maps/SC";
 426   switch (gasType) {
 427     case kNeCO2_9010:
 428       fSpaceChargeFile.Append("_NeCO2");
 429       break;
 430   }
 431   switch (epsilon) {
 432     case kEps5:
 433       fSpaceChargeFile.Append("_eps5");
 434       break;
 435     case kEps10:
 436       fSpaceChargeFile.Append("_eps10");
 437       break;
 438     case kEps20:
 439       fSpaceChargeFile.Append("_eps20");
 440       break;
 441   }
 442   switch (collRate) {
 443     case k50kHz:
 444       fSpaceChargeFile.Append("_50kHz");
 445       break;
 446   }
 447   fSpaceChargeFile.Append("_precal.root");
 448 }
 449
 450 //________________________________________________________________
 451 void AliToyMCEventGenerator::InitSpaceCharge()
 452 {
 453   //
 454   // init the space charge conditions
 455   // this should be called after the tree was connected
 456   //
 457
 458   AliInfo(Form("Using space charge map file: '%s'",fSpaceChargeFile.Data()));
 459
 460   TFile f(fSpaceChargeFile.Data());
 461   fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
 462
 463   if (fOutTree){
 464     AliInfo("Attaching space charge map file name to the tree");
 465     fOutTree->GetUserInfo()->Add(new TObjString(fSpaceChargeFile.Data()));
 466   }
 467 }