o add space charge maps needed for toy event generator
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGenerator.cxx
CommitLineData
526ddf0e 1#include <iostream>
de0014b7 2#include "AliToyMCEventGenerator.h"
526ddf0e 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>
de0014b7 15ClassImp(AliToyMCEventGenerator);
526ddf0e 16
17
de0014b7 18AliToyMCEventGenerator::AliToyMCEventGenerator()
526ddf0e 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//________________________________________________________________
de0014b7 39AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
526ddf0e 40 :TObject(gen)
41 ,fSpaceCharge(gen.fSpaceCharge)
42{
43 //
44}
45//________________________________________________________________
de0014b7 46AliToyMCEventGenerator::~AliToyMCEventGenerator()
526ddf0e 47{
48 delete fSpaceCharge;
49 delete fTPCParam;
50}
51//________________________________________________________________
de0014b7 52Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0) {
526ddf0e 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
de0014b7 72 AliToyMCTrack track(trackIn);
526ddf0e 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}