TPC/Upgrade/ToyMCEventGenerator.cxx

   1 #include <iostream>
   2 #include "ToyMCEventGenerator.h"
   3 #include <AliTPCROC.h>
   4 #include <AliTrackPointArray.h>
   5 #include <TDatabasePDG.h>
   6 #include <TRandom.h>
   7 #include <TH2F.h>
   8 #include <AliTrackerBase.h>
   9 #include <AliCDBManager.h>
  10 #include <AliTPCParam.h>
  11 #include <AliGeomManager.h>
  12 #include <AliTPCcalibDB.h>
  13 #include <TGeoGlobalMagField.h>
  14 #include <AliTPCclusterMI.h>
  15 ClassImp(ToyMCEventGenerator);
  16
  17
  18 ToyMCEventGenerator::ToyMCEventGenerator()
  19   :TObject()
  20   ,fTPCParam(0x0)
  21 {
  22   //TODO: Add method to set custom space charge file
  23   fSpaceCharge = new AliTPCSpaceCharge3D();
  24   fSpaceCharge->SetSCDataFileName("SC_NeCO2_eps5_50kHz.root");
  25   fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
  26   //fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
  27   fSpaceCharge->InitSpaceCharge3DDistortion();
  28   fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
  29   fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
  30   const char* ocdb="local://$ALICE_ROOT/OCDB/";
  31   AliCDBManager::Instance()->SetDefaultStorage(ocdb);
  32   AliCDBManager::Instance()->SetRun(0);
  33   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
  34   AliGeomManager::LoadGeometry("");
  35   fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  36   fTPCParam->ReadGeoMatrices();
  37 }
  38 //________________________________________________________________
  39 ToyMCEventGenerator::ToyMCEventGenerator(const ToyMCEventGenerator &gen)
  40   :TObject(gen)
  41   ,fSpaceCharge(gen.fSpaceCharge)
  42 {
  43   //
  44 }
  45 //________________________________________________________________
  46 ToyMCEventGenerator::~ToyMCEventGenerator()
  47 {
  48   delete fSpaceCharge;
  49   delete fTPCParam;
  50 }
  51 //________________________________________________________________
  52 Bool_t ToyMCEventGenerator::DistortTrack(ToyMCTrack &trackIn, Double_t t0) {
  53
  54   if(!fTPCParam) {
  55     fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  56     fTPCParam->ReadGeoMatrices();
  57     std::cout << "init tpc param" << std::endl;
  58   }
  59   const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36);
  60   // const Double_t kRTPC0  =AliTPCROC::Instance()->GetPadRowRadii(0,0);
  61   // const Double_t kRTPC1  =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1);
  62   const Double_t kMaxSnp = 0.85;
  63   const Double_t kSigmaY=0.1;
  64   const Double_t kSigmaZ=0.1;
  65   const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
  66   const Double_t kMaxZ0=fTPCParam->GetZLength();
  67   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  68
  69   const Double_t iFCRadius=  83.5; //radius constants found in AliTPCCorrection.cxx
  70   const Double_t oFCRadius= 254.5;
  71
  72   ToyMCTrack track(trackIn);
  73   const Int_t nMaxPoints = kNRows+50;
  74   AliTrackPointArray pointArray0(nMaxPoints);
  75   AliTrackPointArray pointArray1(nMaxPoints);
  76   Double_t xyz[3];
  77
  78   if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) return 0;
  79
  80   Int_t npoints=0;
  81   Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ};  //covariance at the local frame
  82
  83
  84   //Simulate track from inner field cage radius to outer and save points along trajectory
  85   Int_t nMinPoints = 40;
  86   for (Double_t radius=iFCRadius; radius<oFCRadius; radius++){
  87
  88     if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
  89     track.GetXYZ(xyz);
  90     xyz[0]+=gRandom->Gaus(0,0.000005);
  91     xyz[1]+=gRandom->Gaus(0,0.000005);
  92     xyz[2]+=gRandom->Gaus(0,0.000005);
  93
  94     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
  95     if (TMath::Abs(track.GetX())<iFCRadius) continue;
  96     if (TMath::Abs(track.GetX())>oFCRadius) continue;
  97
  98     AliTrackPoint pIn0;                               // space point
  99     AliTrackPoint pIn1;
 100     Int_t sector= (xyz[2]>0)? 0:18;
 101     pointArray0.GetPoint(pIn0,npoints);
 102     pointArray1.GetPoint(pIn1,npoints);
 103     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
 104     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
 105     fSpaceCharge->DistortPoint(distPoint, sector);
 106     pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
 107     pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
 108     track.Rotate(alpha);
 109     AliTrackPoint prot0 = pIn0.Rotate(alpha);   // rotate to the local frame - non distoted  point
 110     AliTrackPoint prot1 = pIn1.Rotate(alpha);   // rotate to the local frame -     distorted point
 111     prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
 112     prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
 113     pIn0=prot0.Rotate(-alpha);                       // rotate back to global frame
 114     pIn1=prot1.Rotate(-alpha);                       // rotate back to global frame
 115     pointArray0.AddPoint(npoints, &pIn0);
 116     pointArray1.AddPoint(npoints, &pIn1);
 117     npoints++;
 118     if (npoints>=nMaxPoints) break;
 119   }
 120
 121   if (npoints<nMinPoints) return 0;
 122
 123   //save space points and make clusters of distorted points
 124   Int_t lastRow = 0;
 125   Int_t pntsInCurrentRow = 0;
 126   Double_t xt=0, yt=0, zt=0;
 127   for(Int_t iPoint = 0; iPoint<npoints; iPoint++){
 128
 129     AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());
 130
 131     const Float_t *xArr = pointArray0.GetX();
 132     const Float_t *yArr = pointArray0.GetY();
 133     const Float_t *zArr = pointArray0.GetZ();
 134
 135     const Float_t *xArrDist = pointArray1.GetX();
 136     const Float_t *yArrDist = pointArray1.GetY();
 137     const Float_t *zArrDist = pointArray1.GetZ();
 138
 139     tempCl->SetX(xArr[iPoint]);
 140     tempCl->SetY(yArr[iPoint]);
 141     tempCl->SetZ(zArr[iPoint]);
 142
 143     Float_t xyz[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]};
 144     Int_t index[3] = {0};
 145
 146     Float_t padRow = fTPCParam->GetPadRow(xyz,index);
 147     Int_t sector = index[1];
 148
 149     if(padRow < 0 || (sector < 36 && padRow>62) || (sector > 35 &&padRow > 95) ) continue; //outside sensitive area
 150
 151
 152     if(lastRow!=padRow) {
 153
 154
 155       //make and save distorted cluster
 156       if(pntsInCurrentRow>0){
 157         xt/=pntsInCurrentRow;
 158         yt/=pntsInCurrentRow;
 159         zt/=pntsInCurrentRow;
 160         AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI());
 161         tempDistCl->SetX(xt);
 162         tempDistCl->SetY(yt);
 163         tempDistCl->SetZ(zt);
 164         Float_t clxyz[3] = {xt,yt,zt};
 165         Int_t clindex[3] = {0};
 166         Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex);
 167         Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow);
 168         Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx
 169         tempDistCl->SetPad(pad);
 170         tempDistCl->SetRow(clRow);
 171         tempDistCl->SetDetector(clindex[1]);
 172         tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0  + drift time from dist z coordinate to pad plane
 173       }
 174       xt = 0;
 175       yt = 0;
 176       zt = 0;
 177       pntsInCurrentRow = 0;
 178     }
 179
 180     xt+=xArrDist[iPoint];
 181     yt+=yArrDist[iPoint];
 182     zt+=zArrDist[iPoint];
 183     pntsInCurrentRow++;
 184     lastRow = padRow;
 185   }
 186
 187   return 1;
 188
 189
 190
 191 }