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