o add possibility to switch off Material budget correction
[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
b0c26b55 8#include <AliLog.h>
a1a695e5 9#include <AliTPCROC.h>
10#include <AliTrackPointArray.h>
526ddf0e 11#include <AliTrackerBase.h>
12#include <AliCDBManager.h>
13#include <AliTPCParam.h>
14#include <AliGeomManager.h>
15#include <AliTPCcalibDB.h>
526ddf0e 16#include <AliTPCclusterMI.h>
a1a695e5 17#include <AliTPCSpaceCharge3D.h>
18
19#include "AliToyMCEvent.h"
20#include "AliToyMCTrack.h"
21
22#include "AliToyMCEventGenerator.h"
23
de0014b7 24ClassImp(AliToyMCEventGenerator);
526ddf0e 25
26
de0014b7 27AliToyMCEventGenerator::AliToyMCEventGenerator()
526ddf0e 28 :TObject()
29 ,fTPCParam(0x0)
a1a695e5 30 ,fEvent(0x0)
31 ,fSpaceCharge(0x0)
32 ,fOutputFileName("toyMC.root")
33 ,fOutFile(0x0)
34 ,fOutTree(0x0)
526ddf0e 35{
36 //TODO: Add method to set custom space charge file
37 fSpaceCharge = new AliTPCSpaceCharge3D();
a1a695e5 38 fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz.root");
526ddf0e 39 fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
40 //fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
41 fSpaceCharge->InitSpaceCharge3DDistortion();
b0c26b55 42// fSpaceCharge->CreateHistoSCinZR(0.,50,50)->Draw("surf1");
43// fSpaceCharge->CreateHistoDRPhiinZR(0,100,100)->Draw("colz");
a1a695e5 44 //!!! This should be handled by the CongiOCDB macro
efef400b 45 // const char* ocdb="local://$ALICE_ROOT/OCDB/";
46 // AliCDBManager::Instance()->SetDefaultStorage(ocdb);
47 // AliCDBManager::Instance()->SetRun(0);
48 // TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
49 // AliGeomManager::LoadGeometry("");
526ddf0e 50 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
efef400b 51 //std::cout<<"----------->"<<fTPCParam<<std::endl;
526ddf0e 52 fTPCParam->ReadGeoMatrices();
efef400b 53 //std::cout<<"----------->"<<fTPCParam<<std::endl;
54
526ddf0e 55}
56//________________________________________________________________
de0014b7 57AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
526ddf0e 58 :TObject(gen)
a1a695e5 59 ,fTPCParam(gen.fTPCParam)
60 ,fEvent(0x0)
526ddf0e 61 ,fSpaceCharge(gen.fSpaceCharge)
a1a695e5 62 ,fOutputFileName(gen.fOutputFileName)
63 ,fOutFile(0x0)
64 ,fOutTree(0x0)
526ddf0e 65{
66 //
67}
68//________________________________________________________________
de0014b7 69AliToyMCEventGenerator::~AliToyMCEventGenerator()
526ddf0e 70{
71 delete fSpaceCharge;
a1a695e5 72 //!!! Is fTPCParam not owned by the CDB manager?
526ddf0e 73 delete fTPCParam;
74}
75//________________________________________________________________
a1a695e5 76Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
77{
78 //
79 //
80 //
81
82 //!!! Should this complete function not be moved to AliToyMCTrack, and directly called from the constructor
83 // or setter of the track parameters?
526ddf0e 84
85 if(!fTPCParam) {
86 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
87 fTPCParam->ReadGeoMatrices();
88 std::cout << "init tpc param" << std::endl;
89 }
90 const Int_t kNRows=AliTPCROC::Instance()->GetNRows(0)+AliTPCROC::Instance()->GetNRows(36);
91 // const Double_t kRTPC0 =AliTPCROC::Instance()->GetPadRowRadii(0,0);
92 // const Double_t kRTPC1 =AliTPCROC::Instance()->GetPadRowRadii(36,AliTPCROC::Instance()->GetNRows(36)-1);
93 const Double_t kMaxSnp = 0.85;
b0c26b55 94 const Double_t kSigmaY=0.1;
526ddf0e 95 const Double_t kSigmaZ=0.1;
b0c26b55 96 //!!! why divide by 1000000? This will be in cm/us then. But I guess everything should be in sec;
97 // const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
98 const Double_t kDriftVel = fTPCParam->GetDriftV();
526ddf0e 99 const Double_t kMaxZ0=fTPCParam->GetZLength();
100 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
101
102 const Double_t iFCRadius= 83.5; //radius constants found in AliTPCCorrection.cxx
103 const Double_t oFCRadius= 254.5;
104
de0014b7 105 AliToyMCTrack track(trackIn);
526ddf0e 106 const Int_t nMaxPoints = kNRows+50;
107 AliTrackPointArray pointArray0(nMaxPoints);
108 AliTrackPointArray pointArray1(nMaxPoints);
109 Double_t xyz[3];
110
a1a695e5 111 //!!! when does the propagation not work, how often does it happen?
b0c26b55 112 if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp)) {
113 AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
114 return kFALSE;
115 }
526ddf0e 116
117 Int_t npoints=0;
118 Float_t covPoint[6]={0,0,0, kSigmaY*kSigmaY,0,kSigmaZ*kSigmaZ}; //covariance at the local frame
119
120
121 //Simulate track from inner field cage radius to outer and save points along trajectory
b0c26b55 122 Int_t nMinPoints = 40;
123 //!!! changed to loop over pad rows
124 // an increment of 1cm only might loose pad rows in IROC: 4x7.5 mm2
125 // or we change this to a smaller increment
126 // Do we need this initial creation of space points at all? can't we
127 // simply loop over all rows, crate the clusters in the center of the row,
128 // distort them and fill the distorted clusters
129
130 //for (Double_t radius=iFCRadius; radius<oFCRadius; radius++){
131 for (Int_t irow=0; irow<159; ++irow ){
132 const Int_t isec = irow<63?0:36;
133 const Int_t isecrow = irow<63?irow:irow-63;
134 Double_t radius = fTPCParam->GetPadRowRadii(isec,isecrow);
526ddf0e 135
b0c26b55 136 //!!! changed from return 0 to continue -> Please check
137 if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,5,kTRUE,kMaxSnp)) {
138 AliError(Form("Propagation to %d: %.2f failed\n",irow,radius));
139 continue;
140 }
526ddf0e 141 track.GetXYZ(xyz);
b0c26b55 142
143 //!!! Why is this smeared
526ddf0e 144 xyz[0]+=gRandom->Gaus(0,0.000005);
145 xyz[1]+=gRandom->Gaus(0,0.000005);
146 xyz[2]+=gRandom->Gaus(0,0.000005);
147
148 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
b0c26b55 149 //!!! for future we might want to incude space points outside the active area
150 // since due to the distotions they might get inside
151 // in this case it will of course make sense to make a dense
152 // space point creation first and afterwards calculate the cluster
153 // but then we need to treat cluster and distorted cluster symmetrically (see below)
526ddf0e 154 if (TMath::Abs(track.GetX())<iFCRadius) continue;
155 if (TMath::Abs(track.GetX())>oFCRadius) continue;
156
157 AliTrackPoint pIn0; // space point
158 AliTrackPoint pIn1;
159 Int_t sector= (xyz[2]>0)? 0:18;
160 pointArray0.GetPoint(pIn0,npoints);
161 pointArray1.GetPoint(pIn1,npoints);
162 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
163 Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
164 fSpaceCharge->DistortPoint(distPoint, sector);
165 pIn0.SetXYZ(xyz[0], xyz[1],xyz[2]);
166 pIn1.SetXYZ(distPoint[0], distPoint[1],distPoint[2]);
167 track.Rotate(alpha);
168 AliTrackPoint prot0 = pIn0.Rotate(alpha); // rotate to the local frame - non distoted point
169 AliTrackPoint prot1 = pIn1.Rotate(alpha); // rotate to the local frame - distorted point
170 prot0.SetXYZ(prot0.GetX(),prot0.GetY(), prot0.GetZ(),covPoint);
171 prot1.SetXYZ(prot1.GetX(),prot1.GetY(), prot1.GetZ(),covPoint);
172 pIn0=prot0.Rotate(-alpha); // rotate back to global frame
173 pIn1=prot1.Rotate(-alpha); // rotate back to global frame
174 pointArray0.AddPoint(npoints, &pIn0);
175 pointArray1.AddPoint(npoints, &pIn1);
176 npoints++;
177 if (npoints>=nMaxPoints) break;
178 }
179
180 if (npoints<nMinPoints) return 0;
181
182 //save space points and make clusters of distorted points
b0c26b55 183
184 //!!! Here both, clusters and distorted clusters should be in the center of the pad row
185 // So in case the distorted clusters are calculated as the COG of all space points
186 // in a row, the same should be done for the undistorted clusters! There should
187 // only be one cluster per track per row
526ddf0e 188 Int_t lastRow = 0;
189 Int_t pntsInCurrentRow = 0;
190 Double_t xt=0, yt=0, zt=0;
191 for(Int_t iPoint = 0; iPoint<npoints; iPoint++){
b0c26b55 192
193 //undistorted cluster, currently we assume one space point per pad row, see comments above
526ddf0e 194 AliTPCclusterMI *tempCl = trackIn.AddSpacePoint(AliTPCclusterMI());
195
196 const Float_t *xArr = pointArray0.GetX();
197 const Float_t *yArr = pointArray0.GetY();
198 const Float_t *zArr = pointArray0.GetZ();
199
526ddf0e 200 tempCl->SetX(xArr[iPoint]);
201 tempCl->SetY(yArr[iPoint]);
202 tempCl->SetZ(zArr[iPoint]);
203
b0c26b55 204 Float_t xyzUC[3] = {xArr[iPoint],yArr[iPoint],zArr[iPoint]};
205 Int_t indexUC[3] = {0};
206
207 Float_t padRowUC = fTPCParam->GetPadRow(xyzUC,indexUC);
208 Int_t sectorUC = indexUC[1];
209
210 // use pad row = 255 if outside active area
211 if(padRowUC < 0 || (sectorUC < 36 && padRowUC >62) || (sectorUC > 35 &&padRowUC > 95) ) padRowUC=255; //outside sensitive area
212
213 tempCl->SetPad(xyzUC[1]);
214 tempCl->SetRow(padRowUC);
215 tempCl->SetDetector(sectorUC);
216 tempCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(tempCl->GetZ()))/kDriftVel); // set time as t0 + drift time from dist z
217
218 // Distorted cluster
219 const Float_t *xArrDist = pointArray1.GetX();
220 const Float_t *yArrDist = pointArray1.GetY();
221 const Float_t *zArrDist = pointArray1.GetZ();
222
a1a695e5 223 Float_t xyz2[3] = {xArrDist[iPoint],yArrDist[iPoint],zArrDist[iPoint]};
526ddf0e 224 Int_t index[3] = {0};
225
a1a695e5 226 Float_t padRow = fTPCParam->GetPadRow(xyz2,index);
526ddf0e 227 Int_t sector = index[1];
228
229 if(padRow < 0 || (sector < 36 && padRow>62) || (sector > 35 &&padRow > 95) ) continue; //outside sensitive area
230
b0c26b55 231 //!!! could it be that the last row is missed?
232 // assume row=95 and there is only 1 space point then this part will not be reached I think
526ddf0e 233 if(lastRow!=padRow) {
234
235
236 //make and save distorted cluster
237 if(pntsInCurrentRow>0){
238 xt/=pntsInCurrentRow;
239 yt/=pntsInCurrentRow;
240 zt/=pntsInCurrentRow;
241 AliTPCclusterMI *tempDistCl = trackIn.AddDistortedSpacePoint(AliTPCclusterMI());
242 tempDistCl->SetX(xt);
243 tempDistCl->SetY(yt);
244 tempDistCl->SetZ(zt);
245 Float_t clxyz[3] = {xt,yt,zt};
246 Int_t clindex[3] = {0};
247 Int_t clRow = fTPCParam->GetPadRow(clxyz,clindex);
b0c26b55 248// Int_t nPads = fTPCParam->GetNPads(clindex[1], clRow);
249// Int_t pad = TMath::Nint(clxyz[1] + nPads/2); //taken from AliTPC.cxx
250 //!!! pad should be fractional
251// tempDistCl->SetPad(pad);
252 tempDistCl->SetPad(clxyz[1]);
253 tempDistCl->SetRow(clRow);
526ddf0e 254 tempDistCl->SetDetector(clindex[1]);
255 tempDistCl->SetTimeBin(t0 + (kMaxZ0-TMath::Abs(zt))/kDriftVel); // set time as t0 + drift time from dist z coordinate to pad plane
256 }
257 xt = 0;
258 yt = 0;
259 zt = 0;
260 pntsInCurrentRow = 0;
261 }
262
263 xt+=xArrDist[iPoint];
264 yt+=yArrDist[iPoint];
265 zt+=zArrDist[iPoint];
266 pntsInCurrentRow++;
267 lastRow = padRow;
268 }
269
270 return 1;
271
272
273
274}
a1a695e5 275//________________________________________________________________
276Bool_t AliToyMCEventGenerator::ConnectOutputFile()
277{
278 //
279 // Create the output file name and tree and connect the event
280 //
281
282 fOutFile = new TFile(fOutputFileName.Data(),"recreate");
283
284 if (!fOutFile || !fOutFile->IsOpen()){
285 delete fOutFile;
286 fOutFile=0x0;
287 return kFALSE;
288 }
289
290 fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
291 fOutTree->Branch("event","AliToyMCEvent",&fEvent);
292
293 return kTRUE;
294}
295
296//________________________________________________________________
297Bool_t AliToyMCEventGenerator::CloseOutputFile()
298{
299 //
300 // close the output file
301 //
302 if (!fOutFile) return kFALSE;
303 fOutFile->Write();
304 fOutFile->Close();
305 delete fOutFile;
306 fOutFile=0x0;
307
308 return kTRUE;
309}
310
311//________________________________________________________________
312void AliToyMCEventGenerator::FillTree()
313{
314 // fill the tree
315 if (fOutTree) fOutTree->Fill();
316}
317