]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/Macros/AliTRDComparison.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
next50 trigger mask in AliHLTGlobalEsdConverterComponent
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDComparison.C
CommitLineData
9dc17dc7 1/*
2
3//load mag field
4AliRunLoader *loader = AliRunLoader::Open("galice.root");
5loader->LoadgAlice();
6gAlice = loader->GetAliRun();
7AliMagF * field = gAlice->Field();
e2cd3644 8AliTracker::SetFieldMap(gAlice->Field(),1);
9dc17dc7 9TFile fg("geometry.root");
10gGeoManager = (TGeoManager*)fg.Get("Geo");
7dd16671 11if (!gGeoManager) TGeoManager::Import("geometry.root");
9dc17dc7 12.L AliGenInfo.C+
13.L AliESDComparisonMI.C+
7dd16671 14.L AliTRDComparison.C+g
9dc17dc7 15//.L AliESDComparisonPIC.C+
16.x TDR_style.C
17
18MakeCompTr();
19MakeTree();
20
e2cd3644 21
7dd16671 22
9dc17dc7 23MakeComp();
24MakeAlias();
25TFile fff("TRDdebug.root");
26TTree * tree =(TTree*)fff.Get("tracklet");
27TTree * tree2 =(TTree*)fff.Get("Tracks");
7dd16671 28TTree * treefind =(TTree*)fff.Get("Find");
29TTree * treebug =(TTree*)fff.Get("Bug");
30TTree * treetofs =(TTree*)fff.Get("tofseed");
31AliComparisonDraw compfind;
32compfind.fTree = treefind;
33AliComparisonDraw comptr;
34comptr.fTree = tree;
9dc17dc7 35
7dd16671 36TFile fff2("TOFdebug.root");
37TTree * treetof =(TTree*)fff2.Get("Info");
38TTree * treetoftrd =(TTree*)fff2.Get("TRD");
9dc17dc7 39
7dd16671 40AliComparisonDraw comp2;
41comp2->fTree =tree;
42AliComparisonDraw comptof;
43comptof->fTree =treetof;
9dc17dc7 44
45
46*/
47
48//ROOT includes
49#include "Rtypes.h"
50#include "TFile.h"
51#include "TTree.h"
52#include "TChain.h"
53#include "TCut.h"
54#include "TString.h"
55#include "TBenchmark.h"
56#include "TStopwatch.h"
57#include "TParticle.h"
58#include "TSystem.h"
59#include "TTimer.h"
60#include "TVector3.h"
61#include "TH1F.h"
62#include "TH2F.h"
63#include "TCanvas.h"
64#include "TPad.h"
65#include "TF1.h"
66#include "TView.h"
67#include "TGeometry.h"
68#include "TPolyLine3D.h"
69#include "TLegend.h"
7dd16671 70#include "TGraph.h"
9dc17dc7 71//ALIROOT includes
72#include "AliRun.h"
73#include "AliStack.h"
74#include "AliSimDigits.h"
75#include "AliTPCParam.h"
76#include "AliTPC.h"
77#include "AliTPCLoader.h"
78#include "AliDetector.h"
79#include "AliTrackReference.h"
80#include "AliTPCParamSR.h"
81#include "AliTracker.h"
82#include "AliMagF.h"
83#include "AliHelix.h"
9dc17dc7 84#include "AliESDtrack.h"
85#include "AliTRDtrack.h"
86#include "AliITStrackMI.h"
87#include "AliESDV0MI.h"
88#include "AliESDkink.h"
7014e413 89//#include "AliTRDparameter.h"
9dc17dc7 90#include "AliTRDgeometryFull.h"
91#include "AliTRDcluster.h"
7dd16671 92#include "AliTRDtracker.h"
93#include "AliTOFtrack.h"
94#include "../TOF/AliTOFtrackerMI.h"
9dc17dc7 95// COMPARISON includes
96#include "../STEER/AliGenInfo.h"
97#include "../STEER/AliESDComparisonMI.h"
98
99
100
101AliComparisonDraw comp;
102AliComparisonDraw comptr;
7014e413 103//AliTRDparameter par("pica","vyjebana");
9dc17dc7 104AliTRDgeometryFull geom;
105AliRunLoader *fLoader=0; //= AliRunLoader::Open(fnGalice);
106
7dd16671 107TCut cnesd("cnesd","abs(RC.fLabel)!=MC.fLabel");
9dc17dc7 108TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
109TCut citsin("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<5");
110TCut clength("clength","abs(TRD0.fIntegratedLength-MC.fTOFReferences[0].fLength)<30");
111TCut c1("c1","TRD1.fStopped==0"); // track propagated to refernce
112TCut c0("c0","TRD0.fStopped==0"); // track propagated to reference
113TCut c2("c2","MC.fNTOFRef>0"); //track in tof
114TCut c3("c3","MC.fNTOFRef>0&&MC.fNTRDRef>5"); //track in tof and TRD
115TH1F hpull("hpull","hpull",50,-10,10);
116TH1F hpull6("hpull6","hpull6",50,-6,6);
117
118
119class AliTRDInfoMC: public TObject{
120 public:
121 void Update(AliMCInfo *info); // update TRD info
122 void Reset(); // clear TRD info
123 Short_t fLayer[6]; //indicates layer is on
124 Short_t fLayerA[6]; // indicates layer after is ON
125 Short_t fNLayers; // number of crossed layers
126 Short_t fLast; // last crossed layer
127 AliTrackReference fRef[6]; // trd reference
7dd16671 128 AliTrackReference fRefLocal[6]; // trd local reference
9dc17dc7 129 ClassDef(AliTRDInfoMC,1) // TRD MC information
130};
131
9dc17dc7 132
133class AliTRDInfoCl: public TObject{
134 public:
7dd16671 135 AliTRDInfoCl();
136 void Update();
137 AliTRDcluster fCl[6][25]; // clusters
138 Float_t fCount[6][25]; // number of clusters of given label
139 Int_t fNClusters; // number of clusters
140 Float_t fMeanQ[6]; // mean Q
141 Int_t fMaxBin[6]; // max time bin
142 Int_t fNcl[6]; // number of clusters
143 Float_t fMeanQS[6]; // mean Q
144 Int_t fMaxBinS[6]; // max time bin
145 ClassDef(AliTRDInfoCl,1) // TRD MC information
9dc17dc7 146};
147
148ClassImp(AliTRDInfoCl)
149
150
151
7dd16671 152class AliTRDInfoSeed: public TObject{
153 public:
154 AliTRDInfoSeed(){fLabel=-1; fLabel1=-1;fLabel2=-1;}
155 void Update(AliTRDseed *seeds);
156 AliTRDseed fSeeds[6]; // seeds from seeding
157 Float_t fFakeRatio; //
158 Int_t fLabel; // seed label
159 Int_t fLabel1; // full track label
160 Int_t fLabel2; // second full track label
161 Int_t fNLayers; //
162 Int_t fNUsed; //
163 ClassDef(AliTRDInfoSeed,1) // TRD MC information
164};
165
166
167ClassImp(AliTRDInfoMC)
168ClassImp(AliTRDInfoCl)
169ClassImp(AliTRDInfoSeed)
170
171void AliTRDInfoSeed::Update(AliTRDseed *seeds){
172 //
173 //
174 //
175 fNLayers=0;
176 for (Int_t i=0;i<6;i++){
177 fSeeds[i] = seeds[i];
178 if (seeds[i].fN2>10) fNLayers++;
179 }
180}
181
9dc17dc7 182void AliTRDInfoMC::Reset(){
183 //
184 //
185 fNLayers=0;
186 fLast =0;
187 for (Int_t ilayer = 0; ilayer<6; ilayer++){
188 fLayer[ilayer]=0;
189 fLayerA[ilayer]=0;
190 }
191}
192void AliTRDInfoMC::Update(AliMCInfo *info)
193{
194 //
195 //
196 Reset();
9dc17dc7 197 for (Int_t i=0;i<info->fNTRDRef;i++){
198 AliTrackReference * ref = (AliTrackReference*)info->fTRDReferences->At(i);
7dd16671 199 Float_t localX = ref->LocalX();
9dc17dc7 200 for (Int_t jplane = 0; jplane<6; jplane++){
7dd16671 201 if (TMath::Abs(localX - (300+jplane*14.))<6){
9dc17dc7 202 fLayer[jplane]++;
203 for (Int_t kplane=jplane;kplane>=0;kplane--){
204 fLayerA[kplane]++;
205 fRef[jplane]=*ref;
7dd16671 206 //
207 fRefLocal[jplane].SetTrack(ref->GetTrack());
208 fRefLocal[jplane].SetPosition(ref->LocalX(), ref->LocalY(),ref->Z());
209 Float_t alpha = ref->Alpha();
210 Float_t dir[3] = {ref->Px(), ref->Py(), ref->Pz()};
211 Float_t dirR[3];
212 dirR[0] = dir[0]*TMath::Cos(-alpha) - dir[1]*TMath::Sin(-alpha);
213 dirR[1] = dir[0]*TMath::Sin(-alpha) + dir[1]*TMath::Cos(-alpha);
214 fRefLocal[jplane].SetMomentum(1, dirR[1]/dirR[0],dir[2]/dir[0]);
9dc17dc7 215 }
216 }
217 }
218 }
219 for (Int_t jplane = 0; jplane<6; jplane++){
220 if (fLayer[jplane]>0) fNLayers++;
221 if (fLayer[jplane]>0) fLast = jplane;
222 }
223}
224
7dd16671 225AliTRDInfoCl::AliTRDInfoCl(){
226 //
227 //
228 fNClusters =0;
229 for (Int_t ilayer=0;ilayer<6;ilayer++){
230 fMeanQ[ilayer]=0;
231 fMaxBin[ilayer]=0;
232 fNcl[ilayer]=0;
233 fMeanQS[ilayer]=0;
234 fMaxBinS[ilayer]=0;
235 for (Int_t itime=0;itime<25;itime++){
236 fCount[ilayer][itime]=0;
237 }
238 }
239}
240
241
242void AliTRDInfoCl::Update(){
243 fNClusters = 0;
244 for (Int_t iplane=0;iplane<6;iplane++){
245 fMeanQ[iplane]=-1;
246 fNcl[iplane] = 0;
247 fMaxBin[iplane]=-1;
248 //
249 fMeanQS[iplane]=-1;
250 fMaxBinS[iplane]=-1;
251
252 Float_t maxQ =-10;
253 Float_t maxQS =-10;
254 Int_t sums =0;
255 for (Int_t itime=0;itime<25;itime++){
256 if (fCl[iplane][itime].GetQ()>2){
257 fNcl[iplane]++;
258 fMeanQ[iplane]+=fCl[iplane][itime].GetQ();
259 if (fCl[iplane][itime].GetQ()>=maxQ){
260 fMaxBin[iplane]= itime;
261 maxQ = fCl[iplane][itime].GetQ();
262 }
263 }
264 if (fCl[iplane][itime].GetSumS()>2){
265 sums++;
266 fMeanQS[iplane]+=fCl[iplane][itime].GetSumS();
267 if (fCl[iplane][itime].GetSumS()>=maxQS){
268 fMaxBinS[iplane]= itime;
269 maxQS = fCl[iplane][itime].GetSumS();
270 }
271 }
272 }
273 if (fNcl[iplane]>0) fMeanQ[iplane]/=Float_t(fNcl[iplane]);
274 if (sums>0) fMeanQS[iplane]/=Float_t(sums);
275 fNClusters+=fNcl[iplane];
276 }
277}
278
9dc17dc7 279
280
281
282void MaxBinNorm(const char * expr = "TRD0.fTimBinPlane"){
283 //
284 // max time bin probability distribution
285 //
286 char var[100];
287 sprintf(var,"%s>>hispi",expr);
288 TH1F * hispi= new TH1F("hispi","hispi",20,1,22);
289 comp.fTree->Draw(var,"abs(MC.fPdg)==211&&TRD0.fN>90&&MC.fParticle.P()>3");
290 hispi->Scale(1./(hispi->GetEntries()));
291 hispi->Draw();
292 sprintf(var,"%s>>hisel",expr);
293 TH1F *hisel = new TH1F("hisel","hisel",20,1,22);
294 comp.fTree->Draw(var,"abs(MC.fPdg)==11&&TRD0.fN>90&&MC.fParticle.P()>3");
295 hisel->Scale(1./(hisel->GetEntries()));
296 hispi->SetMarkerStyle(25); hispi->SetLineStyle(2);
297 hisel->SetMarkerStyle(26); hisel->SetLineStyle(1);
298 hisel->SetXTitle("Maximal Time Bin");
299 hisel->SetYTitle("Probability");
300 hisel->Draw();
301 hispi->Draw("same");
302 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Maximal Time Bin position");
303 legend->SetBorderSize(1);
304 legend->AddEntry(hispi ,"Pions");
305 legend->AddEntry(hisel ,"Electrons");
306 legend->Draw();
307
308}
309
310
311void dedxNorm(const char * expr = "TRD0.fdEdxPlane/20"){
312 //
313 // max time bin probability distribution
314 //
315 char var[100];
316 sprintf(var,"%s>>hispi",expr);
317 TH1F * hispi= new TH1F("hispi","hispi",100,1,120);
318 comp.fTree->Draw(var,"abs(MC.fPdg)==211&&TRD0.fN>90&&MC.fParticle.P()>3");
319 hispi->Scale(1./(hispi->GetEntries()));
320 hispi->Draw();
321 sprintf(var,"%s>>hisel",expr);
322 TH1F *hisel = new TH1F("hisel","hisel",100,1,120);
323 comp.fTree->Draw(var,"abs(MC.fPdg)==11&&TRD0.fN>90&&MC.fParticle.P()>3");
324 hisel->Scale(1./(hisel->GetEntries()));
325 hispi->SetMarkerStyle(25); hispi->SetLineStyle(2);
326 hisel->SetMarkerStyle(26); hisel->SetLineStyle(1);
327 hispi->SetXTitle("Normalized dEdx");
328 hispi->SetYTitle("Probability");
329 hispi->Draw();
330 hisel->Draw("same");
331 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Normalized dEdx");
332 legend->SetBorderSize(1);
333 legend->AddEntry(hispi ,"Pions");
334 legend->AddEntry(hisel ,"Electrons");
335 legend->Draw();
336}
337
338
339void EfficiencyPt()
340{
341 //
342 //
343 //
344 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);
345 TH1F * his = (TH1F*)comp.fRes->Clone();
346 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);
347 TH1F * hisf = (TH1F*)comp.fRes->Clone();
348 his->SetMarkerStyle(25); his->SetLineStyle(2);
349 hisf->SetMarkerStyle(26); hisf->SetLineStyle(1);
350 his->SetXTitle("P_{t}(GeV/c)");
351 his->SetYTitle("Tracking efficiency (%)");
352 his->Draw();
353 hisf->Draw("same");
354 TLegend *legend = new TLegend(0.55,0.16,0.85,0.4, "TRD efficciency");
355 legend->SetBorderSize(1);
356 legend->AddEntry(his ,"Efficiency");
357 legend->AddEntry(hisf ,"Fake tracks");
358 legend->Draw();
359
360}
361
362void EfficiencyDip()
363{
364 //
365 //
366 //
367 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);
368 TH1F * his = (TH1F*)comp.fRes->Clone();
369 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);
370 TH1F * hisf = (TH1F*)comp.fRes->Clone();
371 his->SetMarkerStyle(25); his->SetLineStyle(2);
372 hisf->SetMarkerStyle(26); hisf->SetLineStyle(1);
373 his->SetXTitle("#lambda(rad)");
374 his->SetYTitle("Tracking efficiency (%)");
375 his->Draw();
376 hisf->Draw("same");
377 TLegend *legend = new TLegend(0.55,0.16,0.85,0.4, "TRD efficciency");
378 legend->SetBorderSize(1);
379 legend->AddEntry(his ,"Efficiency");
380 legend->AddEntry(hisf ,"Fake tracks");
381 legend->Draw();
382}
383
384void ClusterRatio()
385{
386 TH1F *hisr = new TH1F("hisr","Percentage of found clusters",80,0,120);
387 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);
388 hisr->SetXTitle("(%)");
389 hisr->Draw();
390
391}
392
7dd16671 393void DEestimate(){
394 TH1F * hderel = new TH1F("hderel","hderel",100,-150,150);
395 comp.fTree->Draw("100*derel>>hderel",""+c0+c2+cprim);
396 hderel->Fit("gaus");
397 hderel->SetXTitle("Energy loss correction resolution (%)");
398 hderel->Draw();
399}
400
401void DEestimate2(){
402 //
403 TH1F * hderel0 = new TH1F("hderel0","hderel0",100,-0,30);
404 TH1F * hderel1 = new TH1F("hderel1","hderel1",100,-0,30);
405 TH1F * hderela = new TH1F("hderela","hderela",100,-0,30);
406 comp.fTree->Draw("100*(TPCE-TOFE)/TPCE>>hderela","Pt<1"+c0+c2+cprim);
407 comp.fTree->Draw("100*(TPCE-TOFE)/TPCE>>hderel0","Pt<1&&TRD0.fBudget[0]<10"+c0+c2+cprim);
408 comp.fTree->Draw("100*(TPCE-TOFE)/TPCE>>hderel1","Pt<1&&TRD0.fBudget[0]>10"+c0+c2+cprim);
409 hderela->SetXTitle("(E_{TPC}-E_{TOF})/E_{TPC}(%)");
410 hderela->SetLineStyle(1);
411 hderel0->SetLineStyle(2);
412 hderel1->SetLineStyle(3);
413 hderela->Draw("");
414 hderel0->Draw("same");
415 hderel1->Draw("same");
416 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "Relative Energy Loss Between TPC and TOF");
417 legend->SetBorderSize(1);
418 legend->AddEntry(hderela ,"All tracks");
419 legend->AddEntry(hderela ,"Non crossing of boundary");
420 legend->AddEntry(hderela ,"Crossing boundary");
421 legend->Draw();
422}
423
424
425void EffBudget(){
426 //
427 //
428 //
429 TF1 f1("f1","100*exp(-[0]*x)",0,100);
430 comp.Eff("TRD.fBudget[0]","abs(ThetaM)<1&&abs(MC.fPdg)==211"+cprim,c2,10,1,100);
431 TH1F * heffb = (TH1F*)comp.fRes->Clone();
432 heffb->SetXTitle("Material budget (g/cm^{2})");
433 heffb->SetYTitle("Efficiency to reach TOF (%)");
434 heffb->Fit(&f1);
435 heffb->Draw();
436
437
438}
439
440void TOFbackround(){
441 TH1F * his = new TH1F("his","his",100,-100,10000);
442 TCut cut = "InfoTOF.fNCluster>1&&tofm0";
443 comp.fTree->Draw("InfoTOF.fTime1-33*sqrt(InfoTOF.fCl1.fZ**2+InfoTOF.fCl1.fR**2)>>his","InfoTOF.fTime1<100000"+cut);
444 his->SetXTitle("TOF Time (ns)");
445}
446
447void TOFsignal(){
448 TCut cut = "InfoTOF.fNCluster>0&&abs(MC.fPdg)==2212&&P<2";
449 //
450 TH1F * his0 = new TH1F("his0","his0",100,-1000,1000);
451 comp.fTree->Draw("InfoTOF.fTime0-InfoTOF.fTimes0[4]>>his0","tofm0"+cut);
452 TH1F * hisf = new TH1F("hisf","hisf",100,-1000,1000);
453 comp.fTree->Draw("InfoTOF.fTime0-InfoTOF.fTimes0[4]>>hisf","!(tofm0||tofm0s)"+cut);
454 his0->Draw();
455 hisf->Draw("same");
456 his0->SetXTitle("TOF Time (ns)");
457}
458
459
460
461
462void EffPIDK(Float_t xmin=0.4, Float_t xmax=2){
463 //
464 //
465 //Float_t xmin=0.4; Float_t xmax=2;
466 // TCut cut= "1";
467// // TCut cut = c2;
468// // TCut cut = "tofm"
469// AliTOFtrackInfo::SetPIDMatchCuts();
470// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.1,0)==3",6,xmin,xmax);
471// TH1F *his0 = (TH1F*)comp.fRes->Clone();
472// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.1,1)==3",6,xmin,xmax);
473// TH1F *his1 = (TH1F*)comp.fRes->Clone();
474// //AliTOFtrackInfo::SetPIDMatchCuts(-1,-1,-1,100); // disable pid match cuts
475// //comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.0,0)==3",6,xmin,xmax);
476// //TH1F *hisno = (TH1F*)comp.fRes->Clone();
477// //AliTOFtrackInfo::SetPIDMatchCuts(); //enable again
478// //
479// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.1,2)==3",6,xmin,xmax);
480// TH1F *hisold = (TH1F*)comp.fRes->Clone();
481// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==321"+cut,"InfoTOF.GetCombPID(-1,0.0,-1)==3",6,xmin,xmax);
482// TH1F *histpc = (TH1F*)comp.fRes->Clone();
483// his0->SetMarkerStyle(21);
484// his1->SetMarkerStyle(24);
485// // hisno->SetMarkerStyle(26);
486// hisold->SetMarkerStyle(25);
487// histpc->SetMarkerStyle(23);
488// //
489// his0->Draw();
490// his1->Draw("same");
491// //hisno->Draw("same");
492// hisold->Draw("same");
493// histpc->Draw("same");
494// TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "PID efficiency ");
495// legend->SetBorderSize(1);
496// legend->AddEntry(his0 ,"TPC+New TOF v0");
497// legend->AddEntry(his1 ,"TPC+New TOF v1");
498// //legend->AddEntry(hisno ,"TPC+New TOF v0 - no PID match cuts");
499// legend->AddEntry(hisold ,"TPC+Old TOF");
500// legend->AddEntry(histpc ,"TPC only");
501// legend->Draw();
502 //
503}
504
505void EffPIDP(Float_t xmin=0.4, Float_t xmax=2){
506 //
507 //
508 //Float_t xmin=0.4; Float_t xmax=4;
509 TCut cut= "1";
510 // TCut cut= c2;
511 //TCut cut = "tofm0"
512 //AliTOFtrackInfo::SetPIDMatchCuts(-1,-1,-1,1000); // default cuts
513// AliTOFtrackInfo::SetPIDMatchCuts(); // default cuts
514// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,0)==4",6,xmin,xmax);
515// TH1F *his0 = (TH1F*)comp.fRes->Clone();
516// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,1)==4",6,xmin,xmax);
517// TH1F *his1 = (TH1F*)comp.fRes->Clone();
518// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,2)==4",6,xmin,xmax);
519// TH1F *hisold = (TH1F*)comp.fRes->Clone();
520// comp.Eff("P","TRD0.fIntegratedLength>360&&abs(MC.fPdg)==2212"+cut,"InfoTOF.GetCombPID(-1,0.1,-1)==4",6,xmin,xmax);
521// TH1F *histpc = (TH1F*)comp.fRes->Clone();
522// his0->SetMarkerStyle(21);
523// his1->SetMarkerStyle(24);
524// hisold->SetMarkerStyle(25);
525// histpc->SetMarkerStyle(23);
526// his0->Draw();
527// his1->Draw("same");
528// hisold->Draw("same");
529// histpc->Draw("same");
530// TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "PID efficiency ");
531// legend->SetBorderSize(1);
532// legend->AddEntry(his0 ,"TPC+New TOF v0");
533// legend->AddEntry(his1 ,"TPC+New TOF v1");
534// legend->AddEntry(hisold ,"TPC+Old TOF");
535// legend->AddEntry(histpc ,"TPC only");
536// legend->Draw();
537
538 //
539}
540
541
542
9dc17dc7 543
544
545void PtRes(){
546
547
548}
549
550
7dd16671 551TH1F * MakeCumul(TH1F *his,Bool_t norm = kTRUE)
9dc17dc7 552{
553 TH1F *hcumul = (TH1F*)his->Clone();
7dd16671 554 char name[1000];
555 sprintf(name,"N%f",gRandom->Rndm());
556 hcumul->SetName(name);
557 hcumul->SetTitle(name);
558 //
9dc17dc7 559 Float_t all = (Float_t)his->GetSum();
7dd16671 560 if (!norm) all =1;
9dc17dc7 561 for (Int_t i=0;i<=his->GetNbinsX();i++){
562 hcumul->SetBinContent(i,his->Integral(0,i)/all);
563 }
564 hcumul->SetFillColor(0);
565 return hcumul;
566}
567
7dd16671 568TGraph* His01(TH1F *his0,TH1F *his1){
569 Int_t nbins = his0->GetNbinsX();
570 Double_t *x = new Double_t[nbins];
571 Double_t *y = new Double_t[nbins];
572 for (Int_t i=0;i<nbins;i++) {
573 x[i] = his0->GetBinContent(i+1);
574 y[i] = his1->GetBinContent(i+1);
575 }
576 return new TGraph(nbins,x,y);
577}
578
579
9dc17dc7 580TH1F * PullCumul(const char *var, TCut cut, Int_t div=200, Float_t max=20){
581 TH1F * hpullb = new TH1F("hpullb","hpullb",div,0,max);
582 char v2[1000];
583 sprintf(v2,"%s>>hpullb",var);
584 comp.fTree->Draw(v2, cut);
585 TH1F * res = MakeCumul(hpullb);
586 delete hpullb;
587 return res;
588}
589
590
591
592void MakeComp(){
593 //
594 // Set tree and set aliases
595 //
596 TFile *f = new TFile("trdComp1.root");
597 TTree * tree= (TTree*)f->Get("Comparison");
598 comp.fTree =tree;
599}
600
601void MakeAlias(){
602
603 comp.fTree->SetAlias("Pt","sqrt(MC.fParticle.fPx**2+MC.fParticle.fPy**2)");
604 comp.fTree->SetAlias("P","sqrt(MC.fParticle.fPx**2+MC.fParticle.fPy**2+MC.fParticle.fPz**2)");
605 comp.fTree->SetAlias("Ptpc","sqrt(fTrackRef.fPx**2+fTrackRef.fPy**2+fTrackRef.fPz**2)");
606 comp.fTree->SetAlias("VRadius","sqrt(MC.fParticle.fVx**2+MC.fParticle.fVy**2)");
607 comp.fTree->SetAlias("DRadius","sqrt(MC.fTRdecay.fX**2+MC.fTRdecay.fY**2)");
608 comp.fTree->SetAlias("ThetaM","MC.fParticle.fPz/Pt");
609 comp.fTree->SetAlias("rel0","(MC.fTOFReferences[0].P()-MC.fTPCReferences[0].P())/MC.fTOFReferences[0].P()"); //relative en. loss
610 comp.fTree->SetAlias("rel1","(TR.P()-MC.fTPCReferences[0].P())/MC.fTPCReferences[0].P()"); //relative en. loss - in the last TRD
611 comp.fTree->SetAlias("pt_tof","MC.fTOFReferences[0].Pt()");
612 comp.fTree->SetAlias("dpt_tof","(abs(1/RC.fRp[4])-MC.fTOFReferences[0].Pt())/MC.fTOFReferences[0].Pt()");
613
614
615 comp.fTree->SetAlias("d1pt_tof","(abs(RC.fRp[4])-1/MC.fTOFReferences[0].Pt())");
616 comp.fTree->SetAlias("d1ptp_tof","(abs(RC.fRp[4])-1/MC.fTOFReferences[0].Pt())/sqrt(RC.fRc[14])");
617
618 comp.fTree->SetAlias("dpt_trd","(abs(1/RC.fRp[4])-TR.Pt())/TR.Pt()");
619
620 comp.fTree->SetAlias("tofth","(MC.fTOFReferences[0].fPz/MC.fTOFReferences[0].Pt())");
621 comp.fTree->SetAlias("dth_tof","TRD0.fT[3]-tofth");
7dd16671 622 comp.fTree->SetAlias("pullth0","(TRD0.fT[3]-tofth)/sqrt(TRD0.fCtt)");
9dc17dc7 623 comp.fTree->SetAlias("phi_tof","atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx)");
624 comp.fTree->SetAlias("phi_trd","asin(TRD0.GetSnp())+TRD0.fAlpha");
7dd16671 625 comp.fTree->SetAlias("pullphi0","(sin(asin(TRD0.GetSnp())+TRD0.fAlpha)-sin(atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx)))/sqrt(TRD0.fCee)");
626 comp.fTree->SetAlias("dphir0","(asin(TRD0.GetSnp())+TRD0.fAlpha-atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx))");
9dc17dc7 627 comp.fTree->SetAlias("dz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)");
628 comp.fTree->SetAlias("pully0","TRD0.fY/sqrt(TRD0.fCyy)");
629 comp.fTree->SetAlias("pullz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)/sqrt(TRD0.fCzz)");
630 //
631 comp.fTree->SetAlias("sigmay","sqrt(TRD0.fCyy)*(1+1/(1+abs(4./(RC.fRp[4]))))");
632 comp.fTree->SetAlias("pully3","TRD0.fY/sigmay");
633 comp.fTree->SetAlias("deltaz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)");
634 //
635 comp.fTree->SetAlias("pullc1","(abs(RC.fRp[4])-1/TR.Pt())/sqrt(RC.fRc[14])");
636 comp.fTree->SetAlias("deltac1","(abs(RC.fRp[4])-1/TR.Pt())*TR.Pt()");
637 comp.fTree->SetAlias("dist","sqrt(deltaz0**2+TRD0.fY**2)");
638 comp.fTree->SetAlias("pull0","pully0**2+pullz0**2");
639 comp.fTree->SetAlias("pullc0","(abs(TRD0.fC)-(1/(TRD0.fLocalConvConst*MC.fTOFReferences[0].Pt())))/sqrt(TRD0.fCcc)");
640 //
641 comp.fTree->SetAlias("prob5","exp(-(TRD1.fTracklets[5].fChi2)/16.)*max((TRD1.fTracklets[5].fNFound-6)/9.,0)");
642 comp.fTree->SetAlias("prob4","exp(-(TRD1.fTracklets[4].fChi2)/16.)*max((TRD1.fTracklets[4].fNFound-6)/9.,0)");
643 //
e2cd3644 644 comp.fTree->SetAlias("TPCE","sqrt(MC.fTPCReferences[2].P()**2+MC.fMass**2)");
9dc17dc7 645 comp.fTree->SetAlias("TOFE","sqrt(MC.fTOFReferences[0].P()**2+MC.fMass**2)");
646 comp.fTree->SetAlias("TRDE","sqrt(TR.P()**2+MC.fMass**2)");
647 comp.fTree->SetAlias("dtpi","TRD0.fIntegratedTime[2]-10^12*MC.fTOFReferences.fTime");
648 comp.fTree->SetAlias("dtk","TRD0.fIntegratedTime[3]-10^12*MC.fTOFReferences.fTime");
649 comp.fTree->SetAlias("dtp","TRD0.fIntegratedTime[4]-10^12*MC.fTOFReferences.fTime");
650 comp.fTree->SetAlias("dt","(abs(MC.fPdg)==2212)*dtp+(abs(MC.fPdg)==211)*dtpi+(abs(MC.fPdg)==321)*dtk"); //delta time
651 comp.fTree->SetAlias("Beta","sqrt(TR.P()**2/(TR.P()**2+MC.fMass**2))");
652 comp.fTree->SetAlias("dtsigma","10+(1-Beta)*280");
6c23ffed 653 comp.fTree->SetAlias("dptrel","(abs(TRD0.GetSignedPt())-MC.fTOFReferences[0].Pt())/MC.fTOFReferences[0].Pt()");
9dc17dc7 654 //
655 comp.fTree->SetAlias("dphi","atan2(TR.fY,TR..fX)-atan2(TR.fPy,TR..fPx)");
656 comp.fTree->SetAlias("dphi0","acos((TR.fX*TR.fPx+TR.fY*TR.fPy)/(TR.Pt()*TR.R()))");
7dd16671 657 comp.fTree->SetAlias("derel","(TRD0.fDE-(TPCE-TOFE))/(TPCE-TOFE)");
658 //
659 // TOF stuff
660 //
661 comp.fTree->SetAlias("tofm0","InfoTOF.fCl0.fLab[0]==abs(MC.fLabel)");
662 comp.fTree->SetAlias("tofm1","InfoTOF.fCl1.fLab[0]==abs(MC.fLabel)");
663 comp.fTree->SetAlias("tofm0s","InfoTOF.fCl0.fLab[0]>=abs(MC.fParticle.fDaughter[0]) && InfoTOF.fCl0.fLab[0]<=abs(MC.fParticle.fDaughter[1])");
664 comp.fTree->SetAlias("tofm1s","InfoTOF.fCl1.fLab[0]>=abs(MC.fParticle.fDaughter[0]) && InfoTOF.fCl1.fLab[0]<=abs(MC.fParticle.fDaughter[1])");
665 comp.fTree->SetAlias("tofm","tofm0||tofm1||tofm0s||tofm1s");
666
667 comp.fTree->SetAlias("Yref3","InfoMC.fSeeds[3].fYref[0]+InfoMC.fSeeds[3].fYref[1]*(InfoMC.fRef[3].LocalX()-InfoMC.fSeeds[3].fX0)");
668 comp.fTree->SetAlias("dyref3","Yref3-InfoMC.fRef[3].LocalY()");
9dc17dc7 669
7dd16671 670 comp.fTree->SetAlias("Yref5","InfoMC.fSeeds[5].fYref[0]+InfoMC.fSeeds[5].fYref[1]*(InfoMC.fRef[5].LocalX()-InfoMC.fSeeds[5].fX0)");
671 comp.fTree->SetAlias("dyref5","Yref5-InfoMC.fRef[5].LocalY()");
9dc17dc7 672
673}
674
675
676void MakeCompTr(){
677 //
678 // Set tree and set aliases
679 //
680 TFile *ff = new TFile("TRDdebug.root");
681 TTree * treetr= (TTree*)ff->Get("tracklet");
682 comptr.fTree =treetr;
683
684}
685
686
687void DrawDE(TCut cut){
688 comp.DrawXY("(TRD0.fDE)","TPCE-TOFE",c2,cut,5,0.01,0.1,0,0.2);
689 comp.fRes->Draw();
690}
691
692
7dd16671 693
694
9dc17dc7 695void DrawYResol1(TCut cut){
696 comp.fTree->Draw("TRD1.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c1+cut,"");
697}
698
699void DrawYResol0(TCut cut){
700 comp.fTree->Draw("TRD0.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c0+cut,"");
701}
702
703
e2cd3644 704void 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){
9dc17dc7 705 //
7dd16671 706 //
707 //Float_t ptmin=0.5, ptmax =1.5, dymin=-0.5, dymax=0.5;
708 //Float_t ptmin=0.5, ptmax =5, dymin=-0.3, dymax=0.3;
e2cd3644 709 // Int_t ndiv =5;
7dd16671 710 TCut cl("cl","InfoCl.fNClusters>60");
711 comp.DrawXY("TR.Pt()","TRD0.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
9dc17dc7 712 TH1F * histof = (TH1F*)comp.fRes->Clone();
7dd16671 713 comp.DrawXY("TR.Pt()","TRD1.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
9dc17dc7 714 TH1F * histrd = (TH1F*)comp.fRes->Clone();
715 histof->SetMarkerStyle(24);
716 histrd->SetMarkerStyle(25);
717 histof->SetXTitle("P_{t}(GeV)");
718 histof->SetYTitle("#Delta_{r-#phi}(cm)");
719 histof->Draw();
720 histrd->Draw("same");
721 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD resolution r-#phi (cm)");
722 legend->SetBorderSize(1);
723 legend->AddEntry(histof ,"Resolution in TOF plane");
724 legend->AddEntry(histrd ,"Resolution in last TRD plane");
725 legend->Draw();
726}
727
7dd16671 728void DrawTOFResZ(Float_t ptmin=0.5, Float_t ptmax =1.5, Float_t dzmin=-2.5, Float_t dzmax=2.5){
9dc17dc7 729 //
7dd16671 730 comp.DrawXY("TR.Pt()","TRD0.fZ-MC.fTOFReferences.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
9dc17dc7 731 TH1F * histof = (TH1F*)comp.fRes->Clone();
7dd16671 732 comp.DrawXY("TR.Pt()","TRD1.fZ-TR.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
9dc17dc7 733 TH1F * histrd = (TH1F*)comp.fRes->Clone();
734 histof->SetMarkerStyle(24);
735 histrd->SetMarkerStyle(25);
736 histof->SetXTitle("P_{t}(GeV)");
737 histof->SetYTitle("#Delta_{z}(cm)");
738 histof->Draw();
739 histrd->Draw("same");
740 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD resolution z (cm)");
741 legend->SetBorderSize(1);
742 legend->AddEntry(histof ,"Resolution in TOF plane");
743 legend->AddEntry(histrd ,"Resolution in last TRD plane");
744 legend->Draw();
745}
746
747
748void MakeCumul(TCut & cut)
749{
7dd16671 750 // TCut cut ="1&&TR.Pt()>0.4";
9dc17dc7 751 TH1F * hisy = PullCumul("abs(pully0)",cut+cprim+c0+c2+c3,1000,20);
752 TH1F * hisz = PullCumul("abs(pullz0)",cut+cprim+c0+c2+c3,1000,20);
753 TH1F * hisyz = PullCumul("sqrt(pullz0**2+pully0**2)",cut+cprim+c0+c2+c3,1000,20);
754 hisz->SetXTitle("Pull (unit)");
755 hisz->SetYTitle("Cumulative function");
7dd16671 756 hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
757 hisz->Draw(); hisy->Draw("same");
758 //hisyz->Draw("same");
9dc17dc7 759 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD pull (unit)");
760 legend->SetBorderSize(1);
761 legend->AddEntry(hisy ,"Pull y");
762 legend->AddEntry(hisz,"Pull z");
7dd16671 763 // legend->AddEntry(hisyz,"Pull yz");
9dc17dc7 764 legend->Draw();
765}
766
767void MakeCumulAbs(TCut & cut)
768{
7dd16671 769 // TCut cut = "1";
9dc17dc7 770 TH1F * hisy = PullCumul("abs(TRD0.fY)",cut+cprim+c0+c2+c3,1000,6);
771 TH1F * hisz = PullCumul("abs(deltaz0)",cut+cprim+c0+c2+c3,1000,6);
772 TH1F * hisyz = PullCumul("sqrt(TRD0.fY**2+deltaz0**2)",cut+cprim+c0+c2+c3,1000,6);
773 hisz->SetXTitle("Delta (cm)");
774 hisz->SetYTitle("Cumulative function");
775 //
7dd16671 776 hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
9dc17dc7 777 hisz->Draw(); hisy->Draw("same");hisyz->Draw("same");
778 TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD residulas (cm)");
779 legend->SetBorderSize(1);
780 legend->AddEntry(hisy ,"Delta y (cm)");
781 legend->AddEntry(hisz,"Delta z (cm)");
782 legend->AddEntry(hisyz,"Delta (cm)");
783 legend->Draw();
784}
785
786void LoadClusters(TObjArray *arr, Int_t event)
787{
788 //
789 //
790 // event=0
791 printf("Load cluster - event%d\n",event);
792 fLoader->SetEventNumber(event);
793 AliLoader * trdl = (AliLoader*)fLoader->GetLoader("TRDLoader");
794 trdl->LoadRecPoints();
795
796 TTree * ClusterTree = trdl->TreeR();
797 if (!ClusterTree) return;
798 TBranch * branch = ClusterTree->GetBranch("TRDcluster");
799 TObjArray *clusterArray = new TObjArray(1000);
800 branch->SetAddress(&clusterArray);
801 Int_t nEntries = (Int_t) ClusterTree->GetEntries();
802 //
803 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
804 ClusterTree->GetEvent(iEntry);
805 for (Int_t icl =0; icl<clusterArray->GetEntriesFast();icl++){
806 AliTRDcluster *cl = (AliTRDcluster*)clusterArray->UncheckedAt(icl);
807 if (!cl) continue;
9dc17dc7 808 Int_t detector=cl->GetDetector();
809 Int_t localTimeBin=cl->GetLocalTimeBin();
810 // Int_t sector=geom.GetSector(detector);
811 Int_t plane=geom.GetPlane(detector);
7dd16671 812 for (Int_t ilab=0; ilab<3; ilab++){
813 Int_t label = cl->GetLabel(ilab);
814 Int_t pos = label;
815 if (label<0) continue;
816 if (localTimeBin>21||localTimeBin<0) continue;
817 if (plane<0||plane>5) continue;
818
819 AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(pos);
820 if (!info) {
821 info = new AliTRDInfoCl;
822 arr->AddAt(info,pos);
823 }
824 //
825 if (info->fCount[plane][localTimeBin]>0){
826 if (TMath::Abs(info->fCl[plane][localTimeBin].GetY()-cl->GetY())<20 &&
827 TMath::Abs(info->fCl[plane][localTimeBin].GetZ()-cl->GetZ())<20){
828 info->fCount[plane][localTimeBin]++;
829 }
830 }else{
831 info->fCl[plane][localTimeBin]=*cl;
832 info->fCount[plane][localTimeBin]++;
833 }
834 }
9dc17dc7 835 }
836 }
7dd16671 837 for (Int_t i=0;i<arr->GetEntriesFast();i++){
838 AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(i);
839 if (info) info->Update();
840 }
9dc17dc7 841}
842
843
7dd16671 844
845void ReadSeeds(Int_t eventNr, TObjArray *sarray)
846{
847 //
848 // read seeds form the file
849 //
850 AliTRDseed seeds[6];
851 AliTRDseed *pseeds[6]={&seeds[0],&seeds[1],&seeds[2],&seeds[3],&seeds[4],&seeds[5]};
852 Int_t label, label1, label2, seventNr, nused;
853 Float_t fakeRatio;
854 TFile f("TRDdebug.root");
855 TTree *tree =(TTree*)f.Get("Seeds2");
856 if (!tree) return;
857 TBranch * br0 = tree->GetBranch("S0.");
858 TBranch * br1 = tree->GetBranch("S1.");
859 TBranch * br2 = tree->GetBranch("S2.");
860 TBranch * br3 = tree->GetBranch("S3.");
861 TBranch * br4 = tree->GetBranch("S4.");
862 TBranch * br5 = tree->GetBranch("S5.");
863 TBranch * brL = tree->GetBranch("Label");
864 TBranch * brL1 = tree->GetBranch("Label1");
865 TBranch * brL2 = tree->GetBranch("Label2");
866 TBranch * brE = tree->GetBranch("EventNr");
867 TBranch * brU = tree->GetBranch("NUsed");
868 TBranch * brFR = tree->GetBranch("FakeRatio");
869 brE->SetAddress(&seventNr);
870 brL->SetAddress(&label);
871 brL1->SetAddress(&label1);
872 brL2->SetAddress(&label2);
873 brU->SetAddress(&nused);
874 brFR->SetAddress(&fakeRatio);
875 br0->SetAddress(&pseeds[0]);
876 br1->SetAddress(&pseeds[1]);
877 br2->SetAddress(&pseeds[2]);
878 br3->SetAddress(&pseeds[3]);
879 br4->SetAddress(&pseeds[4]);
880 br5->SetAddress(&pseeds[5]);
881 //
882 // delete all seeds
883 for (Int_t ip=0;ip<sarray->GetEntriesFast();ip++){
884 AliTRDseed * seed = (AliTRDseed*)sarray->UncheckedAt(ip);
885 if (seed){
886 delete seed;
887 }
888 }
889 sarray->Clear();
890 for (Int_t ip=0;ip<tree->GetEntries();ip++){
891 tree->GetEntry(ip);
892 if (seventNr!=eventNr) continue;
893 //if (nused>15) continue;
894 if (!sarray->At(label)){
895 AliTRDInfoSeed * seed = new AliTRDInfoSeed;
896 seed->Update(seeds);
897 seed->fLabel = label;
898 seed->fLabel1 = label1;
899 seed->fLabel2 = label2;
900 seed->fFakeRatio = fakeRatio;
901 seed->fNUsed = nused;
902 sarray->AddAt(seed,label);
903 }
904 }
905}
906
907void ReadTOFs(Int_t eventNr, TObjArray *tofs)
908{
909 //
910 // tofs->Clear();
911// AliTOFtrackInfo *infotof = new AliTOFtrackInfo;
912// TFile ftof("TOFdebug.root");
913// TTree * treetof =(TTree*)ftof.Get("Info");
914// if (!treetof) return;
915// TBranch * br30 = treetof->GetBranch("Info.");
916// br30->SetAddress(&infotof);
917// //
918// // read tof info
919// for (Int_t ip=0;ip<treetof->GetEntries();ip++){
920// treetof->GetEntry(ip);
921// if (eventNr!=infotof->fEventNr) continue;
922// // if (eventNr>info->fEventNr) continue;
923// //
924// Int_t label = TMath::Abs(infotof->fLab);
925// if (label==0) continue;
926// // if (label>=kmaxlabel) continue;
927// if (!(tofs->At(label))){
928// tofs->AddAt(new AliTOFtrackInfo(*infotof),label);
929// }
930// else{
931// if (infotof->fLength0>100)
932// tofs->AddAt(new AliTOFtrackInfo(*infotof),label);
933// }
934// }
935}
936
937void ReadTracks(Int_t eventNr, TObjArray *esds, TObjArray *trds)
938{
939 //
940 // read debug tracks
941 //
942 AliESDtrack * esdp = new AliESDtrack;
943 AliTRDtrack * trdp = new AliTRDtrack;
944 Int_t teventNr;
945 TFile f2("TRDdebug.root");
946 TTree * tree2 =(TTree*)f2.Get("Tracks");
947 if (!tree2) return;
948 TBranch * br20 = tree2->GetBranch("ESD.");
949 TBranch * br21 = tree2->GetBranch("trd.");
950 TBranch * br22 = tree2->GetBranch("EventNr");
951 br20->SetAddress(&esdp);
952 br21->SetAddress(&trdp);
953 br22->SetAddress(&teventNr);
954 //
955 // delete all event
956 for (Int_t ip=0;ip<esds->GetEntriesFast();ip++){
957 if (esds->At(ip)){
958 delete esds->At(ip);
959 delete trds->At(ip);
960 }
961 }
962 esds->Clear();
963 trds->Clear();
964 for (Int_t ip=0;ip<tree2->GetEntries();ip++){
965 tree2->GetEntry(ip);
966 if (eventNr!=teventNr) continue;
967 //
e2cd3644 968 Int_t label = TMath::Abs(esdp->GetLabel());
7dd16671 969 if (label==0) continue;
970 if (!(esds->At(label))){
971 esds->AddAt(new AliESDtrack(*esdp),label);
972 trds->AddAt(new AliTRDtrack(*trdp),label);
973 } else{
974 AliTRDtrack * track =(AliTRDtrack*)trds->At(label);
975 if (track->GetNumberOfClusters() <trdp->GetNumberOfClusters()){
976 esds->AddAt(new AliESDtrack(*esdp),label);
977 trds->AddAt(new AliTRDtrack(*trdp),label);
978 }
979 }
980 }
981}
982
983
984
9dc17dc7 985void MakeTree()
986{
987 fLoader = AliRunLoader::Open("galice.root");
988 AliESDtrack * esd = new AliESDtrack;
7dd16671 989 AliTRDtrack * ptrd = new AliTRDtrack;
990 AliTRDtrack * ptrd0 = new AliTRDtrack;
991 AliTRDtrack * ptrd1 = new AliTRDtrack;
992 AliTRDtrack * ptrdb = new AliTRDtrack;
993 //
994 AliTRDtrack * trd = ptrd;
995 AliTRDtrack * trd0 = ptrd0;
996 AliTRDtrack * trd1 = ptrd1;
997 AliTRDtrack * trdb = ptrdb;
998 //
9dc17dc7 999 AliMCInfo* info = new AliMCInfo;
1000 AliTRDInfoMC *trdinfo = new AliTRDInfoMC;
7dd16671 1001 AliTRDInfoCl *infocl = new AliTRDInfoCl;
1002 AliTRDInfoSeed *infoSeed = new AliTRDInfoSeed;
1003 // AliTOFtrackInfo *infotof = new AliTOFtrackInfo;
1004 // AliTOFtrackInfo dummytof;
1005 // dummytof.fLab=-1;
1006 AliTRDInfoSeed dummySeed;
1007 AliTRDtrack dummytrd;
1008 AliTRDtrack dummytrd0;
1009 AliTRDtrack dummytrd1;
1010 AliTRDtrack dummytrdb;
1011 AliESDtrack dummyesd;
1012 dummySeed.fLabel =-1;
9dc17dc7 1013 AliTrackReference *ref = new AliTrackReference;
7dd16671 1014 TObjArray * clarray = new TObjArray(6*1000000);
1015 const Int_t kmaxlabel=1500000;
9dc17dc7 1016 //
1017 // MC part
1018 //
1019 TFile f("cmpESDTracks.root");
1020 TTree * tree1 = (TTree*) f.Get("ESDcmpTracks");
1021 TBranch * brmc = tree1->GetBranch("MC");
1022 brmc->SetAddress(&info);
1023 //
1024 //
7dd16671 1025 TObjArray *esds = new TObjArray(kmaxlabel);
1026 TObjArray *trds = new TObjArray(kmaxlabel);
1027 TObjArray *tofs = new TObjArray(kmaxlabel);
1028 TObjArray *seeds = new TObjArray(kmaxlabel);
1029 //
9dc17dc7 1030 //
9dc17dc7 1031 //
9dc17dc7 1032 //
1033 // comparison part
1034 //
1035 TFile f3("trdComp1.root","new");
1036 TTree * tree3 = new TTree("Comparison","Comparison");
1037 tree3->Branch("MC.","AliMCInfo",&info);
1038 tree3->Branch("RC.","AliESDtrack",&esd);
1039 tree3->Branch("TRD.","AliTRDtrack",&trd);
1040 tree3->Branch("TRD0.","AliTRDtrack",&trd0); // track at tof ref
1041 tree3->Branch("TRD1.","AliTRDtrack",&trd1); // track at last TRD ref
1042 tree3->Branch("TRDb.","AliTRDtrack",&trdb); // track at tof ref
1043 tree3->Branch("TR.","AliTrackReference",&ref); // track at last TRD ref
1044 tree3->Branch("InfoMC.","AliTRDInfoMC",&trdinfo); // info about TRD track MC
1045 tree3->Branch("InfoCl.","AliTRDInfoCl",&infocl); // info about TRD track MC
7dd16671 1046 tree3->Branch("InfoSeed.","AliTRDInfoSeed",&infoSeed); // info about TRD track MC
1047 // tree3->Branch("InfoTOF.","AliTOFtrackInfo",&infotof); // info about TOF
9dc17dc7 1048 //
1049 //
1050 Int_t lastevent=-1;
1051
1052 for (Int_t i=0;i<tree1->GetEntries();i++){
1053 brmc->GetEntry(i);
7dd16671 1054 // if (info->fEventNr>5) break;
9dc17dc7 1055 if (lastevent!=info->fEventNr){
1056 printf("Read Event\t%d\n", info->fEventNr);
7dd16671 1057 ReadSeeds(info->fEventNr, seeds);
1058 ReadTOFs(info->fEventNr, tofs);
1059 ReadTracks(info->fEventNr,esds,trds);
9dc17dc7 1060 lastevent = info->fEventNr;
9dc17dc7 1061 clarray->Clear();
1062 LoadClusters(clarray,lastevent);
1063 //
9dc17dc7 1064 }
9dc17dc7 1065 // if (info->fNTOFRef==0 && info->fNTRDRef<5) continue;
1066 Int_t label = info->fLabel;
1067 trdinfo->Update(info);
1068 if (clarray->At(label)){
1069 infocl = (AliTRDInfoCl*)clarray->At(label);
7dd16671 1070 infocl->Update();
1071 }
1072 if (seeds->At(label)){
1073 infoSeed = (AliTRDInfoSeed*)seeds->At(label);
1074 }else{
1075 infoSeed = &dummySeed;
9dc17dc7 1076 }
7dd16671 1077 // if (info->fNTPCRef<3) continue;
9dc17dc7 1078 trd0->SetStop(kTRUE);
1079 trdb->SetStop(kTRUE);
1080 trd1->SetStop(kTRUE);
7dd16671 1081 // infotof = (AliTOFtrackInfo*)(tofs->At(label));
1082// if (!infotof) {
1083// infotof = &dummytof;
1084// infotof->fNCluster=0;
1085// esd = (AliESDtrack*)(esds->At(label));
1086// if (esd){
1087// Double_t pid[5];
1088// esd->GetTPCpid(pid);
1089// Double_t suml=0;
1090// for (Int_t ip=0;ip<5;ip++) suml+=pid[ip];
1091// for (Int_t ip=0;ip<5;ip++) {
1092// if (suml>0.0000000001) infotof->fRefPID[ip]= pid[ip]/suml;
1093// else infotof->fRefPID[ip]=0.2;
1094// }
1095// }
1096// }
1097 trd = &dummytrd;
1098 trd0 = &dummytrd0;
1099 trd1 = &dummytrd1;
1100 trdb = &dummytrdb;
1101 esd = &dummyesd;
9dc17dc7 1102 //
1103 if (esds->At(label)) {
7dd16671 1104 trd0 = ptrd0;
1105 trd1 = ptrd1;
1106 trdb = ptrdb;
9dc17dc7 1107 // if (i%1000==0) printf("\n%d\n",i);
1108 trd = (AliTRDtrack*)(trds->At(label));
1109 esd = (AliESDtrack*)(esds->At(label));
1110 //
1111 *trd0 = *trd;
1112 if (trd->GetBackupTrack())
1113 *trdb = *trd->GetBackupTrack();
1114 else
1115 *trdb = *trd;
1116 *trd1 = *trd;
1117 if (info->fNTOFRef>0){
1118 AliTrackReference * ref = ((AliTrackReference*)(info->fTOFReferences->At(0)));
1119 Double_t x = ((AliTrackReference*)(info->fTOFReferences->At(0)))->X();
1120 Double_t y = ((AliTrackReference*)(info->fTOFReferences->At(0)))->Y();
1121 Double_t phi = TMath::ATan2(y,x);
1122 Double_t radius = TMath::Sqrt(x*x+y*y);
1123 //
1124 // rotate and propagate to local system
1125 //
1126 trd0->SetStop(kFALSE);
1127 if (!trd0->Rotate(phi-trd0->GetAlpha())) trd0->SetStop(kTRUE);
1128 else{
1129 Double_t xyz0[3], param[7];
1130 trd0->GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
1131 Double_t xyz1[3]={ref->X(),ref->Y(),ref->Z()};
70ae41ca 1132 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1133 if (!trd0->PropagateTo(radius,param[1],param[0]*param[4]))
1134 trd0->SetStop(kTRUE);
9dc17dc7 1135 }
1136 trdb->SetStop(kFALSE);
1137 if (!trdb->Rotate(phi-trdb->GetAlpha())) trdb->SetStop(kTRUE);
1138 else
1139 if (!trdb->PropagateTo(radius)) trdb->SetStop(kTRUE);
9dc17dc7 1140 }
1141
1142 if (info->fNTRDRef>0){
1143 Double_t cradius=0;
1144 Double_t cphi=0;
1145 for (Int_t itrd =0; itrd<info->fNTRDRef;itrd++){
1146 Double_t x = ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->X();
1147 Double_t y = ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Y();
1148 Double_t phi = TMath::ATan2(y,x);
1149 // Double_t snp = TMath::ATan2(((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Py(),
1150 // ((AliTrackReference*)(info->fTRDReferences->At(itrd)))->Px());
1151 // if (TMath::Abs(snp)>0.9 &&cradius>10) break;
1152 Double_t radius = TMath::Sqrt(x*x+y*y);
1153 if (radius>cradius){
1154 cradius = radius;
1155 cphi = phi;
1156 ref = ((AliTrackReference*)(info->fTRDReferences->At(itrd)));
1157 }else break;
1158 }
1159 //
1160 // rotate and propagate to local system of the last track reference
1161 //
1162 trd1->SetStop(kFALSE);
1163 if (!trd1->Rotate(cphi-trd1->GetAlpha())) trd1->SetStop(kTRUE);
1164 else{
1165 if (!trd1->PropagateTo(cradius)) trd1->SetStop(kTRUE);
1166 }
1167 }
9dc17dc7 1168 }
7dd16671 1169 tree3->Fill();
9dc17dc7 1170 }
1171 f3.cd();
1172 tree3->Write();
1173}
1174
1175