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