TPC/Upgrade/macros/fastToyMCMI.C

   1 /*
   2   Fast toy MC for different purposes. primary goal prepare the motivation figures for the TPC TDR and
   3   expalantion internal note.
   4
   5   testExtrapolationError() - to show track extrapolation errror
   6   testdEdxGEM();  -to show the dEdx respolution detoriation as function of the electron trasnparency
   7
   8
   9   .x $HOME/NimStyle.C
  10   .L $ALICE_ROOT/TPC/Upgrade/macros/fastToyMCMI.C+
  11
  12 */
  13 #include "TStyle.h"
  14 #include "TCut.h"
  15 #include "TCanvas.h"
  16 #include "TMath.h"
  17 #include "TVectorD.h"
  18 #include "TLinearFitter.h"
  19 #include "TGeoGlobalMagField.h"
  20 #include "TRandom.h"
  21 #include "TGraphErrors.h"
  22 #include "TStatToolkit.h"
  23
  24 #include "TTreeStream.h"
  25 #include "AliExternalTrackParam.h"
  26 #include "AliCDBManager.h"
  27 #include "AliMagF.h"
  28 #include "AliTrackerBase.h"
  29 #include "TAxis.h"
  30 #include "TLegend.h"
  31 #include "AliGeomManager.h"
  32 #include "TCut.h"
  33 #include "TCanvas.h"
  34 #include "TH2.h"
  35 #include "TF1.h"
  36 #include "TStyle.h"
  37
  38
  39
  40 void testdEdxGEM(const Int_t ntracks=10000, Double_t clNoise=2, Double_t corrNoiseAdd=0.15, Double_t sigmaLandau=0.26);
  41 void testExtrapolationError(Int_t ntracks, Int_t useGeo=kFALSE, Int_t seed=0);
  42
  43
  44 void fastToyMCMI(Int_t action=1, Float_t arg0=200, Float_t arg1=0, Float_t arg2=0){
  45   //
  46   //
  47   //
  48   if (action==1)  testExtrapolationError(arg0,arg1,arg2);
  49 }
  50
  51
  52 void testdEdxGEM(const Int_t ntracks, Double_t clNoise, Double_t corrNoiseAdd, Double_t sigmaLandau){
  53   //
  54   // test dEdx performance as function of the electron transparency.
  55   //
  56   // For simulation standard gRandom->Landau() generator with mean 15 is used
  57   // Charge between pad-rows are corelated via diffusion - qm = (0.15*q-1+0.7*q0+0.15*q1)
  58   // Electrons are randomly absorbed dependending on the transparency parameter
  59   //
  60   // Parameters:
  61   //   clNoise       - noise integrated by cluster (snoise~1 x 4 bins ~ 2)
  62   //   corrNoiseAdd  - additional correlated noise to be added
  63   //   sigmaLandau   - relative sigma of Landau distirbution
  64   // Setup emulation
  65   //   pp   setup testdEdxGEM(20000,2,0.15,0.26)
  66   //   PbPb setup testdEdxGEM(20000,4,0.5,0.26)
  67   //   PbPb setup testdEdxGEM(20000,4,0.6,0.26)
  68   //   PbPb setup testdEdxGEM(20000,4,0.7,0.26)
  69   //
  70
  71   TTreeSRedirector *pcstream = new TTreeSRedirector("testdEdxResolution.root","update");
  72   TH2F * phisdEdx = (TH2F*)pcstream->GetFile()->Get("hisdEdx");
  73   if (!phisdEdx){
  74
  75     TVectorD qVector(160);
  76     TVectorD qVectorCorr(160);
  77     TVectorD qVectorAcc(160);
  78     TVectorD qVectorSorted(160);
  79     Int_t indexes[160];
  80     Int_t ntrans=20;
  81     TVectorD vecdEdx(ntrans);
  82     TVectorD vecT(ntrans);
  83     //
  84
  85     for (Int_t itrack=0; itrack<ntracks; itrack++){
  86       Double_t meanEl=15*(1.+1*gRandom->Rndm());
  87       Double_t sigmaEl=sigmaLandau*meanEl*TMath::Power(15./meanEl,0.25);
  88       for (Int_t irow=0; irow<160; irow++){
  89         qVector[irow]=gRandom->Landau(meanEl, sigmaEl);
  90       }
  91       qVectorCorr[0]=qVector[0];
  92       qVectorCorr[158]=qVector[158];
  93       for (Int_t irow=1; irow<159; irow++){  //corralte measurement via diffusion
  94         qVectorCorr[irow]= 0.15*(qVector[irow-1]+ qVector[irow+1])+0.7*qVector[irow];
  95       }
  96
  97       for (Int_t itrans=0; itrans<ntrans; itrans++){
  98         Double_t transparency=(itrans+1.)/ntrans;
  99         vecT[itrans]=transparency;
 100         for (Int_t irow=0; irow<160; irow++) {
 101           qVectorAcc[irow]=0;
 102           for (Int_t iel=0; iel<TMath::Nint(qVectorCorr[irow]); iel++) {
 103             if (gRandom->Rndm()<transparency)   qVectorAcc[irow]+=1./transparency;
 104           }
 105         }
 106         for (Int_t irow=0; irow<160; irow++) {
 107           qVectorAcc[irow]+=gRandom->Gaus(0, clNoise);
 108           if (qVectorAcc[irow]<0) qVectorAcc[irow]=0;
 109         }
 110         TMath::Sort(160,qVectorAcc.GetMatrixArray(), indexes, kFALSE);
 111         for (Int_t irow=0; irow<160; irow++) {
 112           qVectorSorted[irow]=qVectorAcc[indexes[irow]];
 113         }
 114         vecdEdx[itrans]=TMath::Mean(0.6*160.,   qVectorSorted.GetMatrixArray());
 115         vecdEdx[itrans]+=gRandom->Gaus(0, corrNoiseAdd);
 116       }
 117       (*pcstream)<<"dEdx"<<
 118         "itrack="<<itrack<<              // itrack
 119         "meanEl="<<meanEl<<              //
 120         "sigmaEl="<<sigmaEl<<            //
 121         "vecT.="<<&vecT<<                //
 122         "vecdEdx.="<<&vecdEdx<<
 123         "qVector.="<<&qVector<<
 124         "\n";
 125     }
 126     TTree * tree = (TTree*)(pcstream->GetFile()->Get("dEdx"));
 127     tree->Draw("vecdEdx.fElements/(meanEl):vecT.fElements>>hisdEdx(16,0.199,1,100,0.70,1.50)","meanEl<100","colz");
 128     phisdEdx= (TH2F*)tree->GetHistogram()->Clone();
 129     gStyle->SetOptFit(1);
 130     gStyle->SetOptStat(0);
 131     phisdEdx->Write("hisdEdx");
 132   }
 133   delete pcstream;
 134   gStyle->SetOptStat(0);
 135   gStyle->SetOptFit(1);
 136   pcstream = new TTreeSRedirector("testdEdxResolution.root","update");
 137   phisdEdx = (TH2F*)pcstream->GetFile()->Get("hisdEdx");
 138   TObjArray *arrFit = new TObjArray(3);
 139   phisdEdx->FitSlicesY(0,0,-1,0,"QNR",arrFit);
 140   TH1D * hisRes = (TH1D*)arrFit->At(2);
 141   hisRes->Divide((TH1D*)arrFit->At(1));
 142   hisRes->Scale(100);
 143   hisRes->SetTitle("dEdx resolution");
 144   hisRes->SetMarkerStyle(21);
 145   hisRes->GetXaxis()->SetTitle("Eff. electron  transparency");
 146   hisRes->GetYaxis()->SetTitle("#sigma_{dEdx}/dEdx (%)");
 147   hisRes->SetMinimum(4);
 148   hisRes->SetMaximum(8);
 149   //
 150   TF1 * ftrans = new TF1("ftrans","[0]*x**(-(abs([1])+0.000001))");
 151   ftrans->SetParName(0,"#sigma_{T1}");
 152   ftrans->SetParName(1,"#sigma slope");
 153   ftrans->SetParameters(0.05,1,0.05);
 154   hisRes->SetMarkerStyle(21);
 155   hisRes->SetMarkerSize(0.9);
 156   TCanvas * canvasTrans = new TCanvas("canvasTrans", "canvasTrans",600,500);
 157   canvasTrans->SetTicks(1,1);
 158   hisRes->Fit(ftrans,"","",0.2,1);
 159   canvasTrans->SaveAs("canvasElTrans.pdf");
 160   //TH
 161
 162
 163 }
 164
 165 void testExtrapolationError(Int_t ntracks, Int_t useGeo, Int_t seed){
 166   //
 167   // check the extrapolation error
 168   // 2 scenarios
 169   //     - ITS extrapolation
 170   //     - ITS + TRD interpolation
 171   //
 172   //
 173   gRandom->SetSeed(seed);
 174   const char *  ocdbpath = "local://$ALICE_ROOT/OCDB/";
 175   AliCDBManager * man = AliCDBManager::Instance();
 176   man->SetDefaultStorage(ocdbpath);
 177   man->SetRun(0);
 178   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,       AliMagF::kBeamTypepp, 2.76/2.));
 179   Double_t   bz=AliTrackerBase::GetBz();
 180   if ( useGeo)   AliGeomManager::LoadGeometry("geometry.root");
 181   if (!useGeo)   AliGeomManager::LoadGeometry();
 182   //
 183   //  Double_t rpoints[13]={2.2,   2.8,   3.6,   20,    22,    41,     43,       300,315,330,345,360,375};
 184   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)
 185   //
 186   Double_t rpoints[13]={2.32,   3.13,   3.91,   19.41,    24.71,    35.33,     40.53,       300,315,330,345,360,375};
 187
 188   //
 189   //
 190   //
 191   const Int_t nbins=40;
 192   TFile * f = TFile::Open("testExtrapolationErr.root","update");
 193   if (f->Get("extrapol")==0){
 194     delete f;
 195     f=0;
 196     Double_t cv[21]={0};
 197     Double_t covarSeed[15]={0};
 198     for (Int_t i=0; i<15; i++) covarSeed[i]=0;
 199     covarSeed[0]=0.1;
 200     covarSeed[2]=0.1;
 201     covarSeed[5]=0.001;
 202     covarSeed[9]=0.001;
 203     covarSeed[14]=0.03;
 204     //
 205     Double_t vertex[3]={0,0,0};
 206     TTreeSRedirector *pcstream = new TTreeSRedirector("testExtrapolationErr.root","update");
 207     //
 208     TVectorD vecR(nbins);
 209     TVectorD vecITSErr0(nbins);
 210     TVectorD vecITSTRDErr0(nbins);
 211     //
 212     for (Int_t itrack=0; itrack<ntracks; itrack++){
 213       //
 214       //Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
 215       if (itrack%20==0) printf("processing\t%d\n", itrack);
 216       Double_t phi = (gRandom->Rndm()-1)*TMath::Pi()/18;
 217       Double_t eta = gRandom->Uniform(-1, 1);
 218       Double_t pt = 1./(gRandom->Rndm()+0.00001);
 219       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
 220       Double_t pxyz[3];
 221       pxyz[0]=pt*TMath::Cos(phi);
 222       pxyz[1]=pt*TMath::Sin(phi);
 223       pxyz[2]=pt*TMath::Tan(theta);
 224       Short_t psign=(gRandom->Rndm()>0.5)? -1.:1.;
 225       AliExternalTrackParam *track= new AliExternalTrackParam(vertex, pxyz, cv, psign);
 226       Double_t alpha=TMath::ATan2(pxyz[1],pxyz[0]);
 227       track->Rotate(alpha);
 228       //
 229       // 0.) Estimate the error using the ITS extrapolation and ITS+TRD interpolation - neglecting the MS -inf. momanta tracks
 230       //
 231       for (Int_t irbin=0; irbin<nbins; irbin++){
 232         Double_t x0 =85+Double_t(irbin)*(245.-85.)/Double_t(nbins);
 233         Double_t x;
 234         //
 235         // parabolic fit
 236         //
 237         TLinearFitter fitterITS(3,"pol2");
 238         TLinearFitter fitterITSTRD(3,"pol2");
 239         vecR[irbin]=x0;
 240         for (Int_t i=0; i<13; i++) {
 241           Double_t y = 0;
 242           if (track->GetYAt(rpoints[irbin],bz,y)){
 243             x=rpoints[i]-x0;
 244             if (i<7) fitterITS.AddPoint(&x,y,spoints[i]);
 245             fitterITSTRD.AddPoint(&x,y,spoints[i]);
 246           }
 247         }
 248         fitterITS.Eval();
 249         fitterITSTRD.Eval();
 250         vecITSErr0[irbin]=fitterITS.GetParError(0);
 251         vecITSTRDErr0[irbin]=fitterITSTRD.GetParError(0);
 252       }
 253       //
 254       //  1.) estimate q/pt resolution for the ITS+TPC, ITS+TPC+TRD and ITS+TRD scenario
 255       //
 256       AliExternalTrackParam * param = new AliExternalTrackParam(*track);
 257       Double_t *covar = (Double_t*)  param->GetCovariance();
 258       for (Int_t i=0; i<15; i++) covar[i]=covarSeed[i];
 259       AliTrackerBase::PropagateTrackToBxByBz(param, 370,0.13,1,kFALSE);
 260
 261       AliExternalTrackParam *trackITSTPC= new AliExternalTrackParam(*param);
 262       AliExternalTrackParam *trackITSTPCTRD= new AliExternalTrackParam(*param);
 263       AliExternalTrackParam *trackITSTRD= new AliExternalTrackParam(*param);
 264       AliExternalTrackParam *trackTPCTRD= new AliExternalTrackParam(*param);
 265       //
 266       Bool_t tStatus=kTRUE;
 267       for (Int_t idet=2;idet>=0; idet--){
 268         Int_t nlayers=7;
 269         if (idet==1) nlayers=159;
 270         if (idet==2) nlayers=6;
 271         for (Int_t ilayer=nlayers; ilayer>=0; ilayer--){
 272           Double_t rlayer=245.-ilayer;
 273           Double_t slayer=0.1;
 274           if (idet==0) {rlayer=rpoints[ilayer];  slayer=spoints[ilayer];}
 275           if (idet==2) {rlayer=rpoints[ilayer+7]; slayer=spoints[ilayer+7];}
 276           param->PropagateTo(rlayer,bz);
 277           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(param, rlayer,0.13,1,kFALSE);
 278           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTPC, rlayer,0.13,1,kFALSE);
 279           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTPCTRD, rlayer,0.13,1,kFALSE);
 280           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTRD, rlayer,0.13,1,kFALSE);
 281           tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackTPCTRD, rlayer,0.13,1,kFALSE);
 282           //if (tStatus==kFALSE) break;
 283           Double_t pointPos[2]={param->GetY(),param->GetZ()};
 284           Double_t pointCov[3]= { slayer*slayer,0,slayer*slayer};
 285           tStatus&=!trackITSTPCTRD->Update(pointPos,pointCov);
 286           if (idet!=2) {
 287             trackITSTPC->Update(pointPos,pointCov);
 288           }
 289           if (idet!=1) {
 290             trackITSTRD->Update(pointPos,pointCov);
 291           }
 292           if (idet!=0) {
 293             trackTPCTRD->Update(pointPos,pointCov);
 294           }
 295         }
 296       }
 297       //
 298       // 2.) Estimate propagation error at given radius
 299       //
 300       //  2.a) Fit ITS and TRD track (gauss smeared point error 0 not MS used)
 301       AliExternalTrackParam *trackITS= new AliExternalTrackParam(*track);
 302       AliExternalTrackParam *trackITS0= new AliExternalTrackParam(*track);
 303       covar = (Double_t*)  trackITS->GetCovariance();
 304       for (Int_t i=0; i<15; i++) covar[i]=covarSeed[i];
 305       AliExternalTrackParam *trackTRD= new AliExternalTrackParam(*param);
 306       AliExternalTrackParam *trackTRD0= new AliExternalTrackParam(*param);
 307       {for (Int_t ilayer=0; ilayer<7; ilayer++){
 308           AliTrackerBase::PropagateTrackToBxByBz(trackITS, rpoints[ilayer],0.13,1,kFALSE);
 309           AliTrackerBase::PropagateTrackToBxByBz(trackITS0, rpoints[ilayer],0.13,1,kFALSE);
 310           Double_t pointPos[2]={trackITS0->GetY()+gRandom->Gaus(0,spoints[ilayer]),trackITS0->GetZ()+gRandom->Gaus(0,spoints[ilayer])};
 311           Double_t pointCov[3]= {spoints[ilayer]*spoints[ilayer],0,spoints[ilayer]*spoints[ilayer]};
 312           trackITS->Update(pointPos,pointCov);
 313         }}
 314       {for (Int_t ilayer=5; ilayer>=0; ilayer--){
 315           AliTrackerBase::PropagateTrackToBxByBz(trackTRD, rpoints[ilayer+7],0.13,1,kFALSE);
 316           AliTrackerBase::PropagateTrackToBxByBz(trackTRD0, rpoints[ilayer+7],0.13,1,kFALSE);
 317           Double_t pointPos[2]={trackTRD0->GetY()+gRandom->Gaus(0,spoints[ilayer+7]),trackTRD0->GetZ()+gRandom->Gaus(0,spoints[ilayer+7])};
 318           Double_t pointCov[3]= {spoints[ilayer+7]*spoints[ilayer+7],0,spoints[ilayer+7]*spoints[ilayer+7]};
 319           trackTRD->Update(pointPos,pointCov);
 320         }}
 321       //
 322       // 2.b) get ITS extrapolation and TRD interpolation errors at random radisues positions
 323       //
 324       const Int_t npoints=20;
 325       TVectorD vecRKalman(npoints);
 326       TVectorD vecITSDeltaY(npoints),    vecITSTRDDeltaY(npoints);
 327       TVectorD vecITSErrY(npoints),      vecITSTRDErrY(npoints);
 328       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 329         Double_t rLayer= 85.+gRandom->Rndm()*(245.-85.);
 330         AliExternalTrackParam *trackITSP= new AliExternalTrackParam(*trackITS);
 331         AliExternalTrackParam *trackITSP0= new AliExternalTrackParam(*trackITS0);
 332         AliExternalTrackParam *trackTRDP= new AliExternalTrackParam(*trackTRD);
 333         AliTrackerBase::PropagateTrackToBxByBz(trackITSP, rLayer,0.13,1,kFALSE);
 334         AliTrackerBase::PropagateTrackToBxByBz(trackITSP0, rLayer,0.13,1,kFALSE);
 335         AliTrackerBase::PropagateTrackToBxByBz(trackTRDP, rLayer,0.13,1,kFALSE);
 336         vecRKalman[ipoint]=rLayer;
 337         trackTRDP->Rotate(trackITSP->GetAlpha());
 338         vecITSDeltaY[ipoint]=trackITSP->GetY()-trackITSP0->GetY();
 339         vecITSErrY[ipoint]=TMath::Sqrt(trackITSP->GetSigmaY2());
 340         AliTrackerBase::UpdateTrack(*trackITSP, *trackTRDP);
 341         vecITSTRDDeltaY[ipoint]=trackITSP->GetY()-trackITSP0->GetY();
 342         vecITSTRDErrY[ipoint]=TMath::Sqrt(trackITSP->GetSigmaY2());
 343       }
 344       //
 345
 346       //vecITSErr0.Print();
 347       (*pcstream)<<"extrapol"<<
 348         "itrack="<<itrack<<
 349         "track.="<<track<<
 350         "vecR.="<<&vecR<<
 351         //
 352         "vecITSErr0.="<<&vecITSErr0<<
 353         "vecITSTRDErr0.="<<&vecITSTRDErr0<<
 354         "tStatus="<<tStatus<<
 355         "trackITSTPC.="<<trackITSTPC<<
 356         "trackITSTPCTRD.="<<trackITSTPCTRD<<
 357         "trackITSTRD.="<<trackITSTRD<<
 358         "trackTPCTRD.="<<trackTPCTRD<<
 359         //
 360         // kalman extapoltation einterplotaltion studies
 361         // extrapolation and ITSTRDitnerpolation at different radiuses
 362         "verRKalman.="<<&vecRKalman<<
 363         "vecITSDeltaY.="<<&vecITSDeltaY<<
 364         "vecITSErrY.="<<&vecITSErrY<<
 365         "vecITSTRDDeltaY.="<<&vecITSTRDDeltaY<<
 366         "vecITSTRDErrY.="<<&vecITSTRDErrY<<
 367         "\n";
 368     }
 369     delete pcstream;
 370   }
 371   delete f;
 372
 373   gStyle->SetOptTitle(0);
 374   f = TFile::Open("testExtrapolationErr.root","update");
 375   TTree * tree = (TTree*)f->Get("extrapol");
 376   tree->SetMarkerStyle(25);
 377   tree->SetMarkerSize(0.3);
 378
 379   //
 380   TGraphErrors * grITS0    = TStatToolkit::MakeGraphErrors(tree,"10*vecITSErr0.fElements:vecR.fElements","",25,4,0);
 381   TGraphErrors * grITSTRD0 = TStatToolkit::MakeGraphErrors(tree,"10*vecITSTRDErr0.fElements:vecR.fElements","",21,2,0);
 382   //
 383   grITS0->GetXaxis()->SetTitle("radius (cm)");
 384   grITS0->GetYaxis()->SetTitle("#sigma_{r#varphi} (mm)");
 385   grITS0->Draw("ap");
 386   grITSTRD0->Draw("p");
 387   TLegend * legend  = new TLegend(0.11,0.65,0.55,0.89,"Track residuals at TPC (q/p_{t}=0)");
 388   legend->AddEntry(grITS0,"ITS extrapolation","p");
 389   legend->AddEntry(grITSTRD0,"ITS-TRD interpolation","p");
 390   legend->SetFillColor(10);
 391   legend->Draw();
 392   //
 393   //
 394   //
 395   TCut cut="abs(track.fP[4])<0.25";
 396   TGraphErrors * grITSTPCqPt    = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTPC.fC[14]):abs(track.fP[4])",cut,20,1,0.8);
 397   TGraphErrors * grITSTRDqPt    = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTRD.fC[14]):abs(track.fP[4])",cut,21,4,0.8);
 398   TGraphErrors * grITSTPCTRDqPt = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTPCTRD.fC[14]):abs(track.fP[4])",cut,24,2,0.8);
 399   TGraphErrors * grTPCTRDqPt    = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackTPCTRD.fC[14]):abs(track.fP[4])",cut,25,3,0.8);
 400   grITSTPCqPt->GetXaxis()->SetTitle("q/p_{t} (c/GeV)");
 401   grITSTPCqPt->GetYaxis()->SetTitle("#sigma_{q/p_{t}} (c/GeV)");
 402   grITSTPCqPt->SetMaximum(0.003);
 403   TCanvas * canvasResol = new TCanvas("canvasResol","canvasResol",800,600);
 404   canvasResol->cd();
 405   canvasResol->SetLeftMargin(0.15);
 406   grITSTPCqPt->GetYaxis()->SetTitleOffset(1.3);
 407
 408   {
 409     grITSTPCqPt->Draw("ap");
 410     grITSTRDqPt->Draw("p");
 411     grITSTPCTRDqPt->Draw("p");
 412     grTPCTRDqPt->Draw("p");
 413     TLegend * legendQpt  = new TLegend(0.41,0.11,0.89,0.4,"q/p_{t} resolution (from covariance matrix) ");
 414     legendQpt->AddEntry(grITSTPCqPt,"ITS+TPC","p");
 415     legendQpt->AddEntry(grITSTRDqPt,"ITS+TRD","p");
 416     legendQpt->AddEntry(grITSTPCTRDqPt,"ITS+TPC+TRD","p");
 417     legendQpt->AddEntry(grTPCTRDqPt,"TPC+TRD","p");
 418     legendQpt->Draw();
 419   }
 420   //
 421   //
 422   TCanvas *canvasResolution = new TCanvas("canvasResolution","canvasResolution",600,600);
 423   canvasResolution->Divide(1,3);
 424   //
 425   canvasResolution->cd(1);
 426   tree->Draw("vecITSErrY.fElements:verRKalman.fElements:abs(track.fP[4])","","colz");
 427   canvasResolution->cd(2);
 428   tree->Draw("vecITSTRDErrY.fElements:verRKalman.fElements:abs(track.fP[4])","","colz");
 429   canvasResolution->cd(3);
 430   tree->Draw("vecITSTRDErrY.fElements/vecITSErrY.fElements:verRKalman.fElements:abs(track.fP[4])","","colz");
 431   canvasResolution->SaveAs("canvasResolution.pdf");
 432   canvasResolution->SaveAs("canvasResolution.png");
 433   //
 434   //
 435   //
 436   TStatToolkit toolkit;
 437   Double_t chi2=0;
 438   Int_t    npoints=0;
 439   TVectorD param;
 440   TMatrixD covar;
 441   //
 442   TCut cut="";
 443   TString  fstringITS="";
 444   fstringITS+="abs(track.fP[4])*(verRKalman.fElements-40)^2++";
 445   fstringITS+="(verRKalman.fElements-40)^2++";
 446   TString * fitErrITS = TStatToolkit::FitPlane(tree,"vecITSErrY.fElements:(0.1+track.fP[4]**2)", fstringITS.Data(),"verRKalman.fElements>0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
 447   tree->SetAlias("fitErrITS",fitErrITS->Data());
 448   fitErrITS->Tokenize("++")->Print();
 449   tree->Draw("vecITSErrY.fElements/fitErrITS:verRKalman.fElements:abs(track.fP[4])","","colz");
 450   //
 451   //
 452   TString  fstringITSTRD="";
 453   fstringITSTRD+="abs(track.fP[4])++";
 454   fstringITSTRD+="abs(track.fP[4]*verRKalman.fElements)++";
 455   fstringITSTRD+="abs(track.fP[4]*verRKalman.fElements**2)++";
 456   for (Int_t iter=0; iter<3; iter++){
 457     TCut cutOut="1";
 458     if (iter>0) cutOut="abs(vecITSTRDErrY.fElements/fitErrITSTRD-1)<0.8";
 459     TString * fitErrITSTRD = TStatToolkit::FitPlane(tree,"vecITSTRDErrY.fElements:(0.1+track.fP[4]**2)", fstringITSTRD.Data(),"verRKalman.fElements>0"+cutOut, chi2,npoints,param,covar,-1,0,180000 , kFALSE);
 460     tree->SetAlias("fitErrITSTRD",fitErrITSTRD->Data());
 461     fitErrITSTRD->Tokenize("++")->Print();
 462   }
 463   tree->Draw("vecITSTRDErrY.fElements/fitErrITSTRD:verRKalman.fElements:abs(track.fP[4])","","colz");
 464
 465
 466   TCanvas *canvasResolutionFit = new TCanvas("canvasResolutionFit","canvasResolutionFit",700,700);
 467   canvasResolutionFit->Divide(1,3);
 468   canvasResolutionFit->cd(1);
 469   tree->Draw("fitErrITS:verRKalman.fElements:abs(track.fP[4])>>hITS","","colz");
 470   htemp = (TH2F*)gPad->GetPrimitive("hITS")->Clone();
 471   htemp->GetXaxis()->SetTitle("#it{r} (cm)");
 472   htemp->GetYaxis()->SetTitle("#sigma_{r#phi} (cm) ITS");
 473   htemp->GetZaxis()->SetTitle("q/#it{p_{t}} (#it{c}/GeV)");
 474   htemp->Draw("colz");
 475   canvasResolutionFit->cd(2);
 476   tree->Draw("fitErrITSTRD:verRKalman.fElements:abs(track.fP[4])>>hITSTRD","","colz");
 477   htemp = (TH2F*)gPad->GetPrimitive("hITSTRD")->Clone();
 478   htemp->GetXaxis()->SetTitle("#it{r} (cm)");
 479   htemp->GetYaxis()->SetTitle("#sigma_{r#phi} (cm) ITS+TRD");
 480   htemp->GetZaxis()->SetTitle("q/#it{p_{t}} (#it{c}/GeV)");
 481   htemp->Draw("colz");
 482   canvasResolutionFit->cd(3);
 483   tree->Draw("fitErrITSTRD:verRKalman.fElements:abs(track.fP[4])>>hITSTRD","","colz");
 484   htemp = (TH2F*)gPad->GetPrimitive("hITSTRD")->Clone();
 485   htemp->GetXaxis()->SetTitle("#it{r} (cm)");
 486   htemp->GetYaxis()->SetTitle("#sigma_{r#phi} (cm) ITS+TRD");
 487   htemp->GetZaxis()->SetTitle("q/#it{p_{t}} (#it{c}/GeV)");
 488   htemp->Draw("colz");
 489   canvasResolutionFit->SaveAs("canvasResolutionFit.pdf");
 490   canvasResolutionFit->SaveAs("canvasResolutionFit.png");
 491
 492 }