4 #include <TDatabasePDG.h>
7 #include <TGeoGlobalMagField.h>
9 #include <TObjString.h>
12 #include <TObjArray.h>
13 #include <TLinearFitter.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>
28 #include "AliToyMCEvent.h"
29 #include "AliToyMCTrack.h"
31 #include "AliToyMCEventGenerator.h"
33 ClassImp(AliToyMCEventGenerator);
36 AliToyMCEventGenerator::AliToyMCEventGenerator()
42 ,fTPCCorrectionAv(0x0)
45 ,fCorrectionFile("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.lookup.root")
46 ,fOutputFileName("toyMC.root")
49 ,fUseStepCorrection(kFALSE)
50 ,fUseMaterialBudget(kFALSE)
52 ,fPrereadSCList(kFALSE)
53 ,fCalculateScaling(kTRUE)
55 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
56 fTPCParam->ReadGeoMatrices();
59 //________________________________________________________________
60 AliToyMCEventGenerator::AliToyMCEventGenerator(const AliToyMCEventGenerator &gen)
62 ,fTPCParam(gen.fTPCParam)
65 ,fTPCCorrection(gen.fTPCCorrection)
66 ,fTPCCorrectionAv(0x0)
68 ,fSCListFile(gen.fSCListFile)
69 ,fCorrectionFile(gen.fCorrectionFile)
70 ,fOutputFileName(gen.fOutputFileName)
73 ,fUseStepCorrection(gen.fUseStepCorrection)
74 ,fUseMaterialBudget(gen.fUseMaterialBudget)
75 ,fIsLaser(gen.fIsLaser)
76 ,fPrereadSCList(gen.fPrereadSCList)
77 ,fCalculateScaling(gen.fCalculateScaling)
82 //________________________________________________________________
83 AliToyMCEventGenerator::~AliToyMCEventGenerator()
85 if (HasSCList() &&!fPrereadSCList) delete fTPCCorrection;
89 //________________________________________________________________
90 Bool_t AliToyMCEventGenerator::DistortTrack(AliToyMCTrack &trackIn, Double_t t0)
97 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
98 fTPCParam->ReadGeoMatrices();
101 MakeITSClusters(trackIn/*,t0*/);
102 MakeTPCClusters(trackIn, t0);
103 MakeTRDClusters(trackIn/*,t0*/);
107 //________________________________________________________________
108 void AliToyMCEventGenerator::MakeITSClusters(AliToyMCTrack &trackIn/*, Double_t t0*/)
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};
115 //resolution of the point is 10um
116 const Float_t sigmaY = 0.001;
117 const Float_t sigmaZ = 0.001;
121 const Double_t kMaxSnp = 0.85;
122 // const Double_t kMaxZ0 = fTPCParam->GetZLength();
123 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
125 AliExternalTrackParam track(trackIn);
126 Double_t xyz[3] = {0.,0.,0.};
128 if (!AliTrackerBase::PropagateTrackTo(&track,ITSRadii[0],kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
129 AliError(Form("Propagation to %.2f failed\n",ITSRadii[0]));
133 for(Int_t iLayer = 0; iLayer<nITSLayers; iLayer++){
135 if (!AliTrackerBase::PropagateTrackTo(&track,ITSRadii[iLayer],kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
136 AliError(Form("Propagation to %.2f failed\n",ITSRadii[iLayer]));
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
145 // if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
146 if (TMath::Abs(track.GetZ())>lengthITS[iLayer]/2) continue;
149 // smear the ideal positions with the cluster resolution
150 xyz[1]+=gRandom->Gaus(0,sigmaY);
151 xyz[2]+=gRandom->Gaus(0,sigmaZ);
153 Float_t xyzf[3]={static_cast<Float_t>(xyz[0]),static_cast<Float_t>(xyz[1]),static_cast<Float_t>(xyz[2])};
155 SetPoint(xyzf,sigmaY,sigmaZ,point);
157 trackIn.AddITSPoint(point)->SetUniqueID(trackIn.GetUniqueID());
161 //________________________________________________________________
162 void AliToyMCEventGenerator::MakeTRDClusters(AliToyMCTrack &trackIn/*, Double_t t0*/)
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};
171 const Float_t sigmaY = 0.06;
172 const Float_t sigmaZ = 0.2;
176 const Double_t kMaxSnp = 0.85;
177 const Double_t kMaxZ0 = fTPCParam->GetZLength();
178 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
180 AliExternalTrackParam track(trackIn);
181 Double_t xyz[3] = {0.,0.,0.};
183 if (!AliTrackerBase::PropagateTrackTo(&track,TRDRadii[0],kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
184 AliError(Form("Propagation to %.2f failed\n",TRDRadii[0]));
188 for(Int_t iLayer = 0; iLayer<nTRDLayers; iLayer++){
190 if (!AliTrackerBase::PropagateTrackTo(&track,TRDRadii[iLayer],kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterialBudget)) {
191 AliError(Form("Propagation to %.2f failed\n",TRDRadii[iLayer]));
196 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
197 if (TMath::Abs(track.GetZ())>lengthTRD[iLayer]/2) continue;
200 // smear the ideal positions with the cluster resolution
201 xyz[1]+=gRandom->Gaus(0,sigmaY);
202 xyz[2]+=gRandom->Gaus(0,sigmaZ);
204 Float_t xyzf[3]={static_cast<Float_t>(xyz[0]),static_cast<Float_t>(xyz[1]),static_cast<Float_t>(xyz[2])};
206 SetPoint(xyzf,sigmaY,sigmaZ,point);
208 trackIn.AddTRDPoint(point)->SetUniqueID(trackIn.GetUniqueID());
212 //________________________________________________________________
213 void AliToyMCEventGenerator::MakeTPCClusters(AliToyMCTrack &trackIn, Double_t t0)
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
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);
231 //________________________________________________________________
232 void AliToyMCEventGenerator::CreateSpacePoints(AliToyMCTrack &trackIn,
233 AliTrackPointArray &arrUdist,
234 AliTrackPointArray &arrDist)
237 // sample the track from the inner to the outer wall of the TPC
238 // a graph is filled in local coordinates for later
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();
246 const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
247 const Double_t oFCRadius = 254.5;
249 const Float_t kSigmaY=0.1;
250 const Float_t kSigmaZ=0.1;
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.};
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));
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()));
274 //!!! Why is this smeared
276 xyzf[0]=Float_t(xyz[0]);
277 xyzf[1]=Float_t(xyz[1]);
278 xyzf[2]=Float_t(xyz[2]);
280 if (TMath::Abs(track.GetZ())>kMaxZ0) continue;
281 // if (TMath::Abs(track.GetX())<iFCRadius) continue;
282 // if (TMath::Abs(track.GetX())>oFCRadius) continue;
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();
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);
298 fTPCCorrection->GetCorrectionIntegralDz(distPoint,sector,dxyz,5);
299 distPoint[0]-=dxyz[0];
300 distPoint[1]-=dxyz[1];
301 distPoint[2]-=dxyz[2];
303 SetPoint(distPoint,kSigmaY,kSigmaZ,pDist);
304 arrDist.AddPoint(npoints, &pDist);
309 arrUdist.SetUniqueID(npoints);
310 arrDist.SetUniqueID(npoints);
313 //________________________________________________________________
314 void AliToyMCEventGenerator::SetPoint(Float_t xyz[3], Float_t sigmaY, Float_t sigmaZ, AliTrackPoint &point)
317 // make AliTrackPoint out of AliTPCclusterMI
320 //covariance at the local frame
321 //assume 1mm distortion in y and z
323 Float_t cov[6]={0,0,0, sigmaY*sigmaY,0,sigmaZ*sigmaZ};
325 const Float_t alpha = -TMath::ATan2(xyz[1],xyz[0]);
326 const Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
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;
336 // voluem ID to add later ....
338 point.SetCov(newcov);
339 // abuse volume ID for the sector number
340 point.SetVolumeID(fTPCParam->Transform0to1(xyz,i));
342 // TODO: Add sampled dE/dx (use SetCharge)
345 //________________________________________________________________
346 void AliToyMCEventGenerator::ConvertTrackPointsToLocalClusters(AliTrackPointArray &arrPoints,
347 AliToyMCTrack &tr, Double_t t0, Int_t type)
353 const Int_t npoints=Int_t(arrPoints.GetUniqueID());
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();
363 //create initial graph
366 //row -> sector mapping
368 for (Int_t irow=0; irow<159; ++irow) rowMap[irow]=-1;
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)
390 arrGraphsXY.AddAt(grXY,sec);
391 arrGraphsXZ.AddAt(grXZ,sec);
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());
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
405 AliTPCclusterMI tempCl;
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
419 grXY=(TGraph*)arrGraphsXY.At(sec);
420 grXZ=(TGraph*)arrGraphsXZ.At(sec);
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);
428 //TODO: make a cluster also in the sector w only one space point?
430 // Double_t tempX=0., tempY = 0., tempZ = 0.;
432 // grXY->GetPoint(0,tempX,localY);
433 // grXZ->GetPoint(0,tempX,localZ);
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;
443 //get interpolated value at the center for the pad row
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)};
451 if (!SetupCluster(tempCl,xyz,sec,t0)) continue;
452 tempCl.SetLabel(tr.GetUniqueID(), 0);
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);
466 //________________________________________________________________
467 Bool_t AliToyMCEventGenerator::SetupCluster(AliTPCclusterMI &tempCl, Float_t xyz[3], Int_t sec, Double_t t0)
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();
480 // smear the ideal positions with the cluster resolution
481 xyz[1]+=gRandom->Gaus(0,kSigmaY);
482 xyz[2]+=gRandom->Gaus(0,kSigmaZ);
488 tempCl.SetSigmaY2(kSigmaY*kSigmaY);
489 tempCl.SetSigmaZ2(kSigmaZ*kSigmaZ);
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);
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;
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);
509 //check if we are in the active area
510 if (pad<0 || pad>=nPads) return kFALSE;
515 //________________________________________________________________
516 Bool_t AliToyMCEventGenerator::ConnectOutputFile()
519 // Create the output file name and tree and connect the event
522 fOutFile = new TFile(fOutputFileName.Data(),"recreate");
524 if (!fOutFile || !fOutFile->IsOpen()){
530 fOutTree = new TTree("toyMCtree","Tree with toyMC simulation");
531 fOutTree->Branch("event","AliToyMCEvent",&fEvent);
538 //________________________________________________________________
539 Bool_t AliToyMCEventGenerator::CloseOutputFile()
542 // close the output file
544 if (!fOutFile) return kFALSE;
553 //________________________________________________________________
554 void AliToyMCEventGenerator::FillTree()
557 if (fOutTree&&fEvent) fOutTree->Fill();
560 //________________________________________________________________
561 void AliToyMCEventGenerator::SetSpaceCharge(EEpsilon epsilon, EGasType gasType/*=kNeCO2_9010*/,
562 ECollRate collRate/*=k50kHz*/, ECorrection corrType/*=kLookup*/)
565 // Set the space charge conditions
567 fCorrectionFile="$ALICE_ROOT/TPC/Calib/maps/SC";
570 fCorrectionFile.Append("_NeCO2");
573 fCorrectionFile.Append("_NeCO2N2");
578 fCorrectionFile.Append("_eps5");
581 fCorrectionFile.Append("_eps10");
584 fCorrectionFile.Append("_eps20");
587 fCorrectionFile.Append("_eps25");
590 fCorrectionFile.Append("_eps30");
593 fCorrectionFile.Append("_eps35");
596 fCorrectionFile.Append("_eps40");
601 fCorrectionFile.Append("_50kHz");
606 fCorrectionFile.Append("_precal.lookup.root");
608 case kSpaceChargeFile:
609 fCorrectionFile.Append("_precal.root");
614 //________________________________________________________________
615 void AliToyMCEventGenerator::InitSpaceCharge()
618 // init the space charge conditions
619 // this should be called after the tree was connected
623 InitSpaceChargeList();
627 AliInfo(Form("Using space charge map file: '%s'",fCorrectionFile.Data()));
629 SetCorrectionFromFile(fCorrectionFile, fTPCCorrection);
632 AliInfo("Attaching space charge map file name to the tree");
633 fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
637 //________________________________________________________________
638 void AliToyMCEventGenerator::InitSpaceChargeList()
641 // init space charge conditions from a list of files
642 // this should be called after the tree was connected
645 std::ifstream file(fSCListFile.Data());
649 TObjArray *arr=list.Tokenize("\n");
650 if (!arr->GetEntriesFast()) {
652 AliFatal(Form("No SC file initialised. SC list '%s' seems empty",fSCListFile.Data()));
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
660 if (fCorrectionFile.IsNull()) {
661 AliFatal("List of SC files set, but no average map is specified. Use 'SetSpaceChargeFile' to do so");
664 AliInfo("Attaching average space charge map file name to the tree");
665 fOutTree->GetUserInfo()->Add(new TObjString(fCorrectionFile.Data()));
668 // check for an average map
669 if (!fCorrectionFile.IsNull()) {
670 SetCorrectionFromFile(fCorrectionFile, fTPCCorrectionAv);
673 // case of non preread
674 // store the names of the files
675 if (!fPrereadSCList) {
676 if (fSCList) delete fSCList;
682 // load all SC files and set them to the list
683 if (!fSCList) fSCList=new TObjArray;
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);
696 //________________________________________________________________
697 void AliToyMCEventGenerator::IterateSC(Int_t ipos)
700 // In case a list of SC files is set iterate over them
703 if (!HasSCList()) return;
706 if (fOutTree) ipos=fOutTree->GetEntries();
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;
718 // case no preread was done
719 TString &file=((TObjString*)sc)->String();
720 delete fTPCCorrection;
723 SetCorrectionFromFile(file, fTPCCorrection);
725 if (!fTPCCorrection) {
726 AliFatal(Form("Could not read SC map from SC file '%s'",file.Data()));
732 //________________________________________________________________
733 Float_t AliToyMCEventGenerator::GetSCScalingFactor(AliTPCCorrection *corr, AliTPCCorrection *averageCorr, Float_t &chi2)
736 // calculate the scaling factor
737 // between the fluctuation map and the average map
740 TLinearFitter fitter(2,"pol1");
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.};
752 corr->GetDistortion(x3,roc,dx3);
753 averageCorr->GetDistortion(x3,roc,dx3av);
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]);
758 fitter.AddPoint(&drav,dr);
764 chi2 = fitter.GetChisquare()/fitter.GetNpoints();
765 return fitter.GetParameter(1);
768 //________________________________________________________________
769 void AliToyMCEventGenerator::SetSCScalingFactor()
772 // if running with a SC list calculate the scaling factor
773 // between the fluctuation map and the average map
776 if ( !(fCalculateScaling && HasSCList() && fTPCCorrection && fTPCCorrectionAv && fEvent) ) return;
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
782 Float_t factor = GetSCScalingFactor(fTPCCorrection, fTPCCorrectionAv, chi2);
784 fEvent->SetSCscale(factor);
785 fEvent->SetSCscaleChi2(chi2);
788 //________________________________________________________________
789 void AliToyMCEventGenerator::SetCorrectionFromFile(TString file, AliTPCCorrection* &corr)
792 // read the correction from file and set it to corr
796 TString corrName("map");
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();
808 TFile f(file.Data());
809 if (!f.IsOpen()||f.IsZombie()) {
810 printf("E-AliToyMCEventGenerator::Could not open SC file '%s'",file.Data());
814 corr=(AliTPCSpaceCharge3D*)f.Get(corrName.Data());
815 if (corr) corr->SetName(f.GetName());