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