New interface for magnetic field setting (M.Ivanov)
[u/mrichter/AliRoot.git] / TRD / AliTRDComparison.C
1 /*
2
3 //load mag field
4 AliRunLoader *loader = AliRunLoader::Open("galice.root");
5 loader->LoadgAlice();
6 gAlice = loader->GetAliRun();
7 AliMagF * field = gAlice->Field();
8 AliTracker::SetFieldMap(gAlice->Field(),1);
9 TFile fg("geometry.root");
10 gGeoManager = (TGeoManager*)fg.Get("Geo");
11 if (!gGeoManager) TGeoManager::Import("geometry.root");
12 .L AliGenInfo.C+
13 .L AliESDComparisonMI.C+
14 .L AliTRDComparison.C+g
15 //.L AliESDComparisonPIC.C+
16 .x TDR_style.C
17
18 MakeCompTr();
19 MakeTree();
20
21  
22  
23 MakeComp();
24 MakeAlias();
25 TFile fff("TRDdebug.root");
26 TTree * tree =(TTree*)fff.Get("tracklet");
27 TTree * tree2 =(TTree*)fff.Get("Tracks");
28 TTree * treefind =(TTree*)fff.Get("Find");
29 TTree * treebug =(TTree*)fff.Get("Bug");
30 TTree * treetofs =(TTree*)fff.Get("tofseed");
31 AliComparisonDraw compfind;
32 compfind.fTree = treefind;
33 AliComparisonDraw comptr;
34 comptr.fTree = tree;
35
36 TFile fff2("TOFdebug.root");
37 TTree * treetof =(TTree*)fff2.Get("Info");
38 TTree * treetoftrd =(TTree*)fff2.Get("TRD");
39
40 AliComparisonDraw comp2;
41 comp2->fTree =tree;
42 AliComparisonDraw comptof;
43 comptof->fTree =treetof;
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"
70 #include "TGraph.h"
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"
93 #include "AliTRDtracker.h"
94 #include "AliTOFtrack.h"
95 #include "../TOF/AliTOFtrackerMI.h"
96 // COMPARISON includes
97 #include "../STEER/AliGenInfo.h"
98 #include "../STEER/AliESDComparisonMI.h"
99
100
101
102 AliComparisonDraw comp;
103 AliComparisonDraw comptr;
104 AliTRDparameter   par("pica","vyjebana");
105 AliTRDgeometryFull    geom;
106 AliRunLoader  *fLoader=0; //= AliRunLoader::Open(fnGalice);
107
108 TCut cnesd("cnesd","abs(RC.fLabel)!=MC.fLabel");
109 TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
110 TCut citsin("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<5");
111 TCut clength("clength","abs(TRD0.fIntegratedLength-MC.fTOFReferences[0].fLength)<30");
112 TCut c1("c1","TRD1.fStopped==0");    // track propagated to refernce
113 TCut c0("c0","TRD0.fStopped==0");    // track propagated to reference
114 TCut c2("c2","MC.fNTOFRef>0");       //track in tof
115 TCut c3("c3","MC.fNTOFRef>0&&MC.fNTRDRef>5");       //track in tof and TRD
116 TH1F hpull("hpull","hpull",50,-10,10);
117 TH1F hpull6("hpull6","hpull6",50,-6,6);
118
119
120 class 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
129   AliTrackReference fRefLocal[6];    // trd local reference
130   ClassDef(AliTRDInfoMC,1)   // TRD MC information 
131 };
132
133
134 class AliTRDInfoCl: public TObject{
135   public:
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 
147 };
148
149 ClassImp(AliTRDInfoCl)
150
151
152
153 class 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
168 ClassImp(AliTRDInfoMC)
169 ClassImp(AliTRDInfoCl)
170 ClassImp(AliTRDInfoSeed)
171
172 void 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
183 void 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 }
193 void  AliTRDInfoMC::Update(AliMCInfo *info)
194 {
195   // 
196   //
197   Reset();
198   for (Int_t i=0;i<info->fNTRDRef;i++){
199     AliTrackReference * ref = (AliTrackReference*)info->fTRDReferences->At(i);
200     Float_t localX = ref->LocalX();
201     for (Int_t jplane = 0; jplane<6; jplane++){
202       if (TMath::Abs(localX - (300+jplane*14.))<6){
203         fLayer[jplane]++;
204         for (Int_t kplane=jplane;kplane>=0;kplane--){
205           fLayerA[kplane]++;
206           fRef[jplane]=*ref;
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]);
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
226 AliTRDInfoCl::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
243 void  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
280
281
282
283 void  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
312 void  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
340 void 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
363 void 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
385 void 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
394 void 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
402 void 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
426 void 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
441 void 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
448 void 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
463 void 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
506 void 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
544
545
546 void PtRes(){
547   
548   
549 }
550
551
552 TH1F * MakeCumul(TH1F *his,Bool_t norm = kTRUE)
553 {
554   TH1F *hcumul = (TH1F*)his->Clone();
555   char name[1000];
556   sprintf(name,"N%f",gRandom->Rndm());
557   hcumul->SetName(name);
558   hcumul->SetTitle(name);
559   //
560   Float_t all = (Float_t)his->GetSum();
561   if (!norm) all =1;
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
569 TGraph* His01(TH1F *his0,TH1F *his1){
570   Int_t nbins = his0->GetNbinsX();
571   Double_t *x = new Double_t[nbins];
572   Double_t *y = new Double_t[nbins];
573   for (Int_t i=0;i<nbins;i++) {
574     x[i] = his0->GetBinContent(i+1);
575     y[i] = his1->GetBinContent(i+1);
576   }
577   return new TGraph(nbins,x,y);
578 }
579
580
581 TH1F * PullCumul(const char *var, TCut cut, Int_t div=200, Float_t max=20){
582   TH1F * hpullb = new TH1F("hpullb","hpullb",div,0,max);
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
593 void 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
602 void 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");
623   comp.fTree->SetAlias("pullth0","(TRD0.fT[3]-tofth)/sqrt(TRD0.fCtt)");
624   comp.fTree->SetAlias("phi_tof","atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx)");
625   comp.fTree->SetAlias("phi_trd","asin(TRD0.GetSnp())+TRD0.fAlpha");
626   comp.fTree->SetAlias("pullphi0","(sin(asin(TRD0.GetSnp())+TRD0.fAlpha)-sin(atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx)))/sqrt(TRD0.fCee)");
627   comp.fTree->SetAlias("dphir0","(asin(TRD0.GetSnp())+TRD0.fAlpha-atan2(MC.fTOFReferences[0].fPy,MC.fTOFReferences[0].fPx))");
628   comp.fTree->SetAlias("dz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)");
629   comp.fTree->SetAlias("pully0","TRD0.fY/sqrt(TRD0.fCyy)");
630   comp.fTree->SetAlias("pullz0","(TRD0.fZ-MC.fTOFReferences[0].fZ)/sqrt(TRD0.fCzz)");
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   //
645   comp.fTree->SetAlias("TPCE","sqrt(MC.fTPCReferences[2].P()**2+MC.fMass**2)");
646   comp.fTree->SetAlias("TOFE","sqrt(MC.fTOFReferences[0].P()**2+MC.fMass**2)");
647   comp.fTree->SetAlias("TRDE","sqrt(TR.P()**2+MC.fMass**2)");
648   comp.fTree->SetAlias("dtpi","TRD0.fIntegratedTime[2]-10^12*MC.fTOFReferences.fTime");
649   comp.fTree->SetAlias("dtk","TRD0.fIntegratedTime[3]-10^12*MC.fTOFReferences.fTime");
650   comp.fTree->SetAlias("dtp","TRD0.fIntegratedTime[4]-10^12*MC.fTOFReferences.fTime");
651   comp.fTree->SetAlias("dt","(abs(MC.fPdg)==2212)*dtp+(abs(MC.fPdg)==211)*dtpi+(abs(MC.fPdg)==321)*dtk"); //delta time
652   comp.fTree->SetAlias("Beta","sqrt(TR.P()**2/(TR.P()**2+MC.fMass**2))");
653   comp.fTree->SetAlias("dtsigma","10+(1-Beta)*280");
654   comp.fTree->SetAlias("dptrel","(abs(TRD0.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()))");
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()");
670
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()");
673
674 }
675
676
677 void 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
688 void 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
694
695
696 void DrawYResol1(TCut cut){
697   comp.fTree->Draw("TRD1.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c1+cut,"");
698 }
699
700 void DrawYResol0(TCut cut){
701   comp.fTree->Draw("TRD0.fY","MC.fNTOFRef>0&&MC.fNTRDRef>5"+c0+cut,"");
702 }
703
704
705 void DrawTOFResY(Int_t ndiv =10,Float_t ptmin=0.5, Float_t ptmax =1.5, Float_t dymin=-0.5, Float_t dymax=0.5){
706   //
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;
710   //  Int_t   ndiv =5;
711   TCut cl("cl","InfoCl.fNClusters>60");
712   comp.DrawXY("TR.Pt()","TRD0.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
713   TH1F * histof = (TH1F*)comp.fRes->Clone();
714   comp.DrawXY("TR.Pt()","TRD1.fY",c0+c2+cl,"1",ndiv,ptmin, ptmax,dymin, dymax,60);
715   TH1F * histrd = (TH1F*)comp.fRes->Clone();
716   histof->SetMarkerStyle(24);
717   histrd->SetMarkerStyle(25);
718   histof->SetXTitle("P_{t}(GeV)");
719   histof->SetYTitle("#Delta_{r-#phi}(cm)");
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
729 void DrawTOFResZ(Float_t ptmin=0.5, Float_t ptmax =1.5, Float_t dzmin=-2.5, Float_t dzmax=2.5){
730   //
731   comp.DrawXY("TR.Pt()","TRD0.fZ-MC.fTOFReferences.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
732   TH1F * histof = (TH1F*)comp.fRes->Clone();
733   comp.DrawXY("TR.Pt()","TRD1.fZ-TR.fZ",c0+c2,"1",10,ptmin,ptmax,dzmin,dzmax,30);
734   TH1F * histrd = (TH1F*)comp.fRes->Clone();
735   histof->SetMarkerStyle(24);
736   histrd->SetMarkerStyle(25);
737   histof->SetXTitle("P_{t}(GeV)");
738   histof->SetYTitle("#Delta_{z}(cm)");
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
749 void MakeCumul(TCut & cut)
750 {
751   //  TCut cut ="1&&TR.Pt()>0.4";
752   TH1F * hisy = PullCumul("abs(pully0)",cut+cprim+c0+c2+c3,1000,20);
753   TH1F * hisz = PullCumul("abs(pullz0)",cut+cprim+c0+c2+c3,1000,20);
754   TH1F * hisyz = PullCumul("sqrt(pullz0**2+pully0**2)",cut+cprim+c0+c2+c3,1000,20);
755   hisz->SetXTitle("Pull (unit)");
756   hisz->SetYTitle("Cumulative function");
757   hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
758   hisz->Draw(); hisy->Draw("same");
759   //hisyz->Draw("same");  
760   TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD pull (unit)");
761   legend->SetBorderSize(1);
762   legend->AddEntry(hisy ,"Pull y");
763   legend->AddEntry(hisz,"Pull z");
764   //  legend->AddEntry(hisyz,"Pull yz");
765   legend->Draw();
766 }
767
768 void MakeCumulAbs(TCut & cut)
769 {
770   // TCut cut = "1";
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   //
777   hisy->SetLineStyle(1);hisz->SetLineStyle(2);hisyz->SetLineStyle(3);
778   hisz->Draw(); hisy->Draw("same");hisyz->Draw("same");  
779   TLegend *legend = new TLegend(0.55,0.12,0.85,0.35, "TRD residulas (cm)");
780   legend->SetBorderSize(1);
781   legend->AddEntry(hisy ,"Delta y (cm)");
782   legend->AddEntry(hisz,"Delta z (cm)");
783   legend->AddEntry(hisyz,"Delta  (cm)");
784   legend->Draw();
785 }
786
787 void 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;
809       Int_t detector=cl->GetDetector();
810       Int_t localTimeBin=cl->GetLocalTimeBin();
811       //      Int_t sector=geom.GetSector(detector);
812       Int_t plane=geom.GetPlane(detector);
813       for (Int_t ilab=0; ilab<3; ilab++){
814         Int_t label = cl->GetLabel(ilab);
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       }
836     }
837   }  
838   for (Int_t i=0;i<arr->GetEntriesFast();i++){
839     AliTRDInfoCl * info =( AliTRDInfoCl*)arr->At(i); 
840     if (info) info->Update();
841   }
842 }
843
844
845
846 void 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
908 void 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
938 void 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     //
969     Int_t label = TMath::Abs(esdp->GetLabel());
970     if (label==0) continue;
971     if (!(esds->At(label))){
972       esds->AddAt(new AliESDtrack(*esdp),label);
973       trds->AddAt(new AliTRDtrack(*trdp),label);
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
986 void MakeTree()
987 {
988   fLoader = AliRunLoader::Open("galice.root");
989   AliESDtrack * esd     = new AliESDtrack;
990   AliTRDtrack * ptrd     = new AliTRDtrack;
991   AliTRDtrack * ptrd0    = new AliTRDtrack;
992   AliTRDtrack * ptrd1    = new AliTRDtrack;
993   AliTRDtrack * ptrdb    = new AliTRDtrack;
994   //
995   AliTRDtrack * trd     = ptrd;
996   AliTRDtrack * trd0    = ptrd0;
997   AliTRDtrack * trd1    = ptrd1;
998   AliTRDtrack * trdb    = ptrdb;
999   //
1000   AliMCInfo* info       = new AliMCInfo;
1001   AliTRDInfoMC *trdinfo = new AliTRDInfoMC;
1002   AliTRDInfoCl *infocl  = new AliTRDInfoCl;
1003   AliTRDInfoSeed *infoSeed = new AliTRDInfoSeed;
1004   //  AliTOFtrackInfo *infotof = new AliTOFtrackInfo;
1005   // AliTOFtrackInfo dummytof;
1006   //  dummytof.fLab=-1;
1007   AliTRDInfoSeed  dummySeed;
1008   AliTRDtrack dummytrd;
1009   AliTRDtrack dummytrd0;
1010   AliTRDtrack dummytrd1;
1011   AliTRDtrack dummytrdb;
1012   AliESDtrack dummyesd;
1013   dummySeed.fLabel =-1;
1014   AliTrackReference *ref = new AliTrackReference;
1015   TObjArray         * clarray  = new TObjArray(6*1000000);
1016   const Int_t kmaxlabel=1500000;
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   // 
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   //
1031   //
1032   //
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
1047   tree3->Branch("InfoSeed.","AliTRDInfoSeed",&infoSeed);  // info about TRD track MC
1048   //  tree3->Branch("InfoTOF.","AliTOFtrackInfo",&infotof);  // info about TOF
1049   //
1050   //
1051   Int_t lastevent=-1;
1052
1053   for (Int_t i=0;i<tree1->GetEntries();i++){
1054     brmc->GetEntry(i);
1055     //    if (info->fEventNr>5) break;
1056     if (lastevent!=info->fEventNr){
1057       printf("Read Event\t%d\n", info->fEventNr);
1058       ReadSeeds(info->fEventNr, seeds);
1059       ReadTOFs(info->fEventNr, tofs);
1060       ReadTracks(info->fEventNr,esds,trds);
1061       lastevent = info->fEventNr;
1062       clarray->Clear();
1063       LoadClusters(clarray,lastevent);
1064       //
1065     }
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);
1071       infocl->Update();
1072     }
1073     if (seeds->At(label)){
1074       infoSeed = (AliTRDInfoSeed*)seeds->At(label);
1075     }else{
1076       infoSeed = &dummySeed;
1077     }
1078     //    if (info->fNTPCRef<3) continue;
1079     trd0->SetStop(kTRUE);
1080     trdb->SetStop(kTRUE);
1081     trd1->SetStop(kTRUE);
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;
1103     //
1104     if (esds->At(label)) {
1105       trd0  = ptrd0; 
1106       trd1  = ptrd1; 
1107       trdb  = ptrdb; 
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);
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       }    
1168     }
1169     tree3->Fill();    
1170   }
1171   f3.cd();
1172   tree3->Write();
1173 }
1174
1175