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