TPC/Upgrade/macros/fastToyMCMI.C

   1 /*
   2   Fast toy MC for differnt purposes.
   3   testExtrapolation() - to show track extrapolation errror
   4   testdEdxGEM();  -to show the dEdx respolution detoraition as function of the electron trasnparency
   5
   6   .x $HOME/NimStyle.C
   7   .L $ALICE_ROOT/TPC/Upgrade/macros/fastToyMCMI.C+
   8   testExtrapolation();
   9
  10 */
  11 #include "TMath.h"
  12 #include "TVectorD.h"
  13 #include "TLinearFitter.h"
  14 #include "TGeoGlobalMagField.h"
  15 #include "TRandom.h"
  16 #include "TGraphErrors.h"
  17 #include "TStatToolkit.h"
  18
  19 #include "TTreeStream.h"
  20 #include "AliExternalTrackParam.h"
  21 #include "AliCDBManager.h"
  22 #include "AliMagF.h"
  23 #include "AliTrackerBase.h"
  24 #include "TAxis.h"
  25 #include "TLegend.h"
  26 #include "AliGeomManager.h"
  27 #include "TCut.h"
  28 #include "TCanvas.h"
  29 #include "TH2.h"
  30 #include "TF1.h"
  31 #include "TStyle.h"
  32
  33 void testdEdxGEM(const Int_t ntracks=10000, Double_t clNoise=2, Double_t corrNoiseAdd=0.15, Double_t sigmaLandau=0.26){
  34   //
  35   // test dEdx performance as function of the electron transparency.
  36   //
  37   // For simulation standard gRandom->Landau() generator with mean 15 is used
  38   // Charge between pad-rows are corelated via diffusion - qm = (0.15*q-1+0.7*q0+0.15*q1)
  39   // Electrons are randomly absorbed dependending on the transparency parameter
  40   //
  41   // Parameters:
  42   //   clNoise       - noise integrated by cluster (snoise~1 x 4 bins ~ 2)
  43   //   corrNoiseAdd  - additional correlated noise to be added
  44   //   sigmaLandau   - relative sigma of Landau distirbution
  45   // Setup emulation
  46   //   pp   setup testdEdxGEM(20000,2,0.15,0.26)
  47   //   PbPb setup testdEdxGEM(20000,4,0.5,0.26)
  48   //   PbPb setup testdEdxGEM(20000,4,0.6,0.26)
  49   //   PbPb setup testdEdxGEM(20000,4,0.7,0.26)
  50   //
  51
  52   TTreeSRedirector *pcstream = new TTreeSRedirector("testdEdxResolution.root","update");
  53   TH2F * phisdEdx = (TH2F*)pcstream->GetFile()->Get("hisdEdx");
  54   if (!phisdEdx){
  55
  56     TVectorD qVector(160);
  57     TVectorD qVectorCorr(160);
  58     TVectorD qVectorAcc(160);
  59     TVectorD qVectorSorted(160);
  60     Int_t indexes[160];
  61     Int_t ntrans=20;
  62     TVectorD vecdEdx(ntrans);
  63     TVectorD vecT(ntrans);
  64     //
  65
  66     for (Int_t itrack=0; itrack<ntracks; itrack++){
  67       Double_t meanEl=15*(1.+1*gRandom->Rndm());
  68       Double_t sigmaEl=sigmaLandau*meanEl*TMath::Power(15./meanEl,0.25);
  69       for (Int_t irow=0; irow<160; irow++){
  70         qVector[irow]=gRandom->Landau(meanEl, sigmaEl);
  71       }
  72       qVectorCorr[0]=qVector[0];
  73       qVectorCorr[158]=qVector[158];
  74       for (Int_t irow=1; irow<159; irow++){  //corralte measurement via diffusion
  75         qVectorCorr[irow]= 0.15*(qVector[irow-1]+ qVector[irow+1])+0.7*qVector[irow];
  76       }
  77
  78       for (Int_t itrans=0; itrans<ntrans; itrans++){
  79         Double_t transparency=(itrans+1.)/ntrans;
  80         vecT[itrans]=transparency;
  81         for (Int_t irow=0; irow<160; irow++) {
  82           qVectorAcc[irow]=0;
  83           for (Int_t iel=0; iel<TMath::Nint(qVectorCorr[irow]); iel++) {
  84             if (gRandom->Rndm()<transparency)   qVectorAcc[irow]+=1./transparency;
  85           }
  86         }
  87         for (Int_t irow=0; irow<160; irow++) {
  88           qVectorAcc[irow]+=gRandom->Gaus(0, clNoise);
  89           if (qVectorAcc[irow]<0) qVectorAcc[irow]=0;
  90         }
  91         TMath::Sort(160,qVectorAcc.GetMatrixArray(), indexes, kFALSE);
  92         for (Int_t irow=0; irow<160; irow++) {
  93           qVectorSorted[irow]=qVectorAcc[indexes[irow]];
  94         }
  95         vecdEdx[itrans]=TMath::Mean(0.6*160.,   qVectorSorted.GetMatrixArray());
  96         vecdEdx[itrans]+=gRandom->Gaus(0, corrNoiseAdd);
  97       }
  98       (*pcstream)<<"dEdx"<<
  99         "itrack="<<itrack<<              // itrack
 100         "meanEl="<<meanEl<<              //
 101         "sigmaEl="<<sigmaEl<<            //
 102         "vecT.="<<&vecT<<                //
 103         "vecdEdx.="<<&vecdEdx<<
 104         "qVector.="<<&qVector<<
 105         "\n";
 106     }
 107     TTree * tree = (TTree*)(pcstream->GetFile()->Get("dEdx"));
 108     tree->Draw("vecdEdx.fElements/(meanEl):vecT.fElements>>hisdEdx(16,0.199,1,100,0.70,1.50)","meanEl<100","colz");
 109     phisdEdx= (TH2F*)tree->GetHistogram()->Clone();
 110     gStyle->SetOptFit(1);
 111     gStyle->SetOptStat(0);
 112     phisdEdx->Write("hisdEdx");
 113   }
 114   delete pcstream;
 115   gStyle->SetOptStat(0);
 116   gStyle->SetOptFit(1);
 117   pcstream = new TTreeSRedirector("testdEdxResolution.root","update");
 118   phisdEdx = (TH2F*)pcstream->GetFile()->Get("hisdEdx");
 119   TObjArray *arrFit = new TObjArray(3);
 120   phisdEdx->FitSlicesY(0,0,-1,0,"QNR",arrFit);
 121   TH1D * hisRes = (TH1D*)arrFit->At(2);
 122   hisRes->Divide((TH1D*)arrFit->At(1));
 123   hisRes->Scale(100);
 124   hisRes->SetTitle("dEdx resolution");
 125   hisRes->SetMarkerStyle(21);
 126   hisRes->GetXaxis()->SetTitle("Eff. electron  transparency");
 127   hisRes->GetYaxis()->SetTitle("#sigma_{dEdx}/dEdx (%)");
 128   hisRes->SetMinimum(4);
 129   hisRes->SetMaximum(8);
 130   //
 131   TF1 * ftrans = new TF1("ftrans","[0]*x**(-(abs([1])+0.000001))");
 132   ftrans->SetParName(0,"#sigma_{T1}");
 133   ftrans->SetParName(1,"#sigma slope");
 134   ftrans->SetParameters(0.05,1,0.05);
 135   hisRes->SetMarkerStyle(21);
 136   hisRes->SetMarkerSize(0.9);
 137   TCanvas * canvasTrans = new TCanvas("canvasTrans", "canvasTrans",600,500);
 138   canvasTrans->SetTicks(1,1);
 139   hisRes->Fit(ftrans,"","",0.2,1);
 140   canvasTrans->SaveAs("canvasElTrans.pdf");
 141   //TH
 142
 143
 144 }
 145
 146 void testExtrapolation(const Int_t ntracks=10){
 147   //
 148   // check the extrapolation error
 149   // 2 scenarios
 150   //     - ITS extrapolation
 151   //     - ITS + TRD interpolation
 152   //
 153   //
 154   const char *  ocdbpath = "local://$ALICE_ROOT/OCDB/";
 155   AliCDBManager * man = AliCDBManager::Instance();
 156   man->SetDefaultStorage(ocdbpath);
 157   man->SetRun(0);
 158   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,       AliMagF::kBeamTypepp, 2.76/2.));
 159   Double_t   bz=AliTrackerBase::GetBz();
 160   AliGeomManager::LoadGeometry();
 161   //
 162   Double_t rpoints[13]={2.2,   2.8,   3.6,   20,    22,    41,     43,       300,315,330,345,360,375};
 163   Double_t spoints[13]={0.0004,0.0004,0.0004,0.0004,0.0004,0.0004, 0.0004,   0.02,0.02,0.02,0.02,0.02,0.02}; // ITS layers R poition (http://arxiv.org/pdf/1304.1306v3.pdf - pixel scenario)
 164   //
 165   //
 166   //
 167   const Int_t nbins=40;
 168   TFile * f = TFile::Open("testExtrapolationErr.root","update");
 169   if (f->Get("extrapol")==0){
 170     delete f;
 171     f=0;
 172     Double_t cv[21]={0};
 173     Double_t covarSeed[15]={0};
 174     for (Int_t i=0; i<15; i++) covarSeed[i]=0;
 175     covarSeed[0]=0.1;
 176     covarSeed[2]=0.1;
 177     covarSeed[5]=0.001;
 178     covarSeed[9]=0.001;
 179     covarSeed[14]=0.02;
 180     //
 181     Double_t vertex[3]={0,0,0};
 182     TTreeSRedirector *pcstream = new TTreeSRedirector("testExtrapolationErr.root","update");
 183     //
 184     TVectorD vecR(nbins);
 185     TVectorD vecITSErr0(nbins);
 186     TVectorD vecITSTRDErr0(nbins);
 187
 188     for (Int_t itrack=0; itrack<ntracks; itrack++){
 189       //
 190       //Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
 191       Double_t phi = 0;// gRandom->Uniform(0.0, 2*TMath::Pi());
 192       Double_t eta = gRandom->Uniform(-1, 1);
 193       Double_t pt = 2./(gRandom->Rndm()+0.00001);
 194       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
 195       Double_t pxyz[3];
 196       pxyz[0]=pt*TMath::Cos(phi);
 197       pxyz[1]=pt*TMath::Sin(phi);
 198       pxyz[2]=pt*TMath::Tan(theta);
 199       Short_t psign=(gRandom->Rndm()>0.5)? -1.:1.;
 200       AliExternalTrackParam *track= new AliExternalTrackParam(vertex, pxyz, cv, psign);
 201       //
 202       //
 203       //
 204       for (Int_t irbin=0; irbin<nbins; irbin++){
 205         Double_t x0 =85+Double_t(irbin)*(245.-85.)/Double_t(nbins);
 206         Double_t x;
 207         //
 208         // parabolic fit
 209         //
 210         TLinearFitter fitterITS(3,"pol2");
 211         TLinearFitter fitterITSTRD(3,"pol2");
 212         vecR[irbin]=x0;
 213         for (Int_t i=0; i<13; i++) {
 214           Double_t y = 0;
 215           if (track->GetYAt(rpoints[irbin],bz,y)){
 216             x=rpoints[i]-x0;
 217             if (i<7) fitterITS.AddPoint(&x,y,spoints[i]);
 218             fitterITSTRD.AddPoint(&x,y,spoints[i]);
 219           }
 220         }
 221         fitterITS.Eval();
 222         fitterITSTRD.Eval();
 223         vecITSErr0[irbin]=fitterITS.GetParError(0);
 224         vecITSTRDErr0[irbin]=fitterITSTRD.GetParError(0);
 225       }
 226       //
 227       //  estimate q/pt resolution for the ITS+TPC, ITS+TPC+TRD and ITS+TRD scenario
 228       //
 229       AliExternalTrackParam * param = new AliExternalTrackParam(*track);
 230       Double_t *covar = (Double_t*)  param->GetCovariance();
 231       for (Int_t i=0; i<15; i++) covar[i]=covarSeed[i];
 232       AliTrackerBase::PropagateTrackToBxByBz(param, 370,0.13,1,kFALSE);
 233
 234       AliExternalTrackParam *trackITSTPC= new AliExternalTrackParam(*param);
 235       AliExternalTrackParam *trackITSTPCTRD= new AliExternalTrackParam(*param);
 236       AliExternalTrackParam *trackITSTRD= new AliExternalTrackParam(*param);
 237       AliExternalTrackParam *trackTPCTRD= new AliExternalTrackParam(*param);
 238       //
 239       Bool_t tStatus=kTRUE;
 240       for (Int_t idet=2;idet>=0; idet--){
 241         Int_t nlayers=7;
 242         if (idet==1) nlayers=159;
 243         if (idet==2) nlayers=6;
 244         for (Int_t ilayer=nlayers; ilayer>=0; ilayer--){
 245           Double_t rlayer=245.-ilayer;
 246           Double_t slayer=0.1;
 247           if (idet==0) {rlayer=rpoints[ilayer];  slayer=spoints[ilayer];}
 248           if (idet==2) {rlayer=rpoints[ilayer+7]; slayer=spoints[ilayer+7];}
 249           param->PropagateTo(rlayer,bz);
 250           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(param, rlayer,0.13,1,kFALSE);
 251           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTPC, rlayer,0.13,1,kFALSE);
 252           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTPCTRD, rlayer,0.13,1,kFALSE);
 253           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTRD, rlayer,0.13,1,kFALSE);
 254           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackTPCTRD, rlayer,0.13,1,kFALSE);
 255           //if (tStatus==kFALSE) break;
 256           Double_t pointPos[2]={param->GetY(),param->GetZ()};
 257           Double_t pointCov[3]= { slayer*slayer,0,slayer*slayer};
 258           tStatus&=!trackITSTPCTRD->Update(pointPos,pointCov);
 259           if (idet!=2) {
 260             trackITSTPC->Update(pointPos,pointCov);
 261           }
 262           if (idet!=1) {
 263             trackITSTRD->Update(pointPos,pointCov);
 264           }
 265           if (idet!=0) {
 266             trackTPCTRD->Update(pointPos,pointCov);
 267           }
 268         }
 269       }
 270
 271
 272       //vecITSErr0.Print();
 273       (*pcstream)<<"extrapol"<<
 274         "itrack="<<itrack<<
 275         "track.="<<track<<
 276         "vecR.="<<&vecR<<
 277         "vecITSErr0.="<<&vecITSErr0<<
 278         "vecITSTRDErr0.="<<&vecITSTRDErr0<<
 279         "tStatus="<<tStatus<<
 280         "trackITSTPC.="<<trackITSTPC<<
 281         "trackITSTPCTRD.="<<trackITSTPCTRD<<
 282         "trackITSTRD.="<<trackITSTRD<<
 283         "trackTPCTRD.="<<trackTPCTRD<<
 284         "\n";
 285     }
 286     delete pcstream;
 287   }
 288   delete f;
 289   f = TFile::Open("testExtrapolationErr.root","update");
 290   TTree * tree = (TTree*)f->Get("extrapol");
 291   //
 292   TGraphErrors * grITS0    = TStatToolkit::MakeGraphErrors(tree,"10*vecITSErr0.fElements:vecR.fElements","",25,4,0);
 293   TGraphErrors * grITSTRD0 = TStatToolkit::MakeGraphErrors(tree,"10*vecITSTRDErr0.fElements:vecR.fElements","",21,2,0);
 294   //
 295   grITS0->GetXaxis()->SetTitle("radius (cm)");
 296   grITS0->GetYaxis()->SetTitle("#sigma_{r#phi} (mm)");
 297   grITS0->Draw("ap");
 298   grITSTRD0->Draw("p");
 299   TLegend * legend  = new TLegend(0.11,0.6,0.5,0.89,"Track residuals at TPC (q/p_{t}=0)");
 300   legend->AddEntry(grITS0,"ITS extrapolation","p");
 301   legend->AddEntry(grITSTRD0,"ITS-TRD interpolation","p");
 302   legend->Draw();
 303   //
 304   //
 305   //
 306   TCut cut="abs(track.fP[4])<0.25";
 307   TGraphErrors * grITSTPCqPt    = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTPC.fC[14]):abs(track.fP[4])",cut,20,1,0.8);
 308   TGraphErrors * grITSTRDqPt    = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTRD.fC[14]):abs(track.fP[4])",cut,21,4,0.8);
 309   TGraphErrors * grITSTPCTRDqPt = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTPCTRD.fC[14]):abs(track.fP[4])",cut,24,2,0.8);
 310   TGraphErrors * grTPCTRDqPt    = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackTPCTRD.fC[14]):abs(track.fP[4])",cut,25,3,0.8);
 311   grITSTPCqPt->GetXaxis()->SetTitle("q/p_{t} (c/GeV)");
 312   grITSTPCqPt->GetYaxis()->SetTitle("#sigma_{q/p_{t}} (c/GeV)");
 313   grITSTPCqPt->SetMaximum(0.003);
 314   TCanvas * canvasResol = new TCanvas("canvasResol","canvasResol",800,600);
 315   canvasResol->cd();
 316   canvasResol->SetLeftMargin(0.15);
 317   grITSTPCqPt->GetYaxis()->SetTitleOffset(1.3);
 318
 319   {
 320     grITSTPCqPt->Draw("ap");
 321     grITSTRDqPt->Draw("p");
 322     grITSTPCTRDqPt->Draw("p");
 323     grTPCTRDqPt->Draw("p");
 324     TLegend * legendQpt  = new TLegend(0.41,0.11,0.89,0.4,"q/p_{t} resolution (from covariance matrix) ");
 325     legendQpt->AddEntry(grITSTPCqPt,"ITS+TPC","p");
 326     legendQpt->AddEntry(grITSTRDqPt,"ITS+TRD","p");
 327     legendQpt->AddEntry(grITSTPCTRDqPt,"ITS+TPC+TRD","p");
 328     legendQpt->AddEntry(grTPCTRDqPt,"TPC+TRD","p");
 329     legendQpt->Draw();
 330   }
 331
 332 }