4 AliRunLoader *loader = AliRunLoader::Open("galice.root");
6 gAlice = loader->GetAliRun();
7 AliMagF * field = gAlice->Field();
8 AliTracker::SetFieldMap(gAlice->Field(),1);
9 TFile fg("geometry.root");
10 gGeoManager = (TGeoManager*)fg.Get("Geo");
11 if (!gGeoManager) TGeoManager::Import("geometry.root");
13 .L AliESDComparisonMI.C+
14 .L AliTRDComparison.C+g
15 //.L AliESDComparisonPIC.C+
25 TFile fff("TRDdebug.root");
26 TTree * tree =(TTree*)fff.Get("tracklet");
27 TTree * tree2 =(TTree*)fff.Get("Tracks");
28 TTree * treefind =(TTree*)fff.Get("Find");
29 TTree * treebug =(TTree*)fff.Get("Bug");
30 TTree * treetofs =(TTree*)fff.Get("tofseed");
31 AliComparisonDraw compfind;
32 compfind.fTree = treefind;
33 AliComparisonDraw comptr;
36 TFile fff2("TOFdebug.root");
37 TTree * treetof =(TTree*)fff2.Get("Info");
38 TTree * treetoftrd =(TTree*)fff2.Get("TRD");
40 AliComparisonDraw comp2;
42 AliComparisonDraw comptof;
43 comptof->fTree =treetof;
55 #include "TBenchmark.h"
56 #include "TStopwatch.h"
57 #include "TParticle.h"
67 #include "TGeometry.h"
68 #include "TPolyLine3D.h"
74 #include "AliSimDigits.h"
75 #include "AliTPCParam.h"
77 #include "AliTPCLoader.h"
78 #include "AliDetector.h"
79 #include "AliTrackReference.h"
80 #include "AliTPCParamSR.h"
81 #include "AliTracker.h"
84 #include "AliPoints.h"
85 #include "AliESDtrack.h"
86 #include "AliTRDtrack.h"
87 #include "AliITStrackMI.h"
88 #include "AliESDV0MI.h"
89 #include "AliESDkink.h"
90 //#include "AliTRDparameter.h"
91 #include "AliTRDgeometryFull.h"
92 #include "AliTRDcluster.h"
93 #include "AliTRDtracker.h"
94 #include "AliTOFtrack.h"
95 #include "../TOF/AliTOFtrackerMI.h"
96 // COMPARISON includes
97 #include "../STEER/AliGenInfo.h"
98 #include "../STEER/AliESDComparisonMI.h"
102 AliComparisonDraw comp;
103 AliComparisonDraw comptr;
104 //AliTRDparameter par("pica","vyjebana");
105 AliTRDgeometryFull geom;
106 AliRunLoader *fLoader=0; //= AliRunLoader::Open(fnGalice);
108 TCut cnesd("cnesd","abs(RC.fLabel)!=MC.fLabel");
109 TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
110 TCut citsin("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<5");
111 TCut clength("clength","abs(TRD0.fIntegratedLength-MC.fTOFReferences[0].fLength)<30");
112 TCut c1("c1","TRD1.fStopped==0"); // track propagated to refernce
113 TCut c0("c0","TRD0.fStopped==0"); // track propagated to reference
114 TCut c2("c2","MC.fNTOFRef>0"); //track in tof
115 TCut c3("c3","MC.fNTOFRef>0&&MC.fNTRDRef>5"); //track in tof and TRD
116 TH1F hpull("hpull","hpull",50,-10,10);
117 TH1F hpull6("hpull6","hpull6",50,-6,6);
120 class AliTRDInfoMC: public TObject{
122 void Update(AliMCInfo *info); // update TRD info
123 void Reset(); // clear TRD info
124 Short_t fLayer[6]; //indicates layer is on
125 Short_t fLayerA[6]; // indicates layer after is ON
126 Short_t fNLayers; // number of crossed layers
127 Short_t fLast; // last crossed layer
128 AliTrackReference fRef[6]; // trd reference
129 AliTrackReference fRefLocal[6]; // trd local reference
130 ClassDef(AliTRDInfoMC,1) // TRD MC information
134 class AliTRDInfoCl: public TObject{
138 AliTRDcluster fCl[6][25]; // clusters
139 Float_t fCount[6][25]; // number of clusters of given label
140 Int_t fNClusters; // number of clusters
141 Float_t fMeanQ[6]; // mean Q
142 Int_t fMaxBin[6]; // max time bin
143 Int_t fNcl[6]; // number of clusters
144 Float_t fMeanQS[6]; // mean Q
145 Int_t fMaxBinS[6]; // max time bin
146 ClassDef(AliTRDInfoCl,1) // TRD MC information
149 ClassImp(AliTRDInfoCl)
153 class AliTRDInfoSeed: public TObject{
155 AliTRDInfoSeed(){fLabel=-1; fLabel1=-1;fLabel2=-1;}
156 void Update(AliTRDseed *seeds);
157 AliTRDseed fSeeds[6]; // seeds from seeding
158 Float_t fFakeRatio; //
159 Int_t fLabel; // seed label
160 Int_t fLabel1; // full track label
161 Int_t fLabel2; // second full track label
164 ClassDef(AliTRDInfoSeed,1) // TRD MC information
168 ClassImp(AliTRDInfoMC)
169 ClassImp(AliTRDInfoCl)
170 ClassImp(AliTRDInfoSeed)
172 void AliTRDInfoSeed::Update(AliTRDseed *seeds){
177 for (Int_t i=0;i<6;i++){
178 fSeeds[i] = seeds[i];
179 if (seeds[i].fN2>10) fNLayers++;
183 void AliTRDInfoMC::Reset(){
188 for (Int_t ilayer = 0; ilayer<6; ilayer++){
193 void AliTRDInfoMC::Update(AliMCInfo *info)
198 for (Int_t i=0;i<info->fNTRDRef;i++){
199 AliTrackReference * ref = (AliTrackReference*)info->fTRDReferences->At(i);
200 Float_t localX = ref->LocalX();
201 for (Int_t jplane = 0; jplane<6; jplane++){
202 if (TMath::Abs(localX - (300+jplane*14.))<6){
204 for (Int_t kplane=jplane;kplane>=0;kplane--){
208 fRefLocal[jplane].SetTrack(ref->GetTrack());
209 fRefLocal[jplane].SetPosition(ref->LocalX(), ref->LocalY(),ref->Z());
210 Float_t alpha = ref->Alpha();
211 Float_t dir[3] = {ref->Px(), ref->Py(), ref->Pz()};
213 dirR[0] = dir[0]*TMath::Cos(-alpha) - dir[1]*TMath::Sin(-alpha);
214 dirR[1] = dir[0]*TMath::Sin(-alpha) + dir[1]*TMath::Cos(-alpha);
215 fRefLocal[jplane].SetMomentum(1, dirR[1]/dirR[0],dir[2]/dir[0]);
220 for (Int_t jplane = 0; jplane<6; jplane++){
221 if (fLayer[jplane]>0) fNLayers++;
222 if (fLayer[jplane]>0) fLast = jplane;
226 AliTRDInfoCl::AliTRDInfoCl(){
230 for (Int_t ilayer=0;ilayer<6;ilayer++){
236 for (Int_t itime=0;itime<25;itime++){
237 fCount[ilayer][itime]=0;
243 void AliTRDInfoCl::Update(){
245 for (Int_t iplane=0;iplane<6;iplane++){
256 for (Int_t itime=0;itime<25;itime++){
257 if (fCl[iplane][itime].GetQ()>2){
259 fMeanQ[iplane]+=fCl[iplane][itime].GetQ();
260 if (fCl[iplane][itime].GetQ()>=maxQ){
261 fMaxBin[iplane]= itime;
262 maxQ = fCl[iplane][itime].GetQ();
265 if (fCl[iplane][itime].GetSumS()>2){
267 fMeanQS[iplane]+=fCl[iplane][itime].GetSumS();
268 if (fCl[iplane][itime].GetSumS()>=maxQS){
269 fMaxBinS[iplane]= itime;
270 maxQS = fCl[iplane][itime].GetSumS();
274 if (fNcl[iplane]>0) fMeanQ[iplane]/=Float_t(fNcl[iplane]);
275 if (sums>0) fMeanQS[iplane]/=Float_t(sums);
276 fNClusters+=fNcl[iplane];
283 void MaxBinNorm(const char * expr = "TRD0.fTimBinPlane"){
285 // max time bin probability distribution
288 sprintf(var,"%s>>hispi",expr);
289 TH1F * hispi= new TH1F("hispi","hispi",20,1,22);
290 comp.fTree->Draw(var,"abs(MC.fPdg)==211&&TRD0.fN>90&&MC.fParticle.P()>3");
291 hispi->Scale(1./(hispi->GetEntries()));
293 sprintf(var,"%s>>hisel",expr);
294 TH1F *hisel = new TH1F("hisel","hisel",20,1,22);
295 comp.fTree->Draw(var,"abs(MC.fPdg)==11&&TRD0.fN>90&&MC.fParticle.P()>3");
296 hisel->Scale(1./(hisel->GetEntries()));
297 hispi->SetMarkerStyle(25); hispi->SetLineStyle(2);
298 hisel->SetMarkerStyle(26); hisel->SetLineStyle(1);
299 hisel->SetXTitle("Maximal Time Bin");
300 hisel->SetYTitle("Probability");
303 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Maximal Time Bin position");
304 legend->SetBorderSize(1);
305 legend->AddEntry(hispi ,"Pions");
306 legend->AddEntry(hisel ,"Electrons");
312 void dedxNorm(const char * expr = "TRD0.fdEdxPlane/20"){
314 // max time bin probability distribution
317 sprintf(var,"%s>>hispi",expr);
318 TH1F * hispi= new TH1F("hispi","hispi",100,1,120);
319 comp.fTree->Draw(var,"abs(MC.fPdg)==211&&TRD0.fN>90&&MC.fParticle.P()>3");
320 hispi->Scale(1./(hispi->GetEntries()));
322 sprintf(var,"%s>>hisel",expr);
323 TH1F *hisel = new TH1F("hisel","hisel",100,1,120);
324 comp.fTree->Draw(var,"abs(MC.fPdg)==11&&TRD0.fN>90&&MC.fParticle.P()>3");
325 hisel->Scale(1./(hisel->GetEntries()));
326 hispi->SetMarkerStyle(25); hispi->SetLineStyle(2);
327 hisel->SetMarkerStyle(26); hisel->SetLineStyle(1);
328 hispi->SetXTitle("Normalized dEdx");
329 hispi->SetYTitle("Probability");
332 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Normalized dEdx");
333 legend->SetBorderSize(1);
334 legend->AddEntry(hispi ,"Pions");
335 legend->AddEntry(hisel ,"Electrons");
345 comp.Eff("TR.Pt()","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab>0",10,0.3,2);
346 TH1F * his = (TH1F*)comp.fRes->Clone();
347 comp.Eff("TR.Pt()","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab<0",10,0.3,2);
348 TH1F * hisf = (TH1F*)comp.fRes->Clone();
349 his->SetMarkerStyle(25); his->SetLineStyle(2);
350 hisf->SetMarkerStyle(26); hisf->SetLineStyle(1);
351 his->SetXTitle("P_{t}(GeV/c)");
352 his->SetYTitle("Tracking efficiency (%)");
355 TLegend *legend = new TLegend(0.55,0.16,0.85,0.4, "TRD efficciency");
356 legend->SetBorderSize(1);
357 legend->AddEntry(his ,"Efficiency");
358 legend->AddEntry(hisf ,"Fake tracks");
368 comp.Eff("atan(TR.fPz/TR.Pt())","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11&&TR.Pt()>0.4"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab>0",10,-0.9,0.9);
369 TH1F * his = (TH1F*)comp.fRes->Clone();
370 comp.Eff("atan(TR.fPz/TR.Pt())","InfoMC.fNLayers>3&&abs(MC.fPdg)!=11&&TR.Pt()>0.4"+c2+cprim,"TRD0.fN/InfoMC.fNLayers>10&&TRD0.fLab<0",10,-0.9,0.9);
371 TH1F * hisf = (TH1F*)comp.fRes->Clone();
372 his->SetMarkerStyle(25); his->SetLineStyle(2);
373 hisf->SetMarkerStyle(26); hisf->SetLineStyle(1);
374 his->SetXTitle("#lambda(rad)");
375 his->SetYTitle("Tracking efficiency (%)");
378 TLegend *legend = new TLegend(0.55,0.16,0.85,0.4, "TRD efficciency");
379 legend->SetBorderSize(1);
380 legend->AddEntry(his ,"Efficiency");
381 legend->AddEntry(hisf ,"Fake tracks");
387 TH1F *hisr = new TH1F("hisr","Percentage of found clusters",80,0,120);
388 comp.fTree->Draw("100*TRD0.fN/(20*InfoMC.fNLayers)>>hisr","InfoMC.fNLayers>3&&TRD0.fLab>0&&abs(MC.fPdg)!=11&&TR.Pt()>0.4"+cprim+c2);
389 hisr->SetXTitle("(%)");
395 TH1F * hderel = new TH1F("hderel","hderel",100,-150,150);
396 comp.fTree->Draw("100*derel>>hderel",""+c0+c2+cprim);
398 hderel->SetXTitle("Energy loss correction resolution (%)");
404 TH1F * hderel0 = new TH1F("hderel0","hderel0",100,-0,30);
405 TH1F * hderel1 = new TH1F("hderel1","hderel1",100,-0,30);
406 TH1F * hderela = new TH1F("hderela","hderela",100,-0,30);
407 comp.fTree->Draw("100*(TPCE-TOFE)/TPCE>>hderela","Pt<1"+c0+c2+cprim);
408 comp.fTree->Draw("100*(TPCE-TOFE)/TPCE>>hderel0","Pt<1&&TRD0.fBudget[0]<10"+c0+c2+cprim);
409 comp.fTree->Draw("100*(TPCE-TOFE)/TPCE>>hderel1","Pt<1&&TRD0.fBudget[0]>10"+c0+c2+cprim);
410 hderela->SetXTitle("(E_{TPC}-E_{TOF})/E_{TPC}(%)");
411 hderela->SetLineStyle(1);
412 hderel0->SetLineStyle(2);
413 hderel1->SetLineStyle(3);
415 hderel0->Draw("same");
416 hderel1->Draw("same");
417 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Relative Energy Loss Between TPC and TOF");
418 legend->SetBorderSize(1);
419 legend->AddEntry(hderela ,"All tracks");
420 legend->AddEntry(hderela ,"Non crossing of boundary");
421 legend->AddEntry(hderela ,"Crossing boundary");
430 TF1 f1("f1","100*exp(-[0]*x)",0,100);
431 comp.Eff("TRD.fBudget[0]","abs(ThetaM)<1&&abs(MC.fPdg)==211"+cprim,c2,10,1,100);
432 TH1F * heffb = (TH1F*)comp.fRes->Clone();
433 heffb->SetXTitle("Material budget (g/cm^{2})");
434 heffb->SetYTitle("Efficiency to reach TOF (%)");
442 TH1F * his = new TH1F("his","his",100,-100,10000);
443 TCut cut = "InfoTOF.fNCluster>1&&tofm0";
444 comp.fTree->Draw("InfoTOF.fTime1-33*sqrt(InfoTOF.fCl1.fZ**2+InfoTOF.fCl1.fR**2)>>his","InfoTOF.fTime1<100000"+cut);
445 his->SetXTitle("TOF Time (ns)");
449 TCut cut = "InfoTOF.fNCluster>0&&abs(MC.fPdg)==2212&&P<2";
451 TH1F * his0 = new TH1F("his0","his0",100,-1000,1000);
452 comp.fTree->Draw("InfoTOF.fTime0-InfoTOF.fTimes0[4]>>his0","tofm0"+cut);
453 TH1F * hisf = new TH1F("hisf","hisf",100,-1000,1000);
454 comp.fTree->Draw("InfoTOF.fTime0-InfoTOF.fTimes0[4]>>hisf","!(tofm0||tofm0s)"+cut);
457 his0->SetXTitle("TOF Time (ns)");
463 void EffPIDK(Float_t xmin=0.4, Float_t xmax=2){
466 //Float_t xmin=0.4; Float_t xmax=2;
469 // // TCut cut = "tofm"
470 // AliTOFtrackInfo::SetPIDMatchCuts();
471 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.1,0)==3",6,xmin,xmax);
472 // TH1F *his0 = (TH1F*)comp.fRes->Clone();
473 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.1,1)==3",6,xmin,xmax);
474 // TH1F *his1 = (TH1F*)comp.fRes->Clone();
475 // //AliTOFtrackInfo::SetPIDMatchCuts(-1,-1,-1,100); // disable pid match cuts
476 // //comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.0,0)==3",6,xmin,xmax);
477 // //TH1F *hisno = (TH1F*)comp.fRes->Clone();
478 // //AliTOFtrackInfo::SetPIDMatchCuts(); //enable again
480 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.1,2)==3",6,xmin,xmax);
481 // TH1F *hisold = (TH1F*)comp.fRes->Clone();
482 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.0,-1)==3",6,xmin,xmax);
483 // TH1F *histpc = (TH1F*)comp.fRes->Clone();
484 // his0->SetMarkerStyle(21);
485 // his1->SetMarkerStyle(24);
486 // // hisno->SetMarkerStyle(26);
487 // hisold->SetMarkerStyle(25);
488 // histpc->SetMarkerStyle(23);
491 // his1->Draw("same");
492 // //hisno->Draw("same");
493 // hisold->Draw("same");
494 // histpc->Draw("same");
495 // TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "PID efficiency ");
496 // legend->SetBorderSize(1);
497 // legend->AddEntry(his0 ,"TPC+New TOF v0");
498 // legend->AddEntry(his1 ,"TPC+New TOF v1");
499 // //legend->AddEntry(hisno ,"TPC+New TOF v0 - no PID match cuts");
500 // legend->AddEntry(hisold ,"TPC+Old TOF");
501 // legend->AddEntry(histpc ,"TPC only");
506 void EffPIDP(Float_t xmin=0.4, Float_t xmax=2){
509 //Float_t xmin=0.4; Float_t xmax=4;
513 //AliTOFtrackInfo::SetPIDMatchCuts(-1,-1,-1,1000); // default cuts
514 // AliTOFtrackInfo::SetPIDMatchCuts(); // default cuts
515 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,0)==4",6,xmin,xmax);
516 // TH1F *his0 = (TH1F*)comp.fRes->Clone();
517 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,1)==4",6,xmin,xmax);
518 // TH1F *his1 = (TH1F*)comp.fRes->Clone();
519 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,2)==4",6,xmin,xmax);
520 // TH1F *hisold = (TH1F*)comp.fRes->Clone();
521 // comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,-1)==4",6,xmin,xmax);
522 // TH1F *histpc = (TH1F*)comp.fRes->Clone();
523 // his0->SetMarkerStyle(21);
524 // his1->SetMarkerStyle(24);
525 // hisold->SetMarkerStyle(25);
526 // histpc->SetMarkerStyle(23);
528 // his1->Draw("same");
529 // hisold->Draw("same");
530 // histpc->Draw("same");
531 // TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "PID efficiency ");
532 // legend->SetBorderSize(1);
533 // legend->AddEntry(his0 ,"TPC+New TOF v0");
534 // legend->AddEntry(his1 ,"TPC+New TOF v1");
535 // legend->AddEntry(hisold ,"TPC+Old TOF");
536 // legend->AddEntry(histpc ,"TPC only");
552 TH1F * MakeCumul(TH1F *his,Bool_t norm = kTRUE)
554 TH1F *hcumul = (TH1F*)his->Clone();
556 sprintf(name,"N%f",gRandom->Rndm());
557 hcumul->SetName(name);
558 hcumul->SetTitle(name);
560 Float_t all = (Float_t)his->GetSum();
562 for (Int_t i=0;i<=his->GetNbinsX();i++){
563 hcumul->SetBinContent(i,his->Integral(0,i)/all);
565 hcumul->SetFillColor(0);
569 TGraph* His01(TH1F *his0,TH1F *his1){
570 Int_t nbins = his0->GetNbinsX();
571 Double_t *x = new Double_t[nbins];
572 Double_t *y = new Double_t[nbins];
573 for (Int_t i=0;i<nbins;i++) {
574 x[i] = his0->GetBinContent(i+1);
575 y[i] = his1->GetBinContent(i+1);
577 return new TGraph(nbins,x,y);
581 TH1F * PullCumul(const char *var, TCut cut, Int_t div=200, Float_t max=20){
582 TH1F * hpullb = new TH1F("hpullb","hpullb",div,0,max);
584 sprintf(v2,"%s>>hpullb",var);
585 comp.fTree->Draw(v2, cut);
586 TH1F * res = MakeCumul(hpullb);
595 // Set tree and set aliases
597 TFile *f = new TFile("trdComp1.root");
598 TTree * tree= (TTree*)f->Get("Comparison");
604 comp.fTree->SetAlias("Pt","sqrt(MC.fParticle.fPx**2+MC.fParticle.fPy**2)");
605 comp.fTree->SetAlias("P","sqrt(MC.fParticle.fPx**2+MC.fParticle.fPy**2+MC.fParticle.fPz**2)");
606 comp.fTree->SetAlias("Ptpc","sqrt(fTrackRef.fPx**2+fTrackRef.fPy**2+fTrackRef.fPz**2)");
607 comp.fTree->SetAlias("VRadius","sqrt(MC.fParticle.fVx**2+MC.fParticle.fVy**2)");
608 comp.fTree->SetAlias("DRadius","sqrt(MC.fTRdecay.fX**2+MC.fTRdecay.fY**2)");
609 comp.fTree->SetAlias("ThetaM","MC.fParticle.fPz/Pt");
610 comp.fTree->SetAlias("rel0","(MC.fTOFReferences[0].P()-MC.fTPCReferences[0].P())/MC.fTOFReferences[0].P()"); //relative en. loss
611 comp.fTree->SetAlias("rel1","(TR.P()-MC.fTPCReferences[0].P())/MC.fTPCReferences[0].P()"); //relative en. loss - in the last TRD
612 comp.fTree->SetAlias("pt_tof","MC.fTOFReferences[0].Pt()");
613 comp.fTree->SetAlias("dpt_tof","(abs(1/RC.fRp[4])-MC.fTOFReferences[0].Pt())/MC.fTOFReferences[0].Pt()");
616 comp.fTree->SetAlias("d1pt_tof","(abs(RC.fRp[4])-1/MC.fTOFReferences[0].Pt())");
617 comp.fTree->SetAlias("d1ptp_tof","(abs(RC.fRp[4])-1/MC.fTOFReferences[0].Pt())/sqrt(RC.fRc[14])");
619 comp.fTree->SetAlias("dpt_trd","(abs(1/RC.fRp[4])-TR.Pt())/TR.Pt()");
621 comp.fTree->SetAlias("tofth","(MC.fTOFReferences[0].fPz/MC.fTOFReferences[0].Pt())");
622 comp.fTree->SetAlias("dth_tof","TRD0.fT[3]-tofth");
623 comp.fTree->SetAlias("pullth0","(TRD0.fT[3]-tofth)/sqrt(TRD0.fCtt)");
624 comp.fTree->SetAlias("phi_tof","atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx)");
625 comp.fTree->SetAlias("phi_trd","asin(TRD0.GetSnp())+TRD0.fAlpha");
626 comp.fTree->SetAlias("pullphi0","(sin(asin(TRD0.GetSnp())+TRD0.fAlpha)-sin(atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx)))/sqrt(TRD0.fCee)");
627 comp.fTree->SetAlias("dphir0","(asin(TRD0.GetSnp())+TRD0.fAlpha-atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx))");
628 comp.fTree->SetAlias("dz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)");
629 comp.fTree->SetAlias("pully0","TRD0.fY/sqrt(TRD0.fCyy)");
630 comp.fTree->SetAlias("pullz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)/sqrt(TRD0.fCzz)");
632 comp.fTree->SetAlias("sigmay","sqrt(TRD0.fCyy)*(1+1/(1+abs(4./(RC.fRp[4]))))");
633 comp.fTree->SetAlias("pully3","TRD0.fY/sigmay");
634 comp.fTree->SetAlias("deltaz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)");
636 comp.fTree->SetAlias("pullc1","(abs(RC.fRp[4])-1/TR.Pt())/sqrt(RC.fRc[14])");
637 comp.fTree->SetAlias("deltac1","(abs(RC.fRp[4])-1/TR.Pt())*TR.Pt()");
638 comp.fTree->SetAlias("dist","sqrt(deltaz0**2+TRD0.fY**2)");
639 comp.fTree->SetAlias("pull0","pully0**2+pullz0**2");
640 comp.fTree->SetAlias("pullc0","(abs(TRD0.fC)-(1/(TRD0.fLocalConvConst*MC.fTOFReferences[0].Pt())))/sqrt(TRD0.fCcc)");
642 comp.fTree->SetAlias("prob5","exp(-(TRD1.fTracklets[5].fChi2)/16.)*max((TRD1.fTracklets[5].fNFound-6)/9.,0)");
643 comp.fTree->SetAlias("prob4","exp(-(TRD1.fTracklets[4].fChi2)/16.)*max((TRD1.fTracklets[4].fNFound-6)/9.,0)");
645 comp.fTree->SetAlias("TPCE","sqrt(MC.fTPCReferences[2].P()**2+MC.fMass**2)");
646 comp.fTree->SetAlias("TOFE","sqrt(MC.fTOFReferences[0].P()**2+MC.fMass**2)");
647 comp.fTree->SetAlias("TRDE","sqrt(TR.P()**2+MC.fMass**2)");
648 comp.fTree->SetAlias("dtpi","TRD0.fIntegratedTime[2]-10^12*MC.fTOFReferences.fTime");
649 comp.fTree->SetAlias("dtk","TRD0.fIntegratedTime[3]-10^12*MC.fTOFReferences.fTime");
650 comp.fTree->SetAlias("dtp","TRD0.fIntegratedTime[4]-10^12*MC.fTOFReferences.fTime");
651 comp.fTree->SetAlias("dt","(abs(MC.fPdg)==2212)*dtp+(abs(MC.fPdg)==211)*dtpi+(abs(MC.fPdg)==321)*dtk"); //delta time
652 comp.fTree->SetAlias("Beta","sqrt(TR.P()**2/(TR.P()**2+MC.fMass**2))");
653 comp.fTree->SetAlias("dtsigma","10+(1-Beta)*280");
654 comp.fTree->SetAlias("dptrel","(abs(TRD0.GetSignedPt())-MC.fTOFReferences[0].Pt())/MC.fTOFReferences[0].Pt()");
656 comp.fTree->SetAlias("dphi","atan2(TR.fY,TR..fX)-atan2(TR.fPy,TR..fPx)");
657 comp.fTree->SetAlias("dphi0","acos((TR.fX*TR.fPx+TR.fY*TR.fPy)/(TR.Pt()*TR.R()))");
658 comp.fTree->SetAlias("derel","(TRD0.fDE-(TPCE-TOFE))/(TPCE-TOFE)");
662 comp.fTree->SetAlias("tofm0","InfoTOF.fCl0.fLab[0]==abs(MC.fLabel)");
663 comp.fTree->SetAlias("tofm1","InfoTOF.fCl1.fLab[0]==abs(MC.fLabel)");
664 comp.fTree->SetAlias("tofm0s","InfoTOF.fCl0.fLab[0]>=abs(MC.fParticle.fDaughter[0]) && InfoTOF.fCl0.fLab[0]<=abs(MC.fParticle.fDaughter[1])");
665 comp.fTree->SetAlias("tofm1s","InfoTOF.fCl1.fLab[0]>=abs(MC.fParticle.fDaughter[0]) && InfoTOF.fCl1.fLab[0]<=abs(MC.fParticle.fDaughter[1])");
666 comp.fTree->SetAlias("tofm","tofm0||tofm1||tofm0s||tofm1s");
668 comp.fTree->SetAlias("Yref3","InfoMC.fSeeds[3].fYref[0]+InfoMC.fSeeds[3].fYref[1]*(InfoMC.fRef[3].LocalX()-InfoMC.fSeeds[3].fX0)");
669 comp.fTree->SetAlias("dyref3","Yref3-InfoMC.fRef[3].LocalY()");
671 comp.fTree->SetAlias("Yref5","InfoMC.fSeeds[5].fYref[0]+InfoMC.fSeeds[5].fYref[1]*(InfoMC.fRef[5].LocalX()-InfoMC.fSeeds[5].fX0)");
672 comp.fTree->SetAlias("dyref5","Yref5-InfoMC.fRef[5].LocalY()");
679 // Set tree and set aliases
681 TFile *ff = new TFile("TRDdebug.root");
682 TTree * treetr= (TTree*)ff->Get("tracklet");
683 comptr.fTree =treetr;
688 void DrawDE(TCut cut){
689 comp.DrawXY("(TRD0.fDE)","TPCE-TOFE",c2,cut,5,0.01,0.1,0,0.2);
696 void DrawYResol1(TCut cut){
697 comp.fTree->Draw("TRD1.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c1+cut,"");
700 void DrawYResol0(TCut cut){
701 comp.fTree->Draw("TRD0.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c0+cut,"");
705 void DrawTOFResY(Int_t ndiv =10,Float_t ptmin=0.5, Float_t ptmax =1.5, Float_t dymin=-0.5, Float_t dymax=0.5){
708 //Float_t ptmin=0.5, ptmax =1.5, dymin=-0.5, dymax=0.5;
709 //Float_t ptmin=0.5, ptmax =5, dymin=-0.3, dymax=0.3;
711 TCut cl("cl","InfoCl.fNClusters>60");
712 comp.DrawXY("TR.Pt()","TRD0.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
713 TH1F * histof = (TH1F*)comp.fRes->Clone();
714 comp.DrawXY("TR.Pt()","TRD1.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
715 TH1F * histrd = (TH1F*)comp.fRes->Clone();
716 histof->SetMarkerStyle(24);
717 histrd->SetMarkerStyle(25);
718 histof->SetXTitle("P_{t}(GeV)");
719 histof->SetYTitle("#Delta_{r-#phi}(cm)");
721 histrd->Draw("same");
722 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD resolution r-#phi (cm)");
723 legend->SetBorderSize(1);
724 legend->AddEntry(histof ,"Resolution in TOF plane");
725 legend->AddEntry(histrd ,"Resolution in last TRD plane");
729 void DrawTOFResZ(Float_t ptmin=0.5, Float_t ptmax =1.5, Float_t dzmin=-2.5, Float_t dzmax=2.5){
731 comp.DrawXY("TR.Pt()","TRD0.fZ-MC.fTOFReferences.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
732 TH1F * histof = (TH1F*)comp.fRes->Clone();
733 comp.DrawXY("TR.Pt()","TRD1.fZ-TR.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
734 TH1F * histrd = (TH1F*)comp.fRes->Clone();
735 histof->SetMarkerStyle(24);
736 histrd->SetMarkerStyle(25);
737 histof->SetXTitle("P_{t}(GeV)");
738 histof->SetYTitle("#Delta_{z}(cm)");
740 histrd->Draw("same");
741 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD resolution z (cm)");
742 legend->SetBorderSize(1);
743 legend->AddEntry(histof ,"Resolution in TOF plane");
744 legend->AddEntry(histrd ,"Resolution in last TRD plane");
749 void MakeCumul(TCut & cut)
751 // TCut cut ="1&&TR.Pt()>0.4";
752 TH1F * hisy = PullCumul("abs(pully0)",cut+cprim+c0+c2+c3,1000,20);
753 TH1F * hisz = PullCumul("abs(pullz0)",cut+cprim+c0+c2+c3,1000,20);
754 TH1F * hisyz = PullCumul("sqrt(pullz0**2+pully0**2)",cut+cprim+c0+c2+c3,1000,20);
755 hisz->SetXTitle("Pull (unit)");
756 hisz->SetYTitle("Cumulative function");
757 hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
758 hisz->Draw(); hisy->Draw("same");
759 //hisyz->Draw("same");
760 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD pull (unit)");
761 legend->SetBorderSize(1);
762 legend->AddEntry(hisy ,"Pull y");
763 legend->AddEntry(hisz,"Pull z");
764 // legend->AddEntry(hisyz,"Pull yz");
768 void MakeCumulAbs(TCut & cut)
771 TH1F * hisy = PullCumul("abs(TRD0.fY)",cut+cprim+c0+c2+c3,1000,6);
772 TH1F * hisz = PullCumul("abs(deltaz0)",cut+cprim+c0+c2+c3,1000,6);
773 TH1F * hisyz = PullCumul("sqrt(TRD0.fY**2+deltaz0**2)",cut+cprim+c0+c2+c3,1000,6);
774 hisz->SetXTitle("Delta (cm)");
775 hisz->SetYTitle("Cumulative function");
777 hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
778 hisz->Draw(); hisy->Draw("same");hisyz->Draw("same");
779 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD residulas (cm)");
780 legend->SetBorderSize(1);
781 legend->AddEntry(hisy ,"Delta y (cm)");
782 legend->AddEntry(hisz,"Delta z (cm)");
783 legend->AddEntry(hisyz,"Delta (cm)");
787 void LoadClusters(TObjArray *arr, Int_t event)
792 printf("Load cluster - event%d\n",event);
793 fLoader->SetEventNumber(event);
794 AliLoader * trdl = (AliLoader*)fLoader->GetLoader("TRDLoader");
795 trdl->LoadRecPoints();
797 TTree * ClusterTree = trdl->TreeR();
798 if (!ClusterTree) return;
799 TBranch * branch = ClusterTree->GetBranch("TRDcluster");
800 TObjArray *clusterArray = new TObjArray(1000);
801 branch->SetAddress(&clusterArray);
802 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
804 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
805 ClusterTree->GetEvent(iEntry);
806 for (Int_t icl =0; icl<clusterArray->GetEntriesFast();icl++){
807 AliTRDcluster *cl = (AliTRDcluster*)clusterArray->UncheckedAt(icl);
809 Int_t detector=cl->GetDetector();
810 Int_t localTimeBin=cl->GetLocalTimeBin();
811 // Int_t sector=geom.GetSector(detector);
812 Int_t plane=geom.GetPlane(detector);
813 for (Int_t ilab=0; ilab<3; ilab++){
814 Int_t label = cl->GetLabel(ilab);
816 if (label<0) continue;
817 if (localTimeBin>21||localTimeBin<0) continue;
818 if (plane<0||plane>5) continue;
820 AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(pos);
822 info = new AliTRDInfoCl;
823 arr->AddAt(info,pos);
826 if (info->fCount[plane][localTimeBin]>0){
827 if (TMath::Abs(info->fCl[plane][localTimeBin].GetY()-cl->GetY())<20 &&
828 TMath::Abs(info->fCl[plane][localTimeBin].GetZ()-cl->GetZ())<20){
829 info->fCount[plane][localTimeBin]++;
832 info->fCl[plane][localTimeBin]=*cl;
833 info->fCount[plane][localTimeBin]++;
838 for (Int_t i=0;i<arr->GetEntriesFast();i++){
839 AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(i);
840 if (info) info->Update();
846 void ReadSeeds(Int_t eventNr, TObjArray *sarray)
849 // read seeds form the file
852 AliTRDseed *pseeds[6]={&seeds[0],&seeds[1],&seeds[2],&seeds[3],&seeds[4],&seeds[5]};
853 Int_t label, label1, label2, seventNr, nused;
855 TFile f("TRDdebug.root");
856 TTree *tree =(TTree*)f.Get("Seeds2");
858 TBranch * br0 = tree->GetBranch("S0.");
859 TBranch * br1 = tree->GetBranch("S1.");
860 TBranch * br2 = tree->GetBranch("S2.");
861 TBranch * br3 = tree->GetBranch("S3.");
862 TBranch * br4 = tree->GetBranch("S4.");
863 TBranch * br5 = tree->GetBranch("S5.");
864 TBranch * brL = tree->GetBranch("Label");
865 TBranch * brL1 = tree->GetBranch("Label1");
866 TBranch * brL2 = tree->GetBranch("Label2");
867 TBranch * brE = tree->GetBranch("EventNr");
868 TBranch * brU = tree->GetBranch("NUsed");
869 TBranch * brFR = tree->GetBranch("FakeRatio");
870 brE->SetAddress(&seventNr);
871 brL->SetAddress(&label);
872 brL1->SetAddress(&label1);
873 brL2->SetAddress(&label2);
874 brU->SetAddress(&nused);
875 brFR->SetAddress(&fakeRatio);
876 br0->SetAddress(&pseeds[0]);
877 br1->SetAddress(&pseeds[1]);
878 br2->SetAddress(&pseeds[2]);
879 br3->SetAddress(&pseeds[3]);
880 br4->SetAddress(&pseeds[4]);
881 br5->SetAddress(&pseeds[5]);
884 for (Int_t ip=0;ip<sarray->GetEntriesFast();ip++){
885 AliTRDseed * seed = (AliTRDseed*)sarray->UncheckedAt(ip);
891 for (Int_t ip=0;ip<tree->GetEntries();ip++){
893 if (seventNr!=eventNr) continue;
894 //if (nused>15) continue;
895 if (!sarray->At(label)){
896 AliTRDInfoSeed * seed = new AliTRDInfoSeed;
898 seed->fLabel = label;
899 seed->fLabel1 = label1;
900 seed->fLabel2 = label2;
901 seed->fFakeRatio = fakeRatio;
902 seed->fNUsed = nused;
903 sarray->AddAt(seed,label);
908 void ReadTOFs(Int_t eventNr, TObjArray *tofs)
912 // AliTOFtrackInfo *infotof = new AliTOFtrackInfo;
913 // TFile ftof("TOFdebug.root");
914 // TTree * treetof =(TTree*)ftof.Get("Info");
915 // if (!treetof) return;
916 // TBranch * br30 = treetof->GetBranch("Info.");
917 // br30->SetAddress(&infotof);
920 // for (Int_t ip=0;ip<treetof->GetEntries();ip++){
921 // treetof->GetEntry(ip);
922 // if (eventNr!=infotof->fEventNr) continue;
923 // // if (eventNr>info->fEventNr) continue;
925 // Int_t label = TMath::Abs(infotof->fLab);
926 // if (label==0) continue;
927 // // if (label>=kmaxlabel) continue;
928 // if (!(tofs->At(label))){
929 // tofs->AddAt(new AliTOFtrackInfo(*infotof),label);
932 // if (infotof->fLength0>100)
933 // tofs->AddAt(new AliTOFtrackInfo(*infotof),label);
938 void ReadTracks(Int_t eventNr, TObjArray *esds, TObjArray *trds)
943 AliESDtrack * esdp = new AliESDtrack;
944 AliTRDtrack * trdp = new AliTRDtrack;
946 TFile f2("TRDdebug.root");
947 TTree * tree2 =(TTree*)f2.Get("Tracks");
949 TBranch * br20 = tree2->GetBranch("ESD.");
950 TBranch * br21 = tree2->GetBranch("trd.");
951 TBranch * br22 = tree2->GetBranch("EventNr");
952 br20->SetAddress(&esdp);
953 br21->SetAddress(&trdp);
954 br22->SetAddress(&teventNr);
957 for (Int_t ip=0;ip<esds->GetEntriesFast();ip++){
965 for (Int_t ip=0;ip<tree2->GetEntries();ip++){
967 if (eventNr!=teventNr) continue;
969 Int_t label = TMath::Abs(esdp->GetLabel());
970 if (label==0) continue;
971 if (!(esds->At(label))){
972 esds->AddAt(new AliESDtrack(*esdp),label);
973 trds->AddAt(new AliTRDtrack(*trdp),label);
975 AliTRDtrack * track =(AliTRDtrack*)trds->At(label);
976 if (track->GetNumberOfClusters() <trdp->GetNumberOfClusters()){
977 esds->AddAt(new AliESDtrack(*esdp),label);
978 trds->AddAt(new AliTRDtrack(*trdp),label);
988 fLoader = AliRunLoader::Open("galice.root");
989 AliESDtrack * esd = new AliESDtrack;
990 AliTRDtrack * ptrd = new AliTRDtrack;
991 AliTRDtrack * ptrd0 = new AliTRDtrack;
992 AliTRDtrack * ptrd1 = new AliTRDtrack;
993 AliTRDtrack * ptrdb = new AliTRDtrack;
995 AliTRDtrack * trd = ptrd;
996 AliTRDtrack * trd0 = ptrd0;
997 AliTRDtrack * trd1 = ptrd1;
998 AliTRDtrack * trdb = ptrdb;
1000 AliMCInfo* info = new AliMCInfo;
1001 AliTRDInfoMC *trdinfo = new AliTRDInfoMC;
1002 AliTRDInfoCl *infocl = new AliTRDInfoCl;
1003 AliTRDInfoSeed *infoSeed = new AliTRDInfoSeed;
1004 // AliTOFtrackInfo *infotof = new AliTOFtrackInfo;
1005 // AliTOFtrackInfo dummytof;
1006 // dummytof.fLab=-1;
1007 AliTRDInfoSeed dummySeed;
1008 AliTRDtrack dummytrd;
1009 AliTRDtrack dummytrd0;
1010 AliTRDtrack dummytrd1;
1011 AliTRDtrack dummytrdb;
1012 AliESDtrack dummyesd;
1013 dummySeed.fLabel =-1;
1014 AliTrackReference *ref = new AliTrackReference;
1015 TObjArray * clarray = new TObjArray(6*1000000);
1016 const Int_t kmaxlabel=1500000;
1020 TFile f("cmpESDTracks.root");
1021 TTree * tree1 = (TTree*) f.Get("ESDcmpTracks");
1022 TBranch * brmc = tree1->GetBranch("MC");
1023 brmc->SetAddress(&info);
1026 TObjArray *esds = new TObjArray(kmaxlabel);
1027 TObjArray *trds = new TObjArray(kmaxlabel);
1028 TObjArray *tofs = new TObjArray(kmaxlabel);
1029 TObjArray *seeds = new TObjArray(kmaxlabel);
1036 TFile f3("trdComp1.root","new");
1037 TTree * tree3 = new TTree("Comparison","Comparison");
1038 tree3->Branch("MC.","AliMCInfo",&info);
1039 tree3->Branch("RC.","AliESDtrack",&esd);
1040 tree3->Branch("TRD.","AliTRDtrack",&trd);
1041 tree3->Branch("TRD0.","AliTRDtrack",&trd0); // track at tof ref
1042 tree3->Branch("TRD1.","AliTRDtrack",&trd1); // track at last TRD ref
1043 tree3->Branch("TRDb.","AliTRDtrack",&trdb); // track at tof ref
1044 tree3->Branch("TR.","AliTrackReference",&ref); // track at last TRD ref
1045 tree3->Branch("InfoMC.","AliTRDInfoMC",&trdinfo); // info about TRD track MC
1046 tree3->Branch("InfoCl.","AliTRDInfoCl",&infocl); // info about TRD track MC
1047 tree3->Branch("InfoSeed.","AliTRDInfoSeed",&infoSeed); // info about TRD track MC
1048 // tree3->Branch("InfoTOF.","AliTOFtrackInfo",&infotof); // info about TOF
1053 for (Int_t i=0;i<tree1->GetEntries();i++){
1055 // if (info->fEventNr>5) break;
1056 if (lastevent!=info->fEventNr){
1057 printf("Read Event\t%d\n", info->fEventNr);
1058 ReadSeeds(info->fEventNr, seeds);
1059 ReadTOFs(info->fEventNr, tofs);
1060 ReadTracks(info->fEventNr,esds,trds);
1061 lastevent = info->fEventNr;
1063 LoadClusters(clarray,lastevent);
1066 // if (info->fNTOFRef==0 && info->fNTRDRef<5) continue;
1067 Int_t label = info->fLabel;
1068 trdinfo->Update(info);
1069 if (clarray->At(label)){
1070 infocl = (AliTRDInfoCl*)clarray->At(label);
1073 if (seeds->At(label)){
1074 infoSeed = (AliTRDInfoSeed*)seeds->At(label);
1076 infoSeed = &dummySeed;
1078 // if (info->fNTPCRef<3) continue;
1079 trd0->SetStop(kTRUE);
1080 trdb->SetStop(kTRUE);
1081 trd1->SetStop(kTRUE);
1082 // infotof = (AliTOFtrackInfo*)(tofs->At(label));
1084 // infotof = &dummytof;
1085 // infotof->fNCluster=0;
1086 // esd = (AliESDtrack*)(esds->At(label));
1089 // esd->GetTPCpid(pid);
1091 // for (Int_t ip=0;ip<5;ip++) suml+=pid[ip];
1092 // for (Int_t ip=0;ip<5;ip++) {
1093 // if (suml>0.0000000001) infotof->fRefPID[ip]= pid[ip]/suml;
1094 // else infotof->fRefPID[ip]=0.2;
1104 if (esds->At(label)) {
1108 // if (i%1000==0) printf("\n%d\n",i);
1109 trd = (AliTRDtrack*)(trds->At(label));
1110 esd = (AliESDtrack*)(esds->At(label));
1113 if (trd->GetBackupTrack())
1114 *trdb = *trd->GetBackupTrack();
1118 if (info->fNTOFRef>0){
1119 AliTrackReference * ref = ((AliTrackReference*)(info->fTOFReferences->At(0)));
1120 Double_t x = ((AliTrackReference*)(info->fTOFReferences->At(0)))->X();
1121 Double_t y = ((AliTrackReference*)(info->fTOFReferences->At(0)))->Y();
1122 Double_t phi = TMath::ATan2(y,x);
1123 Double_t radius = TMath::Sqrt(x*x+y*y);
1125 // rotate and propagate to local system
1127 trd0->SetStop(kFALSE);
1128 if (!trd0->Rotate(phi-trd0->GetAlpha())) trd0->SetStop(kTRUE);
1130 Double_t xyz0[3], param[7];
1131 trd0->GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
1132 Double_t xyz1[3]={ref->X(),ref->Y(),ref->Z()};
1133 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1134 if (!trd0->PropagateTo(radius,param[1],param[0]*param[4]))
1135 trd0->SetStop(kTRUE);
1137 trdb->SetStop(kFALSE);
1138 if (!trdb->Rotate(phi-trdb->GetAlpha())) trdb->SetStop(kTRUE);
1140 if (!trdb->PropagateTo(radius)) trdb->SetStop(kTRUE);
1143 if (info->fNTRDRef>0){
1146 for (Int_t itrd =0; itrd<info->fNTRDRef;itrd++){
1147 Double_t x = ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->X();
1148 Double_t y = ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Y();
1149 Double_t phi = TMath::ATan2(y,x);
1150 // Double_t snp = TMath::ATan2(((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Py(),
1151 // ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Px());
1152 // if (TMath::Abs(snp)>0.9 &&cradius>10) break;
1153 Double_t radius = TMath::Sqrt(x*x+y*y);
1154 if (radius>cradius){
1157 ref = ((AliTrackReference*)(info->fTRDReferences->At(itrd)));
1161 // rotate and propagate to local system of the last track reference
1163 trd1->SetStop(kFALSE);
1164 if (!trd1->Rotate(cphi-trd1->GetAlpha())) trd1->SetStop(kTRUE);
1166 if (!trd1->PropagateTo(cradius)) trd1->SetStop(kTRUE);