o Minor code cleanup
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGenerator.cxx
CommitLineData
526ddf0e 1#include <iostream>
a1a695e5 2
526ddf0e 3#include <TDatabasePDG.h>
4#include <TRandom.h>
5#include <TH2F.h>
a1a695e5 6#include <TGeoGlobalMagField.h>
7
8#include <AliTPCROC.h>
9#include <AliTrackPointArray.h>
526ddf0e 10#include <AliTrackerBase.h>
11#include <AliCDBManager.h>
12#include <AliTPCParam.h>
13#include <AliGeomManager.h>
14#include <AliTPCcalibDB.h>
526ddf0e 15#include <AliTPCclusterMI.h>
a1a695e5 16#include <AliTPCSpaceCharge3D.h>
17
18#include "AliToyMCEvent.h"
19#include "AliToyMCTrack.h"
20
21#include "AliToyMCEventGenerator.h"
22
de0014b7 23ClassImp(AliToyMCEventGenerator);
526ddf0e 24
25
de0014b7 26AliToyMCEventGenerator::AliToyMCEventGenerator()
526ddf0e 27 :TObject()
28 ,fTPCParam(0x0)
a1a695e5 29 ,fEvent(0x0)
30 ,fSpaceCharge(0x0)
31 ,fOutputFileName("toyMC.root")
32 ,fOutFile(0x0)
33 ,fOutTree(0x0)
526ddf0e 34{
35 //TODO: Add method to set custom space charge file
36 fSpaceCharge = new AliTPCSpaceCharge3D();
a1a695e5 37 fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz.root");
526ddf0e 38 fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
39 //fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
40 fSpaceCharge->InitSpaceCharge3DDistortion();
41 fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
a1a695e5 42 fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
43 //!!! This should be handled by the CongiOCDB macro
526ddf0e 44 const char* ocdb="local://$ALICE_ROOT/OCDB/";
45 AliCDBManager::Instance()->SetDefaultStorage(ocdb);
46 AliCDBManager::Instance()->SetRun(0);
47 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
48 AliGeomManager::LoadGeometry("");
49 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
50 fTPCParam->ReadGeoMatrices();
51}
52//________________________________________________________________
de0014b7 53AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
526ddf0e 54 :TObject(gen)
a1a695e5 55 ,fTPCParam(gen.fTPCParam)
56 ,fEvent(0x0)
526ddf0e 57 ,fSpaceCharge(gen.fSpaceCharge)
a1a695e5 58 ,fOutputFileName(gen.fOutputFileName)
59 ,fOutFile(0x0)
60 ,fOutTree(0x0)
526ddf0e 61{
62 //
63}
64//________________________________________________________________
de0014b7 65AliToyMCEventGenerator::~AliToyMCEventGenerator()
526ddf0e 66{
67 delete fSpaceCharge;
a1a695e5 68 //!!! Is fTPCParam not owned by the CDB manager?
526ddf0e 69 delete fTPCParam;
70}
71//________________________________________________________________
a1a695e5 72Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
73{
74 //
75 //
76 //
77
78 //!!! Should this complete function not be moved to AliToyMCTrack, and directly called from the constructor
79 // or setter of the track parameters?
526ddf0e 80
81 if(!fTPCParam) {
82 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
83 fTPCParam->ReadGeoMatrices();
84 std::cout << "init tpc param" << std::endl;
85 }
86 const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36);
87 // const Double_t kRTPC0 =AliTPCROC::Instance()->GetPadRowRadii(0,0);
88 // const Double_t kRTPC1 =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1);
89 const Double_t kMaxSnp = 0.85;
90 const Double_t kSigmaY=0.1;
91 const Double_t kSigmaZ=0.1;
92 const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
93 const Double_t kMaxZ0=fTPCParam->GetZLength();
94 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
95
96 const Double_t iFCRadius= 83.5; //radius constants found in AliTPCCorrection.cxx
97 const Double_t oFCRadius= 254.5;
98
de0014b7 99 AliToyMCTrack track(trackIn);
526ddf0e 100 const Int_t nMaxPoints = kNRows+50;
101 AliTrackPointArray pointArray0(nMaxPoints);
102 AliTrackPointArray pointArray1(nMaxPoints);
103 Double_t xyz[3];
104
a1a695e5 105 //!!! when does the propagation not work, how often does it happen?
526ddf0e 106 if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) return 0;
107
108 Int_t npoints=0;
109 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
110
111
112 //Simulate track from inner field cage radius to outer and save points along trajectory
113 Int_t nMinPoints = 40;
114 for (Double_t radius=iFCRadius; radius<oFCRadius; radius++){
115
116 if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) return 0;
117 track.GetXYZ(xyz);
118 xyz[0]+=gRandom->Gaus(0,0.000005);
119 xyz[1]+=gRandom->Gaus(0,0.000005);
120 xyz[2]+=gRandom->Gaus(0,0.000005);
121
122 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
123 if (TMath::Abs(track.GetX())<iFCRadius) continue;
124 if (TMath::Abs(track.GetX())>oFCRadius) continue;
125
126 AliTrackPoint pIn0; // space point
127 AliTrackPoint pIn1;
128 Int_t sector= (xyz[2]>0)? 0:18;
129 pointArray0.GetPoint(pIn0,npoints);
130 pointArray1.GetPoint(pIn1,npoints);
131 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
132 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
133 fSpaceCharge->DistortPoint(distPoint, sector);
134 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
135 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
136 track.Rotate(alpha);
137 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
138 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
139 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
140 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
141 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
142 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
143 pointArray0.AddPoint(npoints, &pIn0);
144 pointArray1.AddPoint(npoints, &pIn1);
145 npoints++;
146 if (npoints>=nMaxPoints) break;
147 }
148
149 if (npoints<nMinPoints) return 0;
150
151 //save space points and make clusters of distorted points
152 Int_t lastRow = 0;
153 Int_t pntsInCurrentRow = 0;
154 Double_t xt=0, yt=0, zt=0;
155 for(Int_t iPoint = 0; iPoint<npoints; iPoint++){
156
157 AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());
158
159 const Float_t *xArr = pointArray0.GetX();
160 const Float_t *yArr = pointArray0.GetY();
161 const Float_t *zArr = pointArray0.GetZ();
162
163 const Float_t *xArrDist = pointArray1.GetX();
164 const Float_t *yArrDist = pointArray1.GetY();
165 const Float_t *zArrDist = pointArray1.GetZ();
166
167 tempCl->SetX(xArr[iPoint]);
168 tempCl->SetY(yArr[iPoint]);
169 tempCl->SetZ(zArr[iPoint]);
170
a1a695e5 171 Float_t xyz2[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]};
526ddf0e 172 Int_t index[3] = {0};
173
a1a695e5 174 Float_t padRow = fTPCParam->GetPadRow(xyz2,index);
526ddf0e 175 Int_t sector = index[1];
176
177 if(padRow < 0 || (sector < 36 && padRow>62) || (sector > 35 &&padRow > 95) ) continue; //outside sensitive area
178
179
180 if(lastRow!=padRow) {
181
182
183 //make and save distorted cluster
184 if(pntsInCurrentRow>0){
185 xt/=pntsInCurrentRow;
186 yt/=pntsInCurrentRow;
187 zt/=pntsInCurrentRow;
188 AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI());
189 tempDistCl->SetX(xt);
190 tempDistCl->SetY(yt);
191 tempDistCl->SetZ(zt);
192 Float_t clxyz[3] = {xt,yt,zt};
193 Int_t clindex[3] = {0};
194 Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex);
195 Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow);
196 Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx
197 tempDistCl->SetPad(pad);
198 tempDistCl->SetRow(clRow);
199 tempDistCl->SetDetector(clindex[1]);
200 tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0 + drift time from dist z coordinate to pad plane
201 }
202 xt = 0;
203 yt = 0;
204 zt = 0;
205 pntsInCurrentRow = 0;
206 }
207
208 xt+=xArrDist[iPoint];
209 yt+=yArrDist[iPoint];
210 zt+=zArrDist[iPoint];
211 pntsInCurrentRow++;
212 lastRow = padRow;
213 }
214
215 return 1;
216
217
218
219}
a1a695e5 220//________________________________________________________________
221Bool_t AliToyMCEventGenerator::ConnectOutputFile()
222{
223 //
224 // Create the output file name and tree and connect the event
225 //
226
227 fOutFile = new TFile(fOutputFileName.Data(),"recreate");
228
229 if (!fOutFile || !fOutFile->IsOpen()){
230 delete fOutFile;
231 fOutFile=0x0;
232 return kFALSE;
233 }
234
235 fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
236 fOutTree->Branch("event","AliToyMCEvent",&fEvent);
237
238 return kTRUE;
239}
240
241//________________________________________________________________
242Bool_t AliToyMCEventGenerator::CloseOutputFile()
243{
244 //
245 // close the output file
246 //
247 if (!fOutFile) return kFALSE;
248 fOutFile->Write();
249 fOutFile->Close();
250 delete fOutFile;
251 fOutFile=0x0;
252
253 return kTRUE;
254}
255
256//________________________________________________________________
257void AliToyMCEventGenerator::FillTree()
258{
259 // fill the tree
260 if (fOutTree) fOutTree->Fill();
261}
262