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