]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDComparison.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add pad gain correction + the x position 0 <-> time bin 0 - x increases TOWARDS inter...
[u/mrichter/AliRoot.git] / TRD / 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"
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"
7dd16671 93#include "AliTRDtracker.h"
94#include "AliTOFtrack.h"
95#include "../TOF/AliTOFtrackerMI.h"
9dc17dc7 96// COMPARISON includes
97#include "../STEER/AliGenInfo.h"
98#include "../STEER/AliESDComparisonMI.h"
99
100
101
102AliComparisonDraw comp;
103AliComparisonDraw comptr;
104AliTRDparameter par("pica","vyjebana");
105AliTRDgeometryFull geom;
106AliRunLoader *fLoader=0; //= AliRunLoader::Open(fnGalice);
107
7dd16671 108TCut cnesd("cnesd","abs(RC.fLabel)!=MC.fLabel");
9dc17dc7 109TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
110TCut citsin("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<5");
111TCut clength("clength","abs(TRD0.fIntegratedLength-MC.fTOFReferences[0].fLength)<30");
112TCut c1("c1","TRD1.fStopped==0"); // track propagated to refernce
113TCut c0("c0","TRD0.fStopped==0"); // track propagated to reference
114TCut c2("c2","MC.fNTOFRef>0"); //track in tof
115TCut c3("c3","MC.fNTOFRef>0&&MC.fNTRDRef>5"); //track in tof and TRD
116TH1F hpull("hpull","hpull",50,-10,10);
117TH1F hpull6("hpull6","hpull6",50,-6,6);
118
119
120class AliTRDInfoMC: public TObject{
121 public:
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
7dd16671 129 AliTrackReference fRefLocal[6]; // trd local reference
9dc17dc7 130 ClassDef(AliTRDInfoMC,1) // TRD MC information
131};
132
9dc17dc7 133
134class AliTRDInfoCl: public TObject{
135 public:
7dd16671 136 AliTRDInfoCl();
137 void Update();
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
9dc17dc7 147};
148
149ClassImp(AliTRDInfoCl)
150
151
152
7dd16671 153class AliTRDInfoSeed: public TObject{
154 public:
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
162 Int_t fNLayers; //
163 Int_t fNUsed; //
164 ClassDef(AliTRDInfoSeed,1) // TRD MC information
165};
166
167
168ClassImp(AliTRDInfoMC)
169ClassImp(AliTRDInfoCl)
170ClassImp(AliTRDInfoSeed)
171
172void AliTRDInfoSeed::Update(AliTRDseed *seeds){
173 //
174 //
175 //
176 fNLayers=0;
177 for (Int_t i=0;i<6;i++){
178 fSeeds[i] = seeds[i];
179 if (seeds[i].fN2>10) fNLayers++;
180 }
181}
182
9dc17dc7 183void AliTRDInfoMC::Reset(){
184 //
185 //
186 fNLayers=0;
187 fLast =0;
188 for (Int_t ilayer = 0; ilayer<6; ilayer++){
189 fLayer[ilayer]=0;
190 fLayerA[ilayer]=0;
191 }
192}
193void AliTRDInfoMC::Update(AliMCInfo *info)
194{
195 //
196 //
197 Reset();
9dc17dc7 198 for (Int_t i=0;i<info->fNTRDRef;i++){
199 AliTrackReference * ref = (AliTrackReference*)info->fTRDReferences->At(i);
7dd16671 200 Float_t localX = ref->LocalX();
9dc17dc7 201 for (Int_t jplane = 0; jplane<6; jplane++){
7dd16671 202 if (TMath::Abs(localX - (300+jplane*14.))<6){
9dc17dc7 203 fLayer[jplane]++;
204 for (Int_t kplane=jplane;kplane>=0;kplane--){
205 fLayerA[kplane]++;
206 fRef[jplane]=*ref;
7dd16671 207 //
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()};
212 Float_t dirR[3];
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]);
9dc17dc7 216 }
217 }
218 }
219 }
220 for (Int_t jplane = 0; jplane<6; jplane++){
221 if (fLayer[jplane]>0) fNLayers++;
222 if (fLayer[jplane]>0) fLast = jplane;
223 }
224}
225
7dd16671 226AliTRDInfoCl::AliTRDInfoCl(){
227 //
228 //
229 fNClusters =0;
230 for (Int_t ilayer=0;ilayer<6;ilayer++){
231 fMeanQ[ilayer]=0;
232 fMaxBin[ilayer]=0;
233 fNcl[ilayer]=0;
234 fMeanQS[ilayer]=0;
235 fMaxBinS[ilayer]=0;
236 for (Int_t itime=0;itime<25;itime++){
237 fCount[ilayer][itime]=0;
238 }
239 }
240}
241
242
243void AliTRDInfoCl::Update(){
244 fNClusters = 0;
245 for (Int_t iplane=0;iplane<6;iplane++){
246 fMeanQ[iplane]=-1;
247 fNcl[iplane] = 0;
248 fMaxBin[iplane]=-1;
249 //
250 fMeanQS[iplane]=-1;
251 fMaxBinS[iplane]=-1;
252
253 Float_t maxQ =-10;
254 Float_t maxQS =-10;
255 Int_t sums =0;
256 for (Int_t itime=0;itime<25;itime++){
257 if (fCl[iplane][itime].GetQ()>2){
258 fNcl[iplane]++;
259 fMeanQ[iplane]+=fCl[iplane][itime].GetQ();
260 if (fCl[iplane][itime].GetQ()>=maxQ){
261 fMaxBin[iplane]= itime;
262 maxQ = fCl[iplane][itime].GetQ();
263 }
264 }
265 if (fCl[iplane][itime].GetSumS()>2){
266 sums++;
267 fMeanQS[iplane]+=fCl[iplane][itime].GetSumS();
268 if (fCl[iplane][itime].GetSumS()>=maxQS){
269 fMaxBinS[iplane]= itime;
270 maxQS = fCl[iplane][itime].GetSumS();
271 }
272 }
273 }
274 if (fNcl[iplane]>0) fMeanQ[iplane]/=Float_t(fNcl[iplane]);
275 if (sums>0) fMeanQS[iplane]/=Float_t(sums);
276 fNClusters+=fNcl[iplane];
277 }
278}
279
9dc17dc7 280
281
282
283void MaxBinNorm(const char * expr = "TRD0.fTimBinPlane"){
284 //
285 // max time bin probability distribution
286 //
287 char var[100];
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()));
292 hispi->Draw();
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");
301 hisel->Draw();
302 hispi->Draw("same");
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");
307 legend->Draw();
308
309}
310
311
312void dedxNorm(const char * expr = "TRD0.fdEdxPlane/20"){
313 //
314 // max time bin probability distribution
315 //
316 char var[100];
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()));
321 hispi->Draw();
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");
330 hispi->Draw();
331 hisel->Draw("same");
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");
336 legend->Draw();
337}
338
339
340void EfficiencyPt()
341{
342 //
343 //
344 //
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 (%)");
353 his->Draw();
354 hisf->Draw("same");
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");
359 legend->Draw();
360
361}
362
363void EfficiencyDip()
364{
365 //
366 //
367 //
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 (%)");
376 his->Draw();
377 hisf->Draw("same");
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");
382 legend->Draw();
383}
384
385void ClusterRatio()
386{
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("(%)");
390 hisr->Draw();
391
392}
393
7dd16671 394void DEestimate(){
395 TH1F * hderel = new TH1F("hderel","hderel",100,-150,150);
396 comp.fTree->Draw("100*derel>>hderel",""+c0+c2+cprim);
397 hderel->Fit("gaus");
398 hderel->SetXTitle("Energy loss correction resolution (%)");
399 hderel->Draw();
400}
401
402void DEestimate2(){
403 //
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);
414 hderela->Draw("");
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");
422 legend->Draw();
423}
424
425
426void EffBudget(){
427 //
428 //
429 //
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 (%)");
435 heffb->Fit(&f1);
436 heffb->Draw();
437
438
439}
440
441void TOFbackround(){
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)");
446}
447
448void TOFsignal(){
449 TCut cut = "InfoTOF.fNCluster>0&&abs(MC.fPdg)==2212&&P<2";
450 //
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);
455 his0->Draw();
456 hisf->Draw("same");
457 his0->SetXTitle("TOF Time (ns)");
458}
459
460
461
462
463void EffPIDK(Float_t xmin=0.4, Float_t xmax=2){
464 //
465 //
466 //Float_t xmin=0.4; Float_t xmax=2;
467 // TCut cut= "1";
468// // TCut cut = c2;
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
479// //
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);
489// //
490// his0->Draw();
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");
502// legend->Draw();
503 //
504}
505
506void EffPIDP(Float_t xmin=0.4, Float_t xmax=2){
507 //
508 //
509 //Float_t xmin=0.4; Float_t xmax=4;
510 TCut cut= "1";
511 // TCut cut= c2;
512 //TCut cut = "tofm0"
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);
527// his0->Draw();
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");
537// legend->Draw();
538
539 //
540}
541
542
543
9dc17dc7 544
545
546void PtRes(){
547
548
549}
550
551
7dd16671 552TH1F * MakeCumul(TH1F *his,Bool_t norm = kTRUE)
9dc17dc7 553{
554 TH1F *hcumul = (TH1F*)his->Clone();
7dd16671 555 char name[1000];
556 sprintf(name,"N%f",gRandom->Rndm());
557 hcumul->SetName(name);
558 hcumul->SetTitle(name);
559 //
9dc17dc7 560 Float_t all = (Float_t)his->GetSum();
7dd16671 561 if (!norm) all =1;
9dc17dc7 562 for (Int_t i=0;i<=his->GetNbinsX();i++){
563 hcumul->SetBinContent(i,his->Integral(0,i)/all);
564 }
565 hcumul->SetFillColor(0);
566 return hcumul;
567}
568
7dd16671 569TGraph* 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);
576 }
577 return new TGraph(nbins,x,y);
578}
579
580
9dc17dc7 581TH1F * PullCumul(const char *var, TCut cut, Int_t div=200, Float_t max=20){
582 TH1F * hpullb = new TH1F("hpullb","hpullb",div,0,max);
583 char v2[1000];
584 sprintf(v2,"%s>>hpullb",var);
585 comp.fTree->Draw(v2, cut);
586 TH1F * res = MakeCumul(hpullb);
587 delete hpullb;
588 return res;
589}
590
591
592
593void MakeComp(){
594 //
595 // Set tree and set aliases
596 //
597 TFile *f = new TFile("trdComp1.root");
598 TTree * tree= (TTree*)f->Get("Comparison");
599 comp.fTree =tree;
600}
601
602void MakeAlias(){
603
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()");
614
615
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])");
618
619 comp.fTree->SetAlias("dpt_trd","(abs(1/RC.fRp[4])-TR.Pt())/TR.Pt()");
620
621 comp.fTree->SetAlias("tofth","(MC.fTOFReferences[0].fPz/MC.fTOFReferences[0].Pt())");
622 comp.fTree->SetAlias("dth_tof","TRD0.fT[3]-tofth");
7dd16671 623 comp.fTree->SetAlias("pullth0","(TRD0.fT[3]-tofth)/sqrt(TRD0.fCtt)");
9dc17dc7 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");
7dd16671 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))");
9dc17dc7 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)");
631 //
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)");
635 //
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)");
641 //
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)");
644 //
e2cd3644 645 comp.fTree->SetAlias("TPCE","sqrt(MC.fTPCReferences[2].P()**2+MC.fMass**2)");
9dc17dc7 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.GetPt())-MC.fTOFReferences[0].Pt())/MC.fTOFReferences[0].Pt()");
655 //
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()))");
7dd16671 658 comp.fTree->SetAlias("derel","(TRD0.fDE-(TPCE-TOFE))/(TPCE-TOFE)");
659 //
660 // TOF stuff
661 //
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");
667
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()");
9dc17dc7 670
7dd16671 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()");
9dc17dc7 673
674}
675
676
677void MakeCompTr(){
678 //
679 // Set tree and set aliases
680 //
681 TFile *ff = new TFile("TRDdebug.root");
682 TTree * treetr= (TTree*)ff->Get("tracklet");
683 comptr.fTree =treetr;
684
685}
686
687
688void DrawDE(TCut cut){
689 comp.DrawXY("(TRD0.fDE)","TPCE-TOFE",c2,cut,5,0.01,0.1,0,0.2);
690 comp.fRes->Draw();
691}
692
693
7dd16671 694
695
9dc17dc7 696void DrawYResol1(TCut cut){
697 comp.fTree->Draw("TRD1.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c1+cut,"");
698}
699
700void DrawYResol0(TCut cut){
701 comp.fTree->Draw("TRD0.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c0+cut,"");
702}
703
704
e2cd3644 705void 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 706 //
7dd16671 707 //
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;
e2cd3644 710 // Int_t ndiv =5;
7dd16671 711 TCut cl("cl","InfoCl.fNClusters>60");
712 comp.DrawXY("TR.Pt()","TRD0.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
9dc17dc7 713 TH1F * histof = (TH1F*)comp.fRes->Clone();
7dd16671 714 comp.DrawXY("TR.Pt()","TRD1.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
9dc17dc7 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)");
720 histof->Draw();
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");
726 legend->Draw();
727}
728
7dd16671 729void DrawTOFResZ(Float_t ptmin=0.5, Float_t ptmax =1.5, Float_t dzmin=-2.5, Float_t dzmax=2.5){
9dc17dc7 730 //
7dd16671 731 comp.DrawXY("TR.Pt()","TRD0.fZ-MC.fTOFReferences.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
9dc17dc7 732 TH1F * histof = (TH1F*)comp.fRes->Clone();
7dd16671 733 comp.DrawXY("TR.Pt()","TRD1.fZ-TR.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
9dc17dc7 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)");
739 histof->Draw();
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");
745 legend->Draw();
746}
747
748
749void MakeCumul(TCut & cut)
750{
7dd16671 751 // TCut cut ="1&&TR.Pt()>0.4";
9dc17dc7 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");
7dd16671 757 hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
758 hisz->Draw(); hisy->Draw("same");
759 //hisyz->Draw("same");
9dc17dc7 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");
7dd16671 764 // legend->AddEntry(hisyz,"Pull yz");
9dc17dc7 765 legend->Draw();
766}
767
768void MakeCumulAbs(TCut & cut)
769{
7dd16671 770 // TCut cut = "1";
9dc17dc7 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");
776 //
7dd16671 777 hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
9dc17dc7 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)");
784 legend->Draw();
785}
786
787void LoadClusters(TObjArray *arr, Int_t event)
788{
789 //
790 //
791 // event=0
792 printf("Load cluster - event%d\n",event);
793 fLoader->SetEventNumber(event);
794 AliLoader * trdl = (AliLoader*)fLoader->GetLoader("TRDLoader");
795 trdl->LoadRecPoints();
796
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();
803 //
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);
808 if (!cl) continue;
9dc17dc7 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);
7dd16671 813 for (Int_t ilab=0; ilab<3; ilab++){
814 Int_t label = cl->GetLabel(ilab);
815 Int_t pos = label;
816 if (label<0) continue;
817 if (localTimeBin>21||localTimeBin<0) continue;
818 if (plane<0||plane>5) continue;
819
820 AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(pos);
821 if (!info) {
822 info = new AliTRDInfoCl;
823 arr->AddAt(info,pos);
824 }
825 //
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]++;
830 }
831 }else{
832 info->fCl[plane][localTimeBin]=*cl;
833 info->fCount[plane][localTimeBin]++;
834 }
835 }
9dc17dc7 836 }
837 }
7dd16671 838 for (Int_t i=0;i<arr->GetEntriesFast();i++){
839 AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(i);
840 if (info) info->Update();
841 }
9dc17dc7 842}
843
844
7dd16671 845
846void ReadSeeds(Int_t eventNr, TObjArray *sarray)
847{
848 //
849 // read seeds form the file
850 //
851 AliTRDseed seeds[6];
852 AliTRDseed *pseeds[6]={&seeds[0],&seeds[1],&seeds[2],&seeds[3],&seeds[4],&seeds[5]};
853 Int_t label, label1, label2, seventNr, nused;
854 Float_t fakeRatio;
855 TFile f("TRDdebug.root");
856 TTree *tree =(TTree*)f.Get("Seeds2");
857 if (!tree) return;
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]);
882 //
883 // delete all seeds
884 for (Int_t ip=0;ip<sarray->GetEntriesFast();ip++){
885 AliTRDseed * seed = (AliTRDseed*)sarray->UncheckedAt(ip);
886 if (seed){
887 delete seed;
888 }
889 }
890 sarray->Clear();
891 for (Int_t ip=0;ip<tree->GetEntries();ip++){
892 tree->GetEntry(ip);
893 if (seventNr!=eventNr) continue;
894 //if (nused>15) continue;
895 if (!sarray->At(label)){
896 AliTRDInfoSeed * seed = new AliTRDInfoSeed;
897 seed->Update(seeds);
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);
904 }
905 }
906}
907
908void ReadTOFs(Int_t eventNr, TObjArray *tofs)
909{
910 //
911 // tofs->Clear();
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);
918// //
919// // read tof info
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;
924// //
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);
930// }
931// else{
932// if (infotof->fLength0>100)
933// tofs->AddAt(new AliTOFtrackInfo(*infotof),label);
934// }
935// }
936}
937
938void ReadTracks(Int_t eventNr, TObjArray *esds, TObjArray *trds)
939{
940 //
941 // read debug tracks
942 //
943 AliESDtrack * esdp = new AliESDtrack;
944 AliTRDtrack * trdp = new AliTRDtrack;
945 Int_t teventNr;
946 TFile f2("TRDdebug.root");
947 TTree * tree2 =(TTree*)f2.Get("Tracks");
948 if (!tree2) return;
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);
955 //
956 // delete all event
957 for (Int_t ip=0;ip<esds->GetEntriesFast();ip++){
958 if (esds->At(ip)){
959 delete esds->At(ip);
960 delete trds->At(ip);
961 }
962 }
963 esds->Clear();
964 trds->Clear();
965 for (Int_t ip=0;ip<tree2->GetEntries();ip++){
966 tree2->GetEntry(ip);
967 if (eventNr!=teventNr) continue;
968 //
e2cd3644 969 Int_t label = TMath::Abs(esdp->GetLabel());
7dd16671 970 if (label==0) continue;
971 if (!(esds->At(label))){
972 esds->AddAt(new AliESDtrack(*esdp),label);
973 trds->AddAt(new AliTRDtrack(*trdp),label);
974 } else{
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);
979 }
980 }
981 }
982}
983
984
985
9dc17dc7 986void MakeTree()
987{
988 fLoader = AliRunLoader::Open("galice.root");
989 AliESDtrack * esd = new AliESDtrack;
7dd16671 990 AliTRDtrack * ptrd = new AliTRDtrack;
991 AliTRDtrack * ptrd0 = new AliTRDtrack;
992 AliTRDtrack * ptrd1 = new AliTRDtrack;
993 AliTRDtrack * ptrdb = new AliTRDtrack;
994 //
995 AliTRDtrack * trd = ptrd;
996 AliTRDtrack * trd0 = ptrd0;
997 AliTRDtrack * trd1 = ptrd1;
998 AliTRDtrack * trdb = ptrdb;
999 //
9dc17dc7 1000 AliMCInfo* info = new AliMCInfo;
1001 AliTRDInfoMC *trdinfo = new AliTRDInfoMC;
7dd16671 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;
9dc17dc7 1014 AliTrackReference *ref = new AliTrackReference;
7dd16671 1015 TObjArray * clarray = new TObjArray(6*1000000);
1016 const Int_t kmaxlabel=1500000;
9dc17dc7 1017 //
1018 // MC part
1019 //
1020 TFile f("cmpESDTracks.root");
1021 TTree * tree1 = (TTree*) f.Get("ESDcmpTracks");
1022 TBranch * brmc = tree1->GetBranch("MC");
1023 brmc->SetAddress(&info);
1024 //
1025 //
7dd16671 1026 TObjArray *esds = new TObjArray(kmaxlabel);
1027 TObjArray *trds = new TObjArray(kmaxlabel);
1028 TObjArray *tofs = new TObjArray(kmaxlabel);
1029 TObjArray *seeds = new TObjArray(kmaxlabel);
1030 //
9dc17dc7 1031 //
9dc17dc7 1032 //
9dc17dc7 1033 //
1034 // comparison part
1035 //
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
7dd16671 1047 tree3->Branch("InfoSeed.","AliTRDInfoSeed",&infoSeed); // info about TRD track MC
1048 // tree3->Branch("InfoTOF.","AliTOFtrackInfo",&infotof); // info about TOF
9dc17dc7 1049 //
1050 //
1051 Int_t lastevent=-1;
1052
1053 for (Int_t i=0;i<tree1->GetEntries();i++){
1054 brmc->GetEntry(i);
7dd16671 1055 // if (info->fEventNr>5) break;
9dc17dc7 1056 if (lastevent!=info->fEventNr){
1057 printf("Read Event\t%d\n", info->fEventNr);
7dd16671 1058 ReadSeeds(info->fEventNr, seeds);
1059 ReadTOFs(info->fEventNr, tofs);
1060 ReadTracks(info->fEventNr,esds,trds);
9dc17dc7 1061 lastevent = info->fEventNr;
9dc17dc7 1062 clarray->Clear();
1063 LoadClusters(clarray,lastevent);
1064 //
9dc17dc7 1065 }
9dc17dc7 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);
7dd16671 1071 infocl->Update();
1072 }
1073 if (seeds->At(label)){
1074 infoSeed = (AliTRDInfoSeed*)seeds->At(label);
1075 }else{
1076 infoSeed = &dummySeed;
9dc17dc7 1077 }
7dd16671 1078 // if (info->fNTPCRef<3) continue;
9dc17dc7 1079 trd0->SetStop(kTRUE);
1080 trdb->SetStop(kTRUE);
1081 trd1->SetStop(kTRUE);
7dd16671 1082 // infotof = (AliTOFtrackInfo*)(tofs->At(label));
1083// if (!infotof) {
1084// infotof = &dummytof;
1085// infotof->fNCluster=0;
1086// esd = (AliESDtrack*)(esds->At(label));
1087// if (esd){
1088// Double_t pid[5];
1089// esd->GetTPCpid(pid);
1090// Double_t suml=0;
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;
1095// }
1096// }
1097// }
1098 trd = &dummytrd;
1099 trd0 = &dummytrd0;
1100 trd1 = &dummytrd1;
1101 trdb = &dummytrdb;
1102 esd = &dummyesd;
9dc17dc7 1103 //
1104 if (esds->At(label)) {
7dd16671 1105 trd0 = ptrd0;
1106 trd1 = ptrd1;
1107 trdb = ptrdb;
9dc17dc7 1108 // if (i%1000==0) printf("\n%d\n",i);
1109 trd = (AliTRDtrack*)(trds->At(label));
1110 esd = (AliESDtrack*)(esds->At(label));
1111 //
1112 *trd0 = *trd;
1113 if (trd->GetBackupTrack())
1114 *trdb = *trd->GetBackupTrack();
1115 else
1116 *trdb = *trd;
1117 *trd1 = *trd;
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);
1124 //
1125 // rotate and propagate to local system
1126 //
1127 trd0->SetStop(kFALSE);
1128 if (!trd0->Rotate(phi-trd0->GetAlpha())) trd0->SetStop(kTRUE);
1129 else{
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 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1134 if (!trd0->PropagateTo(radius,param[1],param[0])) trd0->SetStop(kTRUE);
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