Fixed ClassDef() accidentally removed
[u/mrichter/AliRoot.git] / TPC / TPCupgrade / 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