o mods
[u/mrichter/AliRoot.git] / TPC / Upgrade / 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 {
54   fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
55   fTPCParam->ReadGeoMatrices();
56   gRandom->SetSeed();
57 }
58 //________________________________________________________________
59 AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
60   :TObject(gen)
61   ,fTPCParam(gen.fTPCParam)
62   ,fEvent(0x0)
63   ,fCurrentTrack(0)
64   ,fTPCCorrection(gen.fTPCCorrection)
65   ,fTPCCorrectionAv(0x0)
66   ,fSCList(gen.fSCList)
67   ,fSCListFile(gen.fSCListFile)
68   ,fCorrectionFile(gen.fCorrectionFile)
69   ,fOutputFileName(gen.fOutputFileName)
70   ,fOutFile(0x0)
71   ,fOutTree(0x0)
72   ,fUseStepCorrection(gen.fUseStepCorrection)
73   ,fUseMaterialBudget(gen.fUseMaterialBudget)
74   ,fIsLaser(gen.fIsLaser)
75   ,fPrereadSCList(gen.fPrereadSCList)
76 {
77   //
78   gRandom->SetSeed();
79 }
80 //________________________________________________________________
81 AliToyMCEventGenerator::~AliToyMCEventGenerator() 
82 {
83   if (HasSCList() &&!fPrereadSCList) delete fTPCCorrection;
84   delete fSCList;
85 }
86
87 //________________________________________________________________
88 Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
89 {
90   //
91   //
92   //
93
94   if(!fTPCParam) {
95     fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
96     fTPCParam->ReadGeoMatrices();
97   }
98
99   MakeITSClusters(trackIn/*,t0*/);
100   MakeTPCClusters(trackIn, t0);
101   MakeTRDClusters(trackIn/*,t0*/);
102
103   return kTRUE;
104 }
105 //________________________________________________________________
106 void AliToyMCEventGenerator::MakeITSClusters(AliToyMCTrack &trackIn/*, Double_t t0*/)
107 {
108   //Upgrade ITS parameters
109   const Int_t nITSLayers = 7;
110   const Double_t ITSRadii[nITSLayers] = {2.2, 2.8, 3.6, 20.0, 22.0, 41.0, 43.0};
111   const Double_t lengthITS[nITSLayers] = {22.4, 24.2, 26.8, 78.0, 83.6, 142.4, 148.6};
112
113   //resolution of the point is 10um
114   const Float_t sigmaY = 0.001;
115   const Float_t sigmaZ = 0.001;
116
117   AliTrackPoint point;
118   
119   const Double_t kMaxSnp = 0.85;
120 //   const Double_t kMaxZ0  = fTPCParam->GetZLength();
121   const Double_t kMass   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
122   
123   AliExternalTrackParam track(trackIn);
124   Double_t xyz[3]  = {0.,0.,0.};
125       
126   if (!AliTrackerBase::PropagateTrackTo(&track,ITSRadii[0],kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
127     AliError(Form("Propagation to %.2f failed\n",ITSRadii[0]));
128     return;
129   }
130
131   for(Int_t iLayer = 0; iLayer<nITSLayers; iLayer++){
132
133     if (!AliTrackerBase::PropagateTrackTo(&track,ITSRadii[iLayer],kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
134       AliError(Form("Propagation to %.2f failed\n",ITSRadii[iLayer]));
135       continue;
136     }
137     // since I don't know how to extract the volumeid of the ITS layers I use the following strategy:
138     //  - rotate the track to the next integer angle
139     //  - use coordinates there and set as volumeid the integer angle
140     //  - in the reco one can then rotate the cluster to the global frame using the angle stored in the volume id
141     track.GetXYZ(xyz);
142     
143 //     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
144     if (TMath::Abs(track.GetZ())>lengthITS[iLayer]/2) continue;
145
146
147     // smear the ideal positions with the cluster resolution
148     xyz[1]+=gRandom->Gaus(0,sigmaY);
149     xyz[2]+=gRandom->Gaus(0,sigmaZ);
150
151     Float_t xyzf[3]={xyz[0],xyz[1],xyz[2]};
152
153     SetPoint(xyzf,sigmaY,sigmaZ,point);
154     
155     trackIn.AddITSPoint(point)->SetUniqueID(trackIn.GetUniqueID());
156   }
157
158 }
159 //________________________________________________________________
160 void AliToyMCEventGenerator::MakeTRDClusters(AliToyMCTrack &trackIn/*, Double_t t0*/)
161
162 {
163   //Uses current TRD parameters
164   const Int_t nTRDLayers = 6;
165   const Double_t distToMid = 3.2 + 30./2; //dist to middle of drift region (radiator + half drift region)
166   const Double_t TRDRadii[nTRDLayers] = {294.5 + distToMid, 307.1 + distToMid, 319.7 + distToMid, 332.3 + distToMid, 344.9 + distToMid, 357.5 + distToMid};
167   const Double_t lengthTRD[nTRDLayers] = {604.0, 634.0, 656.0, 686.0, 700.0, 700.0};
168   
169   const Float_t sigmaY = 0.06;
170   const Float_t sigmaZ = 0.2;
171
172   AliTrackPoint point;
173   
174   const Double_t kMaxSnp = 0.85;
175   const Double_t kMaxZ0  = fTPCParam->GetZLength();
176   const Double_t kMass   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
177   
178   AliExternalTrackParam track(trackIn);
179   Double_t xyz[3]  = {0.,0.,0.};
180       
181   if (!AliTrackerBase::PropagateTrackTo(&track,TRDRadii[0],kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
182     AliError(Form("Propagation to %.2f failed\n",TRDRadii[0]));
183     return;
184   }
185
186   for(Int_t iLayer = 0; iLayer<nTRDLayers; iLayer++){
187
188     if (!AliTrackerBase::PropagateTrackTo(&track,TRDRadii[iLayer],kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
189       AliError(Form("Propagation to %.2f failed\n",TRDRadii[iLayer]));
190       continue;
191     }
192     track.GetXYZ(xyz);
193     
194     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
195     if (TMath::Abs(track.GetZ())>lengthTRD[iLayer]/2) continue;
196     
197
198     // smear the ideal positions with the cluster resolution
199     xyz[1]+=gRandom->Gaus(0,sigmaY);
200     xyz[2]+=gRandom->Gaus(0,sigmaZ);
201     
202     Float_t xyzf[3]={xyz[0],xyz[1],xyz[2]};
203
204     SetPoint(xyzf,sigmaY,sigmaZ,point);
205
206     trackIn.AddTRDPoint(point)->SetUniqueID(trackIn.GetUniqueID());
207   }
208
209 }
210 //________________________________________________________________
211 void AliToyMCEventGenerator::MakeTPCClusters(AliToyMCTrack &trackIn, Double_t t0)
212 {
213
214   // make it big enough to hold all points
215   // store real number of generated points in the unique id
216   const Int_t nMaxPoints=3000;
217   static AliTrackPointArray pointArray0(nMaxPoints);  //undistorted
218   static AliTrackPointArray pointArray1(nMaxPoints);  //distorted
219   
220   //Create space point of undistorted and distorted clusters along the propagated track trajectory
221   CreateSpacePoints(trackIn,pointArray0,pointArray1);
222   //Convert the space points into clusters in the local frame
223   //for undistorted and distorted clusters using the same function
224   ConvertTrackPointsToLocalClusters(pointArray0,trackIn,t0,0);
225   ConvertTrackPointsToLocalClusters(pointArray1,trackIn,t0,1);
226   
227
228 }
229 //________________________________________________________________
230 void AliToyMCEventGenerator::CreateSpacePoints(AliToyMCTrack &trackIn,
231                                                AliTrackPointArray &arrUdist,
232                                                AliTrackPointArray &arrDist)
233 {
234   //
235   // sample the track from the inner to the outer wall of the TPC
236   // a graph is filled in local coordinates for later
237   //
238   
239   Double_t kMaxSnp = 0.85;
240   if (fIsLaser) kMaxSnp=0.99;
241   const Double_t kMaxZ0  = fTPCParam->GetZLength();
242   const Double_t kMass   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
243   
244   const Double_t iFCRadius =  83.5; //radius constants found in AliTPCCorrection.cxx
245   const Double_t oFCRadius = 254.5;
246
247   const Float_t kSigmaY=0.1;
248   const Float_t kSigmaZ=0.1;
249   
250   AliExternalTrackParam track(trackIn);
251   //!!! TODO: make this adjustable perhaps
252   const Double_t stepSize=0.1;
253   Double_t xyz[3]  = {0.,0.,0.};
254   Float_t  xyzf[3] = {0.,0.,0.};
255   
256   //!!! when does the propagation not work, how often does it happen?
257   if (!AliTrackerBase::PropagateTrackTo(&track,iFCRadius,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget) && !fIsLaser) {
258     AliError(Form("Propagation to IFC: %.2f failed\n",iFCRadius));
259     return;
260   }
261   
262   Int_t npoints=0;
263   
264   for (Double_t radius=iFCRadius; radius<oFCRadius; radius+=stepSize){
265     //!!! changed from return 0 to continue -> Please check
266     if (!AliTrackerBase::PropagateTrackTo(&track,radius,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
267       AliError(Form("Propagation to r=%.2f (snp=%.2f) failed\n",radius,track.GetSnp()));
268       continue;
269     }
270     track.GetXYZ(xyz);
271     
272     //!!! Why is this smeared
273     
274     xyzf[0]=Float_t(xyz[0]);
275     xyzf[1]=Float_t(xyz[1]);
276     xyzf[2]=Float_t(xyz[2]);
277     
278     if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
279     //       if (TMath::Abs(track.GetX())<iFCRadius) continue;
280     //       if (TMath::Abs(track.GetX())>oFCRadius) continue;
281     
282     
283     AliTrackPoint pUdist;                               // undistorted space point
284     AliTrackPoint pDist;                                // distorted space point
285     // Set undistorted point
286     SetPoint(xyzf,kSigmaY,kSigmaZ,pUdist);
287     arrUdist.AddPoint(npoints, &pUdist);
288     Int_t sector=pUdist.GetVolumeID();    
289     
290     // set distorted point
291     Float_t distPoint[3]={xyz[0],xyz[1],xyz[2]};
292     Float_t dxyz[3]={0.,0.,0.};
293     if (!fUseStepCorrection){
294       fTPCCorrection->DistortPoint(distPoint, sector);
295     } else {
296       fTPCCorrection->GetCorrectionIntegralDz(distPoint,sector,dxyz,5);
297       distPoint[0]-=dxyz[0];
298       distPoint[1]-=dxyz[1];
299       distPoint[2]-=dxyz[2];
300     }
301     SetPoint(distPoint,kSigmaY,kSigmaZ,pDist);
302     arrDist.AddPoint(npoints, &pDist);
303     
304     ++npoints;
305   }
306   
307   arrUdist.SetUniqueID(npoints);
308   arrDist.SetUniqueID(npoints);
309 }
310
311 //________________________________________________________________
312 void AliToyMCEventGenerator::SetPoint(Float_t xyz[3], Float_t sigmaY, Float_t sigmaZ, AliTrackPoint &point)
313 {
314   //
315   // make AliTrackPoint out of AliTPCclusterMI
316   //
317   
318   //covariance at the local frame
319   //assume 1mm distortion in y and z
320   Int_t i[3]={0,0,0};
321   Float_t cov[6]={0,0,0, sigmaY*sigmaY,0,sigmaZ*sigmaZ};
322   
323   const Float_t alpha = -TMath::ATan2(xyz[1],xyz[0]);
324   const Float_t sin   = TMath::Sin(alpha), cos = TMath::Cos(alpha);
325   
326   Float_t newcov[6];
327   newcov[0] = cov[0]*cos*cos           + 2*cov[1]*sin*cos         + cov[3]*sin*sin;
328   newcov[1] = cov[1]*(cos*cos-sin*sin) + (cov[3]-cov[0])*sin*cos;
329   newcov[2] = cov[2]*cos               + cov[4]*sin;
330   newcov[3] = cov[0]*sin*sin           - 2*cov[1]*sin*cos          +  cov[3]*cos*cos;
331   newcov[4] = cov[4]*cos               - cov[2]*sin;
332   newcov[5] = cov[5];
333   
334   // voluem ID to add later ....
335   point.SetXYZ(xyz);
336   point.SetCov(newcov);
337   // abuse volume ID for the sector number
338   point.SetVolumeID(fTPCParam->Transform0to1(xyz,i));
339
340   // TODO: Add sampled dE/dx  (use SetCharge)
341 }
342
343 //________________________________________________________________
344 void AliToyMCEventGenerator::ConvertTrackPointsToLocalClusters(AliTrackPointArray &arrPoints,
345                                                                AliToyMCTrack &tr, Double_t t0, Int_t type)
346 {
347   //
348   //
349   //
350
351   const Int_t npoints=Int_t(arrPoints.GetUniqueID());
352   Int_t secOld=-1;
353
354   // create an array for graphs which are used for local interpolation
355   // we need a new graph if the sector changed since then the local frame will also change
356   TObjArray arrGraphsXY(72);
357   arrGraphsXY.SetOwner();
358   TObjArray arrGraphsXZ(72);
359   arrGraphsXZ.SetOwner();
360   AliTrackPoint p;
361   //create initial graph
362   TGraph *grXY=0x0;
363   TGraph *grXZ=0x0;
364   //row -> sector mapping
365   Int_t rowMap[159];
366   for (Int_t irow=0; irow<159; ++irow) rowMap[irow]=-1;
367
368   // 1. Step
369   // Make from the list of global space points graphs in the local frame
370   // one graph per sector is needed
371   for (Int_t ipoint=0; ipoint<npoints; ++ipoint){
372     arrPoints.GetPoint(p,ipoint);
373     Float_t xyz[3] = {p.GetX(),p.GetY(),p.GetZ()};
374     Int_t index[3] = {0,0,0};
375     Int_t row = fTPCParam->GetPadRow(xyz,index);
376     // rotate space point to local frame
377     // the angle is given by the VolumeID which was set in CrateSpacePoints
378     const Int_t sec=p.GetVolumeID();
379     if (row<0 || (sec<36 && row>62) || row>95 ) continue;
380     Double_t angle=((sec%18)*20.+10.)/TMath::RadToDeg();
381     AliTrackPoint pRot=p.Rotate(angle);
382     const Int_t secrow=row+(sec>35)*63;
383     if (rowMap[secrow]==-1) rowMap[secrow]=sec;
384     // check if we need a new graph (sector change)
385     if (secOld!=sec){
386       grXY=new TGraph;
387       grXZ=new TGraph;
388       arrGraphsXY.AddAt(grXY,sec);
389       arrGraphsXZ.AddAt(grXZ,sec);
390     }
391     
392     //add coordinates in local frame for later interpolation
393     grXY->SetPoint(grXY->GetN(), pRot.GetX(), pRot.GetY());
394     grXZ->SetPoint(grXZ->GetN(), pRot.GetX(), pRot.GetZ());
395     secOld=sec;
396   }
397
398   // 2. Step
399   // create in the center of each row a space point by using the graph to interpolate
400   // the the center of the row. This is done in xy and xz
401   TSpline3 *splXY=0x0;
402   TSpline3 *splXZ=0x0;
403   AliTPCclusterMI tempCl;
404   secOld=-1;
405   for (Int_t irow=0; irow<159; ++irow ){
406     const Int_t sec     = rowMap[irow];
407     if (sec==-1) continue;
408     const Int_t secrow = irow<63?irow:irow-63;
409     Double_t localX = fTPCParam->GetPadRowRadii(sec,secrow);
410     // get graph for the current row
411     if (sec!=secOld){
412       delete splXY;
413       splXY=0x0;
414       delete splXZ;
415       splXZ=0x0;
416       
417       grXY=(TGraph*)arrGraphsXY.At(sec);
418       grXZ=(TGraph*)arrGraphsXZ.At(sec);
419       if (!grXY) continue;
420
421 //       if(grXY->GetN()>1 && grXZ->GetN()>1) { //causes segmentation violation if N==1
422 //              splXY=new TSpline3("splXY",grXY);
423 //      splXZ=new TSpline3("splXZ",grXZ);
424 //       }
425 //       else {
426         //TODO: make a cluster also in the sector w only one space point?
427 //      continue;
428         // Double_t tempX=0., tempY = 0., tempZ = 0.;
429         
430         // grXY->GetPoint(0,tempX,localY);
431         // grXZ->GetPoint(0,tempX,localZ);
432 //       }
433
434     }
435     secOld=sec;
436
437     // check we are in an active area
438 //     if (splXY->FindX(localX)<1 || splXZ->FindX(localX)<1) continue;
439     if ( localX<grXY->GetX()[0] || localX>grXY->GetX()[grXY->GetN()-1] || localX<grXZ->GetX()[0] || localX>grXZ->GetX()[grXZ->GetN()-1]) continue;
440     
441     //get interpolated value at the center for the pad row
442     //  using splines
443 //     const Double_t localY=splXY->Eval(localX/*,0x0,"S"*/);
444 //     const Double_t localZ=splXZ->Eval(localX/*,0x0,"S"*/);
445     const Double_t localY=grXY->Eval(localX/*,0x0,"S"*/);
446     const Double_t localZ=grXZ->Eval(localX/*,0x0,"S"*/);
447     Float_t xyz[3]={localX,localY,localZ};
448
449     if (!SetupCluster(tempCl,xyz,sec,t0)) continue;
450     tempCl.SetLabel(tr.GetUniqueID(), 0);
451
452     if (type==0) tr.AddSpacePoint(tempCl);
453     else tr.AddDistortedSpacePoint(tempCl);
454 //     printf("SetupCluster %3d: (%.2f, %.2f, %.2f), %d, %.2f\n",irow,xyz[0],xyz[1],xyz[2],sec,t0);
455   }
456
457   delete splXY;
458   splXY=0x0;
459   delete splXZ;
460   splXZ=0x0;
461   
462 }
463
464 //________________________________________________________________
465 Bool_t AliToyMCEventGenerator::SetupCluster(AliTPCclusterMI &tempCl, Float_t xyz[3], Int_t sec, Double_t t0)
466 {
467   //
468   //
469   //
470
471   // intrinsic cluster resolution is 1mm
472   const Double_t kSigmaY   = 0.1;
473   const Double_t kSigmaZ   = 0.1;
474   const Double_t kMaxZ0    = fTPCParam->GetZLength();
475   //TODO: Get this from the OCDB at some point?
476   const Double_t kDriftVel = fTPCParam->GetDriftV();
477
478   // smear the ideal positions with the cluster resolution
479   xyz[1]+=gRandom->Gaus(0,kSigmaY);
480   xyz[2]+=gRandom->Gaus(0,kSigmaZ);
481   
482   tempCl.SetX(xyz[0]);
483   tempCl.SetY(xyz[1]);
484   tempCl.SetZ(xyz[2]);
485   
486   tempCl.SetSigmaY2(kSigmaY*kSigmaY);
487   tempCl.SetSigmaZ2(kSigmaZ*kSigmaZ);
488
489   // transform from the local coordinates to the coordinates expressed in pad coordinates
490   Int_t index[3] = {0,sec,0};
491   fTPCParam->Transform2to3(xyz,index);
492   fTPCParam->Transform3to4(xyz,index);
493   fTPCParam->Transform4to8(xyz,index);
494
495   const Int_t   row   = index[2];
496   const Int_t   nPads = fTPCParam->GetNPads(sec, row);
497   // pad is fractional, but it needs to be shifted from the center
498   // to the edge of the row
499   const Float_t pad   = xyz[1] + nPads/2;
500
501   tempCl.SetRow(row);
502   tempCl.SetPad(pad);
503   Float_t timeBin=Float_t(t0 + (kMaxZ0-TMath::Abs(tempCl.GetZ()))/kDriftVel);
504   tempCl.SetTimeBin(timeBin); // set time as t0  + drift time from dist z
505   tempCl.SetDetector(sec);
506
507   //check if we are in the active area
508   if (pad<0 || pad>=nPads) return kFALSE;
509
510   return kTRUE;
511 }
512
513 //________________________________________________________________
514 Bool_t AliToyMCEventGenerator::ConnectOutputFile()
515 {
516   //
517   // Create the output file name and tree and connect the event
518   //
519
520   fOutFile = new TFile(fOutputFileName.Data(),"recreate");
521
522   if (!fOutFile || !fOutFile->IsOpen()){
523     delete fOutFile;
524     fOutFile=0x0;
525     return kFALSE;
526   }
527   
528   fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
529   fOutTree->Branch("event","AliToyMCEvent",&fEvent);
530
531   gROOT->cd();
532
533   return kTRUE;
534 }
535
536 //________________________________________________________________
537 Bool_t AliToyMCEventGenerator::CloseOutputFile()
538 {
539   //
540   // close the output file
541   //
542   if (!fOutFile) return kFALSE;
543   fOutFile->Write();
544   fOutFile->Close();
545   delete fOutFile;
546   fOutFile=0x0;
547
548   return kTRUE;
549 }
550
551 //________________________________________________________________
552 void AliToyMCEventGenerator::FillTree()
553 {
554   // fill the tree
555   if (fOutTree&&fEvent) fOutTree->Fill();
556 }
557
558 //________________________________________________________________
559 void AliToyMCEventGenerator::SetSpaceCharge(EEpsilon epsilon, EGasType gasType/*=kNeCO2_9010*/,
560                                             ECollRate collRate/*=k50kHz*/, ECorrection corrType/*=kLookup*/)
561 {
562   //
563   // Set the space charge conditions
564   //
565   fCorrectionFile="$ALICE_ROOT/TPC/Calib/maps/SC";
566   switch (gasType) {
567     case kNeCO2_9010:
568       fCorrectionFile.Append("_NeCO2");
569       break;
570     case kNeCO2N2_90105:
571       fCorrectionFile.Append("_NeCO2N2");
572       break;
573   }
574   switch (epsilon) {
575     case kEps5:
576       fCorrectionFile.Append("_eps5");
577       break;
578     case kEps10:
579       fCorrectionFile.Append("_eps10");
580       break;
581     case kEps20:
582       fCorrectionFile.Append("_eps20");
583       break;
584     case kEps25:
585       fCorrectionFile.Append("_eps25");
586       break;
587     case kEps30:
588       fCorrectionFile.Append("_eps30");
589       break;
590     case kEps35:
591       fCorrectionFile.Append("_eps35");
592       break;
593     case kEps40:
594       fCorrectionFile.Append("_eps40");
595       break;
596   }
597   switch (collRate) {
598     case k50kHz:
599       fCorrectionFile.Append("_50kHz");
600       break;
601   }
602   switch (corrType) {
603     case kLookup:
604       fCorrectionFile.Append("_precal.lookup.root");
605       break;
606     case kSpaceChargeFile:
607       fCorrectionFile.Append("_precal.root");
608       break;
609   }
610 }
611
612 //________________________________________________________________
613 void AliToyMCEventGenerator::InitSpaceCharge()
614 {
615   //
616   // init the space charge conditions
617   // this should be called after the tree was connected
618   //
619
620   if (HasSCList()) {
621     InitSpaceChargeList();
622     return;
623   }
624   
625   AliInfo(Form("Using space charge map file: '%s'",fCorrectionFile.Data()));
626
627   SetCorrectionFromFile(fCorrectionFile, fTPCCorrection);
628
629   if (fOutTree){
630     AliInfo("Attaching space charge map file name to the tree");
631     fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
632   }
633 }
634
635 //________________________________________________________________
636 void AliToyMCEventGenerator::InitSpaceChargeList()
637 {
638   //
639   // init space charge conditions from a list of files
640   // this should be called after the tree was connected
641   //
642
643   std::ifstream file(fSCListFile.Data());
644   TString list;
645   list.ReadFile(file);
646   
647   TObjArray *arr=list.Tokenize("\n");
648   if (!arr->GetEntriesFast()) {
649     delete arr;
650     AliFatal(Form("No SC file initialised. SC list '%s' seems empty",fSCListFile.Data()));
651     return;
652   }
653
654   // it is assumed that in case of an input list
655   // fCorrectionFile contains the name of the average correction
656   // to be then used in the reconstruction
657   if (fOutTree){
658     if (fCorrectionFile.IsNull()) {
659       AliFatal("List of SC files set, but no average map is specified. Use 'SetSpaceChargeFile' to do so");
660       return;
661     }
662     AliInfo("Attaching average space charge map file name to the tree");
663     fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
664   }
665   
666   // check for an average map
667   if (!fCorrectionFile.IsNull()) {
668     SetCorrectionFromFile(fCorrectionFile, fTPCCorrectionAv);
669   }
670   
671   // case of non preread
672   // store the names of the files
673   if (!fPrereadSCList) {
674     if (fSCList) delete fSCList;
675     fSCList=arr;
676     return;
677   }
678
679   // case of preread
680   // load all SC files and set them to the list
681   if (!fSCList) fSCList=new TObjArray;
682   fSCList->SetOwner();
683   fSCList->Delete();
684   
685   for (Int_t ifile=0; ifile<arr->GetEntriesFast(); ++ifile) {
686     TString scname=arr->At(ifile)->GetName();
687     AliTPCCorrection *sc=0x0;
688     SetCorrectionFromFile(scname, sc);
689     if (!sc) continue;
690     fSCList->Add(sc);
691   }
692 }
693
694 //________________________________________________________________
695 void AliToyMCEventGenerator::IterateSC(Int_t ipos)
696 {
697   //
698   // In case a list of SC files is set iterate over them
699   //
700
701   if (!HasSCList()) return;
702
703   if (ipos<0) {
704     if (fOutTree) ipos=fOutTree->GetEntries();
705     else ipos=0;
706   }
707
708   TObject *sc=fSCList->At(ipos%fSCList->GetEntriesFast());
709   AliInfo(Form("Event: %d - SC: %s",ipos,sc->GetName()));
710   // case SC files have been preread
711   if (fPrereadSCList) {
712     fTPCCorrection=(AliTPCCorrection*)sc;
713     return;
714   }
715
716   // case no preread was done
717   TString &file=((TObjString*)sc)->String();
718   delete fTPCCorrection;
719   fTPCCorrection=0x0;
720
721   SetCorrectionFromFile(file, fTPCCorrection);
722   
723   if (!fTPCCorrection) {
724     AliFatal(Form("Could not read SC map from SC file '%s'",file.Data()));
725     return;
726   }
727
728   SetSCScalingFactor();
729 }
730
731 //________________________________________________________________
732 Float_t AliToyMCEventGenerator::GetSCScalingFactor(AliTPCCorrection *corr, AliTPCCorrection *averageCorr, Float_t &chi2)
733 {
734   //
735   // calculate the scaling factor
736   // between the fluctuation map and the average map
737   //
738
739   TLinearFitter fitter(2,"pol1");
740   
741   for (Float_t iz=-245; iz<=245; iz+=10) {
742     Short_t roc=(iz>=0)?0:18;
743     for (Float_t ir=86; ir<250; ir+=10) {
744       for (Float_t iphi=0; iphi<TMath::TwoPi(); iphi+=10*TMath::DegToRad()){
745         Float_t x=ir*(Float_t)TMath::Cos(iphi);
746         Float_t y=ir*(Float_t)TMath::Sin(iphi);
747         Float_t x3[3]    = {x,y,iz};
748         Float_t dx3[3]   = {0.,0.,0.};
749         Float_t dx3av[3] = {0.,0.,0.};
750         
751         corr->GetDistortion(x3,roc,dx3);
752         averageCorr->GetDistortion(x3,roc,dx3av);
753
754         Double_t dr   = TMath::Sqrt(dx3[0]*dx3[0]     +  dx3[1]*dx3[1]);
755         Double_t drav = TMath::Sqrt(dx3av[0]*dx3av[0] +  dx3av[1]*dx3av[1]);
756
757         fitter.AddPoint(&drav,dr);
758       }
759     }
760   }
761   
762   fitter.Eval();
763   chi2 = fitter.GetChisquare()/fitter.GetNpoints();
764   return fitter.GetParameter(1);
765 }
766
767 //________________________________________________________________
768 void AliToyMCEventGenerator::SetSCScalingFactor()
769 {
770   //
771   // if running with a SC list calculate the scaling factor
772   // between the fluctuation map and the average map
773   //
774
775   if ( !(fTPCCorrection && fTPCCorrectionAv && fEvent) ) return;
776
777   // loop over several z, r and phi bins and find a factor that minimises
778   // the distortions between the current map and the average map
779
780   Float_t chi2   = 0.;
781   Float_t factor = GetSCScalingFactor(fTPCCorrection, fTPCCorrectionAv, chi2);
782
783   fEvent->SetSCscale(factor);
784   fEvent->SetSCscaleChi2(chi2);
785 }
786
787 //________________________________________________________________
788 void AliToyMCEventGenerator::SetCorrectionFromFile(TString file, AliTPCCorrection* &corr)
789 {
790   //
791   // read the correction from file and set it to corr
792   //
793
794   corr=0x0;
795   TString corrName("map");
796   
797   // allow for specifying an object name for the AliTPCCorrection in the file name
798   // separated by a ':'
799   TObjArray *arr=file.Tokenize(":");
800   if (arr->GetEntriesFast()>1) {
801     file=arr->At(0)->GetName();
802     corrName=arr->At(1)->GetName();
803   }
804   delete arr;
805   
806   
807   TFile f(file.Data());
808   if (!f.IsOpen()||f.IsZombie()) {
809     printf("E-AliToyMCEventGenerator::Could not open SC file '%s'",file.Data());
810     return;
811   }
812   
813   corr=(AliTPCSpaceCharge3D*)f.Get(corrName.Data());
814   if (corr) corr->SetName(f.GetName());
815 }
816
817