Add to repository
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Jan 2003 15:25:46 +0000 (15:25 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Jan 2003 15:25:46 +0000 (15:25 +0000)
TRD/AliTRDsaveTrackableSeeds.C [new file with mode: 0644]
TRD/AliTRDseedsAnalysis.C [new file with mode: 0644]

diff --git a/TRD/AliTRDsaveTrackableSeeds.C b/TRD/AliTRDsaveTrackableSeeds.C
new file mode 100644 (file)
index 0000000..db5099e
--- /dev/null
@@ -0,0 +1,321 @@
+#ifndef __CINT__
+  #include <iostream.h>
+  #include "AliTRDtracker.h"
+
+  #include "AliTRDcluster.h" 
+  #include "AliTRDhit.h" 
+  #include "AliTRDv1.h"
+  #include "AliTRDgeometry.h"    
+  #include "AliTRDparameter.h"    
+  #include "alles.h"  
+  #include "AliTRDmcTrack.h"
+
+  #include "AliTRDtrack.h"
+  #include "AliTrackReference.h"
+  
+  #include "TFile.h"
+  #include "TParticle.h"
+  #include "TStopwatch.h"
+#endif
+
+
+void AliTRDsaveTrackableSeeds() {
+
+  TObjArray mctracks(2000);
+  TObjArray *TRDmcTracks = &mctracks; 
+
+  TFile *geofile =TFile::Open("AliTRDclusters.root");
+  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
+  Int_t nEvent = 0;
+  Tracker->SetEventNumber(nEvent);
+
+  AliTRDgeometry *fGeo   = (AliTRDgeometry*) geofile->Get("TRDgeometry");
+  //AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter");  
+
+  Char_t *alifile = "galice.root";
+
+  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
+  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
+  if (!gafl) {
+    cout << "Open the ALIROOT-file " << alifile << endl;
+    gafl = new TFile(alifile);
+  }
+  else {
+    cout << alifile << " is already open" << endl;
+  }
+
+  // Get AliRun object from file or create it if not on file
+  gAlice = (AliRun*) gafl->Get("gAlice");
+  if (gAlice)
+    cout << "AliRun object found on file" << endl;
+  else
+    gAlice = new AliRun("gAlice","Alice test program");
+
+  AliTRDv1       *fTRD     = (AliTRDv1*) gAlice->GetDetector("TRD");   
+  
+  // Import the Trees for the event nEvent in the file
+  const Int_t nparticles = gAlice->GetEvent(nEvent);
+
+  printf("found %d particles in event %d \n", nparticles, nEvent);
+    
+  // Create TRDmcTracks for each tpc seed
+  Int_t label = -1, charge = 0;
+  Bool_t primary;
+  Float_t mass;
+  Int_t pdg_code;
+
+  TDirectory *savedir=gDirectory;
+
+  TFile *in=TFile::Open("AliTPCBackTracks.root");
+  if (!in->IsOpen()) {
+    cerr<<"can't open file AliTPCBackTracks.root  !\n"; return;
+  }
+
+  char   tname[100];
+  sprintf(tname,"seedsTPCtoTRD_%d",nEvent);
+  TTree *seedTree=(TTree*)in->Get(tname);
+  if (!seedTree) {
+     cerr<<"AliTRDtracker::PropagateBack(): ";
+     cerr<<"can't get a tree with seeds from TPC !\n";
+  }
+
+  AliTPCtrack *seed=new AliTPCtrack;
+  seedTree->SetBranchAddress("tracks",&seed);
+
+  Int_t nSeeds = (Int_t) seedTree->GetEntries();
+  for (Int_t is=0; is<nSeeds; is++) {
+     seedTree->GetEvent(is);
+     Int_t lbl = seed->GetLabel();
+     if(TMath::Abs(lbl) > nparticles) { 
+       printf("extra high seed label %d \n", lbl);
+       continue;
+     }
+     TParticle *p = gAlice->Particle(TMath::Abs(lbl));
+
+     primary = kTRUE;
+     if (p->GetFirstMother()>=0) primary = kFALSE;
+    
+     pdg_code = (Int_t) p->GetPdgCode();
+
+     if ((pdg_code == 10010020) ||
+        (pdg_code == 10010030) || 
+        (pdg_code == 50000050) || 
+        (pdg_code == 50000051) || 
+        (pdg_code == 10020040)) {
+       mass = 0.;
+       charge = 0;
+     }
+     else {
+       TParticlePDG *pdg = p->GetPDG();
+       charge = (Int_t) pdg->Charge(); 
+       mass=pdg->Mass();
+     }
+
+     AliTRDmcTrack *t = new AliTRDmcTrack(TMath::Abs(lbl), lbl, primary, 
+                                         mass, charge, pdg_code);
+     TRDmcTracks->AddLast(t);
+  }
+  delete seed;
+  delete seedTree;
+
+  savedir->cd();                         
+
+  Int_t i, mctIndex[nparticles];
+  for(i=0; i<nparticles; i++) mctIndex[i]=-1;
+
+  Int_t nMCtracks = TRDmcTracks->GetEntriesFast();
+  AliTRDmcTrack *mct = 0;
+  for(i = 0; i < nMCtracks; i++) {
+    mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(i);
+    mctIndex[mct->GetTrackIndex()] = i;
+  }
+
+  // Loop through the TRD rec points and set XYZ and Pin and Pout
+
+  Double_t xin[6], xout[6];
+  Double_t Px, Py, Pz;
+
+  Float_t  pos[3], rot[3];
+
+  for(Int_t j = 0; j < 6; j++) {
+    xin[j] = Tracker->GetX(0,j,0)-3;
+    xout[j] = Tracker->GetX(0,j,0);
+  }
+
+  TTree *trRefTree = gAlice->TreeTR();
+  if (!trRefTree) {
+    cout << "<AliTRDreadTrackRef> No TR tree found" << endl;
+    return;
+  }       
+
+  Int_t nBytes   = 0;
+  Bool_t fEntrance, fExit;
+
+  Int_t nTrack = (Int_t) trRefTree->GetEntries();
+  cout << "<AliTRDreadTrackRef> Found " << nTrack
+       << " primary particles with track refs" << endl;
+
+  for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+
+    printf("TrackRefs for track %d out of %d \n",iTrack,nTrack);
+    gAlice->ResetTrackReferences();
+    nBytes += trRefTree->GetEvent(iTrack);  
+
+    AliTrackReference *tr = 0;
+    tr = (AliTrackReference*) fTRD->FirstTrackReference(-1);
+    
+    while (tr) {
+
+      Int_t   track = tr->GetTrack();
+      if(mctIndex[track] >= 0) { 
+       mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(mctIndex[track]);
+       
+       pos[0] = tr->X();
+       pos[1] = tr->Y();
+       pos[2] = tr->Z();
+       
+       Double_t phi = TMath::ATan2(pos[1],pos[0]);
+       if(phi < 0) phi = 2 * TMath::Pi() + phi;
+       
+       Int_t sec = ((Int_t) (phi*180/TMath::Pi())) / 20;
+       fGeo->Rotate(fGeo->GetDetector(0,0,17-sec), pos, rot);
+       
+       fEntrance = kFALSE; fExit = kFALSE;
+       for(i = 0; i < 6; i++) {
+         if(TMath::Abs(xin[i] - rot[0]) < 1.4) { fEntrance = kTRUE; break; }
+         if(TMath::Abs(xout[i] - rot[0]) < 1.4) { fExit = kTRUE; break; }
+       }
+      
+       Px = tr->Px();
+       Py = tr->Py();
+       Pz = tr->Pz();
+       
+       if(fEntrance) {
+         mct->SetXYZin(i, rot[0], rot[1], rot[2]);
+         mct->SetPin(i, Px, Py, Pz);
+         /*
+           printf("entr plane %d: rp_x = %f vs %f = xin \n",i,rot[0],xin[i]); 
+           printf("             : Px, Py, Pz = %f, %f, %f \n",Px,Py,Pz); 
+           printf("             : y, z = %f, %f\n",rot[1],rot[2]); 
+         */
+       }
+       if(fExit) {
+         mct->SetXYZout(i, rot[0], rot[1], rot[2]);
+         mct->SetPout(i, Px, Py, Pz);
+         /*
+           printf("exit plane %d: rp_x = %f vs %f = xout \n",i,rot[0],xout[i]); 
+           printf("             : Px, Py, Pz = %f, %f, %f \n",Px,Py,Pz); 
+           printf("             : y, z = %f, %f\n",rot[1],rot[2]);
+         */ 
+       }       
+      }                                   
+      tr = (AliTrackReference *) fTRD->NextTrackReference();
+    }
+  }  
+
+  // Loop through TRD clusters and assign indexes to MC tracks
+  
+  Char_t *clusterfile = "AliTRDclusters.root";
+  TObjArray carray(2000);
+  TObjArray *ClustersArray = &carray;  
+  Tracker->ReadClusters(ClustersArray,clusterfile);
+
+  printf("done with ReadClusters()\n");
+
+  Int_t nClusters = carray.GetEntriesFast();
+
+  Int_t track, det, det0, det1, ltb, plane, ind0, ind1, ncl;
+  Double_t y, y0, y1, q, q0, q1;
+
+  AliTRDcluster *c0 = NULL, *c1 = NULL;
+  printf("nClusters = %d \n", nClusters);
+
+  for (Int_t i = 0; i < nClusters; i++) {
+
+    printf("\r assigning cluster %d out of %d", i, nClusters);
+    AliTRDcluster *cl = (AliTRDcluster *) ClustersArray->UncheckedAt(i);
+
+    for(Int_t j=0; j<3; j++) {
+      track = cl->GetLabel(j);
+      if(track < 0) continue;
+      if(track >= nparticles) { 
+       printf("Track index %d is larger than total number of particles %d\n"
+               ,track,nparticles);
+       continue;
+      }     
+      if(mctIndex[track] < 0) continue;
+      mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(mctIndex[track]);
+
+      label = mct->GetTrackIndex();
+      if(label != track) { 
+       printf("Got label %d while expected %d !\n", label, track);
+       continue;
+      }
+
+      ncl = mct->GetNumberOfClusters();
+      mct->SetNumberOfClusters(ncl+1);
+
+      det=cl->GetDetector();
+      plane = fGeo->GetPlane(det);
+      ltb=cl->GetLocalTimeBin();
+      if((ltb < 0) || (ltb >= kMAX_TB)) continue;
+
+      ind0 = mct->GetClusterIndex(ltb, plane, 0);
+      ind1 = mct->GetClusterIndex(ltb, plane, 1);
+      if(ind0 < 0) {
+       mct->Update(ltb,plane,0,i);            
+      } else {
+       c0 = (AliTRDcluster *) ClustersArray->UncheckedAt(ind0);
+       det0 = c0->GetDetector();
+       y =  cl->GetY();
+       y0 = c0->GetY();
+       q  = TMath::Abs(cl->GetQ());
+       q0 = TMath::Abs(c0->GetQ());
+       if((det == det0) && (TMath::Abs(y-y0) < 1.5)) {
+         if(q > q0) mct->Update(ltb,plane,0,i);
+       } else {
+         if (ind1 < 0) {
+           mct->Update(ltb,plane,1,i);
+         } 
+         else {
+           c1 = (AliTRDcluster *) ClustersArray->UncheckedAt(ind1);
+           det1 = c1->GetDetector();
+           y1 = c1->GetY();
+           q1 = TMath::Abs(c1->GetQ());
+           if((det == det0) && (TMath::Abs(y-y0) < 1.5)) {
+             if(q > q1) mct->Update(ltb,plane,1,i);
+           }     
+         }       
+       }
+      }
+    }
+  }        
+       
+  // Write TRDmcTracks in output file  
+
+  savedir=gDirectory; 
+  Char_t *filename   = "AliTRDtrackableSeeds.root";
+  TFile *out = new TFile(filename,"RECREATE");
+
+  TTree *tracktree = new TTree("MCtracks","TRD MC tracks");
+
+  AliTRDmcTrack *iotrack=0;
+  tracktree->Branch("MCtracks","AliTRDmcTrack",&iotrack,32000,0);
+  
+  Int_t ntracks = TRDmcTracks->GetEntriesFast();
+  
+  for (Int_t i=0; i<ntracks; i++) {
+    AliTRDmcTrack *pt=(AliTRDmcTrack*)TRDmcTracks->UncheckedAt(i);    
+    iotrack=pt;
+    tracktree->Fill();
+    printf("Put track with label %d and %d clusters in the output tree \n",
+          pt->GetTrackIndex(),pt->GetNumberOfClusters());
+  }
+
+  tracktree->Write();
+  out->Close();     
+  savedir->cd(); 
+
+  return;
+}
+
diff --git a/TRD/AliTRDseedsAnalysis.C b/TRD/AliTRDseedsAnalysis.C
new file mode 100644 (file)
index 0000000..682b3f9
--- /dev/null
@@ -0,0 +1,618 @@
+#ifndef __CINT__
+  #include <iostream.h>
+  #include "AliTRDtracker.h"
+  #include "AliTRDcluster.h" 
+  #include "AliTRDv1.h"
+  #include "AliTRDgeometry.h"    
+  #include "AliTRDparameter.h"    
+  #include "alles.h"  
+  #include "AliTRDmcTrack.h"
+  #include "AliTRDtrack.h"
+  
+  #include "TFile.h"
+  #include "TParticle.h"
+  #include "TStopwatch.h"
+#endif
+
+
+void AliTRDseedsAnalysis() {
+
+  const Int_t nPrimaries = (Int_t) 3*86030/4;
+  //  const Int_t nPrimaries = 500;
+
+  Bool_t page[6]; for(Int_t i=0; i<6; i++) page[i]=kTRUE;
+
+  const Int_t nPtSteps = 30;
+  const Float_t maxPt = 3.;
+  const Float_t minPt = 0.;
+  
+  const Int_t nEtaSteps = 40;
+  const Float_t maxEta = 1.;
+  const Float_t minEta = -1.;
+  
+  // page[0]
+  TH1F *hNcl=new TH1F("hNcl","Seeds No of Clusters",255,-9.5,500.5); 
+  hNcl->SetFillColor(4);
+  hNcl->SetXTitle("No of clusters"); 
+  hNcl->SetYTitle("counts"); 
+  TH1F *hNtb=new TH1F("hNtb","Seeds No of TB with clusters",160,-9.5,150.5); 
+  hNtb->SetFillColor(4);
+  hNtb->SetXTitle("No of timebins with clusters"); 
+  hNtb->SetYTitle("counts"); 
+  TH1F *hNep=new TH1F("hNep","Seeds end point",160, -9.5, 150.5); 
+  hNep->SetFillColor(4);
+  hNep->SetXTitle("outermost timebin with cluster"); 
+  hNep->SetYTitle("counts"); 
+  TH1F *hNmg=new TH1F("hNmg","Seeds max gap",160, -9.5, 150.5); 
+  hNmg->SetFillColor(4);
+  hNmg->SetXTitle("max gap (tb)"); 
+  hNmg->SetYTitle("counts"); 
+
+  // page[1]
+  Float_t fMin = -5.5, fMax = 134.5;
+  Int_t   iChan = 140; 
+  TH2F *h2ep = new TH2F("h2ep","MC vs RT (end point)",iChan,fMin,fMax,iChan,fMin,fMax); 
+  h2ep->SetMarkerColor(4);
+  h2ep->SetMarkerSize(2);
+  TH2F *h2ntb = new TH2F("h2ntb","MC vs RT (TB with Clusters)",iChan,fMin,fMax,iChan,fMin,fMax); 
+  h2ntb->SetMarkerColor(4);
+  h2ntb->SetMarkerSize(2);
+  TH2F *h2mg = new TH2F("h2mg","MC vs RT (max gap)",iChan/2,fMin,fMax/2,iChan/2,fMin,fMax/2); 
+  h2mg->SetMarkerColor(4);
+  h2mg->SetMarkerSize(2);
+
+  // page[2]
+  TH1F *hPt_all=new TH1F("hPt_all","Seeds Pt",nPtSteps,minPt,maxPt);
+  hPt_all->SetLineColor(4);
+  hPt_all->SetXTitle("Pt (GeV/c)"); 
+  hPt_all->SetYTitle("counts"); 
+
+  TH1F *hPt_short=new TH1F("hPt_short","Short seeds Pt",nPtSteps,minPt,maxPt);
+  hPt_short->SetLineColor(8);
+  TH1F *hPt_long=new TH1F("hPt_long","Long seeds Pt",nPtSteps,minPt,maxPt);
+  hPt_long->SetLineColor(2);
+  
+  TH1F *hrtPt_short=new TH1F("hrtPt_short","RT from short seeds",nPtSteps,minPt,maxPt);
+  hrtPt_short->SetLineColor(8);
+  TH1F *hrtPt_long=new TH1F("hrtPt_long","RT from long seeds",nPtSteps,minPt,maxPt);
+  hrtPt_long->SetLineColor(2);
+
+  TH1F *hEta_all=new TH1F("hEta_all","Seeds Eta",nEtaSteps,minEta,maxEta);
+  hEta_all->SetLineColor(4);
+  hEta_all->SetXTitle("Eta"); 
+  hEta_all->SetYTitle("counts"); 
+  TH1F *hEta_short=new TH1F("hEta_short","Short seeds Eta",nEtaSteps,minEta,maxEta);
+  hEta_short->SetLineColor(8);
+  TH1F *hEta_long=new TH1F("hEta_long","Long seeds Eta",nEtaSteps,minEta,maxEta);
+  hEta_long->SetLineColor(2);
+  
+  TH1F *hrtEta_short=new TH1F("hrtEta_short","RT from short seeds",nEtaSteps,minEta,maxEta);
+  hrtEta_short->SetLineColor(8);
+  TH1F *hrtEta_long=new TH1F("hrtEta_long","RT from long seeds",nEtaSteps,minEta,maxEta);
+  hrtPt_long->SetLineColor(2);
+
+  // page[3]
+  TH1F *hEff_long = new TH1F("hEff_long","Efficiency vs Pt",nPtSteps,minPt,maxPt); 
+  hEff_long->SetLineColor(2); hEff_long->SetLineWidth(2);  
+  hEff_long->SetXTitle("Pt, GeV/c"); 
+  hEff_long->SetYTitle("Efficiency"); 
+  
+  TH1F *hEff_short = new TH1F("hEff_short","Efficiency short",nPtSteps,minPt,maxPt); 
+  hEff_short->SetLineColor(8); hEff_short->SetLineWidth(2);  
+  hEff_short->SetXTitle("Pt, GeV/c"); 
+  hEff_short->SetYTitle("Efficiency"); 
+  
+  TH1F *hEff_long_eta = new TH1F("hEff_long_eta","Efficiency vs Eta",nEtaSteps,minEta,maxEta); 
+  hEff_long_eta->SetLineColor(2); hEff_long_eta->SetLineWidth(2);  
+  hEff_long_eta->SetXTitle("Pt, GeV/c"); 
+  hEff_long_eta->SetYTitle("Efficiency"); 
+  
+  TH1F *hEff_short_eta = new TH1F("hEff_short_eta","Efficiency short",nEtaSteps,minEta,maxEta); 
+  hEff_short_eta->SetLineColor(8); hEff_short_eta->SetLineWidth(2);  
+  hEff_short_eta->SetXTitle("Pt, GeV/c"); 
+  hEff_short_eta->SetYTitle("Efficiency"); 
+  
+  // page[4]
+  TH1F *hY=new TH1F("hY","Y resolution",50,-0.5,0.5);hY->SetLineColor(4);
+  hY->SetLineWidth(2);
+  hY->SetXTitle("(cm)");
+  TH1F *hZ=new TH1F("hZ","Z resolution",80,-4,4);hZ->SetLineColor(4);
+  hZ->SetLineWidth(2);
+  hZ->SetXTitle("(cm)");
+  TH1F *hPhi=new TH1F("hPhi","PHI resolution",100,-20.,20.); 
+  hPhi->SetFillColor(4);
+  hPhi->SetXTitle("(mrad)");
+  TH1F *hLambda=new TH1F("hLambda","Lambda resolution",100,-100,100);
+  hLambda->SetFillColor(17);
+  hLambda->SetXTitle("(mrad)");
+  TH1F *hPt=new TH1F("hPt","Relative Pt resolution",30,-10.,10.);  
+  
+  // page[5]
+  TH1F *hY_n=new TH1F("hY_n","Normalized Y resolution",50,-0.5,0.5); hY_n->SetLineColor(4);
+  hY_n->SetLineWidth(2);
+  hY_n->SetXTitle("  ");
+  TH1F *hZ_n=new TH1F("hZ_n","Normalized Z resolution",50,-0.5,0.5); hZ_n->SetLineColor(2);
+  hZ_n->SetLineWidth(2);
+  hZ_n->SetXTitle("   ");
+  TH1F *hLambda_n=new TH1F("hLambda_n","Normalized Lambda resolution",50,-0.5,0.5);
+  hLambda_n->SetFillColor(17);
+  TH1F *hPt_n=new TH1F("hPt_n","Normalized Pt resolution",50,-0.5,0.5); 
+  
+
+  Int_t nEvent = 0;
+
+  // Load Seeds saved as MC tracks 
+  TObjArray mctarray(2000);
+  TFile *mctf=TFile::Open("AliTRDtrackableSeeds.root");
+  if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDtrackableSeeds.root !\n"; return;}
+  TTree *mctracktree=(TTree*)mctf->Get("MCtracks");
+  TBranch *mctbranch=mctracktree->GetBranch("MCtracks");
+  Int_t const nMCtracks = (Int_t) mctracktree->GetEntries();
+  cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
+  for (Int_t i=0; i<nMCtracks; i++) {
+    AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack();
+    mctbranch->SetAddress(&ioMCtrack);
+    mctracktree->GetEvent(i);
+    mctarray.AddLast(ioMCtrack);
+  }
+  mctf->Close();                 
+
+  TFile *tf=TFile::Open("AliTRDtracks.root");
+
+  if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
+  TObjArray rtarray(2000);
+
+  char   tname[100];
+  sprintf(tname,"TRDb_%d",nEvent);     
+  TTree *tracktree=(TTree*)tf->Get(tname);
+
+  TBranch *tbranch=tracktree->GetBranch("tracks");
+
+  Int_t nRecTracks = (Int_t) tracktree->GetEntries();
+  cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;
+
+  for (Int_t i=0; i<nRecTracks; i++) {
+    AliTRDtrack *iotrack=new AliTRDtrack();
+    tbranch->SetAddress(&iotrack);
+    tracktree->GetEvent(i);
+    rtarray.AddLast(iotrack);
+  }
+  tf->Close();                 
+
+  //  return;
+
+  AliTRDmcTrack *mct;
+  AliTRDtrack *rt;
+  Int_t rtIndex[nMCtracks];
+  for(Int_t j = 0; j < nMCtracks; j++) {
+    mct = (AliTRDmcTrack*) mctarray.UncheckedAt(j);
+    rtIndex[j] = 0;
+    for (Int_t i=0; i<nRecTracks; i++) {
+      rt = (AliTRDtrack*) rtarray.UncheckedAt(i);
+      Int_t label = rt->GetSeedLabel();
+      if(mct->GetTrackIndex() == label) rtIndex[j] = i;
+    }
+  } 
+
+  // Load clusters
+
+  TFile *geofile =TFile::Open("AliTRDclusters.root");   
+  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
+  Tracker->SetEventNumber(nEvent);
+
+  AliTRDgeometry *fGeom   = (AliTRDgeometry*) geofile->Get("TRDgeometry"); 
+  AliTRDparameter *fPar   = (AliTRDparameter*) geofile->Get("TRDparameter"); 
+  
+  Char_t *alifile = "AliTRDclusters.root";
+  TObjArray carray(2000);
+  TObjArray *ClustersArray = &carray;
+  Tracker->ReadClusters(ClustersArray,alifile);   
+
+
+  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
+  alifile = "galice.root";
+  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
+  if (!gafl) {
+    cout << "Open the ALIROOT-file " << alifile << endl;
+    gafl = new TFile(alifile);
+  }
+  else {
+    cout << alifile << " is already open" << endl;
+  }
+  gAlice = (AliRun*) gafl->Get("gAlice");
+  if (gAlice)
+    cout << "AliRun object found on file" << endl;
+  else
+    gAlice = new AliRun("gAlice","Alice test program");
+
+  AliTRDv1       *TRD = (AliTRDv1*) gAlice->GetDetector("TRD");    
+  const Int_t nparticles = gAlice->GetEvent(nEvent);
+  if (nparticles <= 0) return;
+
+  TParticle *p;
+  Bool_t electrons[300000];
+  for(Int_t i = 0; i < 300000; i++) electrons[i] = kFALSE;
+
+  Bool_t mark_electrons = kFALSE;
+  if(mark_electrons) {
+    printf("mark electrons\n");
+
+    for(Int_t i = nPrimaries; i < nparticles; i++) {
+      p = gAlice->Particle(i);
+      if(p->GetMass() > 0.01) continue;
+      if(p->GetMass() < 0.00001) continue;
+      electrons[i] = kTRUE;
+    }
+  }
+
+  AliTRDcluster *cl = 0;
+  Int_t nw, label, index, ti[3];
+  Int_t mct_tbwc, rt_tbwc, mct_max_gap, rt_max_gap, mct_sector, rt_sector; 
+  Int_t mct_max_tb, rt_max_tb, mct_min_tb, rt_min_tb;
+
+  Int_t det, plane, ltb, gtb, gap, max_gap, sector;
+  Double_t Pt, Px, Py, Pz;
+  Double_t mcPt, mcPx, mcPy, mcPz;
+  Double_t x,y,z, mcX;
+  Int_t rtClusters, rtCorrect;
+
+  Int_t nToFind_long, nFound_long, nToFind_short, nFound_short;
+
+  printf("\n");
+
+  Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
+  Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
+  Double_t dx = (Double_t) fPar->GetTimeBinSize();   
+
+  Int_t tbAmp = fPar->GetTimeBefore();
+  Int_t maxAmp = 0;
+  Int_t tbDrift = fPar->GetTimeMax();
+  Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
+
+  tbDrift = TMath::Min(tbDrift,maxDrift);
+  tbAmp = TMath::Min(tbAmp,maxAmp);
+
+  const Int_t nPlanes = fGeom->Nplan();
+  const Int_t tbpp = Tracker->GetTimeBinsPerPlane();
+  const Int_t nTB = tbpp * nPlanes;
+  Int_t mask[nTB][2];
+
+  nToFind_long = 0; nFound_long = 0;
+  nToFind_short = 0; nFound_short = 0;
+
+  for(Int_t i=0; i < nMCtracks; i++) {
+    mct = (AliTRDmcTrack *) mctarray.UncheckedAt(i);
+    label = mct->GetTrackIndex();
+
+    // Waveform of the MC track
+    for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
+    Int_t mct_ncl = mct->GetNumberOfClusters();
+
+    for(ltb = 0; ltb < kMAX_TB; ltb++) {
+      for(plane = 0; plane < nPlanes; plane++) {
+       for(Int_t n = 0; n < 2; n++) {
+         index = mct->GetClusterIndex(ltb,plane,n);
+         if(index < 0) continue; 
+         cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
+
+         for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
+
+         if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
+           printf("mc wrong track label: %d, %d, %d != %d \n",
+                  ti[0], ti[1], ti[2], label);
+         } 
+         det=cl->GetDetector();
+         if(fGeom->GetPlane(det) != plane) 
+           printf("cluster plane = %d != %d expected plane\n", 
+                  fGeom->GetPlane(det), plane);
+         if(cl->GetLocalTimeBin() != ltb) 
+           printf("cluster ltb = %d != %d expected ltb\n", 
+                  cl->GetLocalTimeBin(), ltb);
+         gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+         mask[gtb][n] = index;
+       }
+      }
+    }
+    
+    for(plane = 0; plane < nPlanes; plane++) {
+      for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
+       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+       if(mask[gtb][0] > -1) break;
+      }
+      if(ltb >= -tbAmp) break;
+    }  
+    if((plane == nPlanes) && (ltb == -tbAmp-1)) {
+      // printf("warning: for track %d min tb is not found and set to %d!\n",
+      //        label, nTB-1);
+      mct_min_tb = nTB-1;
+      // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
+    }
+    else {
+      mct_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+    }
+
+    for(plane = nPlanes-1 ; plane>=0; plane--) {
+      for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
+       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+       if(mask[gtb][0] > -1) break;
+      }
+      if(ltb < tbDrift) break;
+    }  
+    if((plane == -1) && (ltb == tbDrift)) {
+      //      printf("warning: for track %d max tb is not found and set to 0!\n",label);
+      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
+      mct_max_tb = 0;
+      //      mcX = Tracker->GetX(0,0,tbDrift-1);
+    }
+    else {
+      mct_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
+      //      mcX = Tracker->GetX(0,plane,ltb);
+      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
+      det=cl->GetDetector();
+      sector = fGeom->GetSector(det);  
+      mct_sector = sector;
+    }
+
+    mct_tbwc = 0;
+    mct_max_gap = 0;
+    gap = -1;
+    
+    for(nw = mct_min_tb; nw < mct_max_tb+1; nw++) {
+      gap++;
+      if(mask[nw][0] > -1) {
+       mct_tbwc++;
+       if(gap > mct_max_gap) mct_max_gap = gap;
+       gap = 0;
+      }
+    }
+
+    //  Waveform of the reconstructed track
+    if(rtIndex[i] >= 0) rt = (AliTRDtrack *) rtarray.UncheckedAt(rtIndex[i]);
+
+    for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; }
+    Int_t rt_ncl = rt->GetNumberOfClusters();
+
+    for(Int_t n = 0; n < rt_ncl; n++) {
+      index = rt->GetClusterIndex(n);
+      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
+      
+      for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
+      
+      if((ti[0] != label) && (ti[1] != label) &&  (ti[2] != label)) {
+       // printf("rt wrong track label: %d, %d, %d != %d \n", ti[0], ti[1], ti[2], label);
+       continue;
+      } 
+      
+      det=cl->GetDetector();
+      plane = fGeom->GetPlane(det); 
+      ltb = cl->GetLocalTimeBin(); 
+      gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+      mask[gtb][0] = index;
+    }
+
+    for(plane = 0; plane < nPlanes; plane++) {
+      for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
+       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+       if(mask[gtb][0] > -1) break;
+      }
+      if(ltb >= -tbAmp) break;
+    }  
+    if((plane == nPlanes) && (ltb == -tbAmp-1)) {
+      // printf("warning: for track %d min tb is not found and set to %d!\n",
+      //            label, nTB-1);
+      rt_min_tb = nTB-1;
+      // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
+    }
+    else {
+      rt_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+    }
+
+    for(plane = nPlanes-1 ; plane>=0; plane--) {
+      for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
+       gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
+       if(mask[gtb][0] > -1) break;
+      }
+      if(ltb < tbDrift) break;
+    }  
+    if((plane == -1) && (ltb == tbDrift)) {
+      // printf("warning: for track %d max tb is not found and set to 0!\n",label);
+      //      for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
+      rt_max_tb = 0;
+      //      mcX = Tracker->GetX(0,0,tbDrift-1);
+    }
+    else {
+      rt_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
+      //      mcX = Tracker->GetX(0,plane,ltb);
+      cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]);
+      det=cl->GetDetector();
+      sector = fGeom->GetSector(det);  
+      rt_sector = sector;
+    }
+
+    rt_tbwc = 0;
+    rt_max_gap = 0;
+    gap = -1;
+    
+    for(nw = rt_min_tb; nw < rt_max_tb+1; nw++) {
+      gap++;
+      if(mask[nw][0] > -1) {
+       rt_tbwc++;
+       if(gap > rt_max_gap) rt_max_gap = gap;
+       gap = 0;
+      }
+    }
+
+    // Fill the histoes
+
+    if(page[0]) {
+      hNcl->Fill((Float_t) mct_ncl);
+      hNtb->Fill((Float_t) mct_tbwc);
+      hNep->Fill((Float_t) mct_max_tb);
+      hNmg->Fill((Float_t) mct_max_gap);
+    }
+    if(page[1]) {
+      h2ep->Fill((Float_t) rt_max_tb, (Float_t) mct_max_tb);
+      h2ntb->Fill((Float_t) rt_tbwc, (Float_t) mct_tbwc);
+      h2mg->Fill((Float_t) rt_max_gap, (Float_t) mct_max_gap);
+    }
+    if(page[2]) {
+      p = gAlice->Particle(label);
+      hPt_all->Fill(p->Pt());
+      hEta_all->Fill(p->Eta());
+      if(mct_max_tb > 60) {
+       nToFind_long++;
+       hPt_long->Fill(p->Pt());
+       hEta_long->Fill(p->Eta());
+       if(((mct_max_tb - rt_max_tb) < 10) && 
+          (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
+         nFound_long++;
+         hrtPt_long->Fill(p->Pt());
+         hrtEta_long->Fill(p->Eta());
+       }         
+      }
+      if((mct_max_tb < 60) && (mct_max_tb > 10)) {
+       nToFind_short++;
+       hPt_short->Fill(p->Pt());
+       hEta_short->Fill(p->Eta());
+       if(((mct_max_tb - rt_max_tb) < 10) && 
+          (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) {
+         nFound_short++;
+         hrtPt_short->Fill(p->Pt());
+         hrtEta_short->Fill(p->Eta());
+       }         
+      }
+    }
+    if(page[4] && page[5]) {
+      if((mct_tbwc > 50) && (rt_tbwc > 50)) {
+       if(rt->GetSeedLabel() != label) {
+         printf("mct vs rt index mismatch: %d != %d \n",
+                label, rt->GetSeedLabel());
+         return;
+       }
+       Pt = TMath::Abs(rt->GetPt()); 
+       Double_t cc = TMath::Abs(rt->GetSigmaC2()); 
+       mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1);
+       rt->GetPxPyPz(Px,Py,Pz);      
+
+       printf("\n\ntrack %d \n", label);
+       printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz); 
+       printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz); 
+    
+       mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy);
+       if(mcPt < 0.0001) mcPt = 0.0001;
+    
+       Float_t lamg=TMath::ATan2(mcPz,mcPt);
+       Float_t lam =TMath::ATan2(Pz,Pt);
+       if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
+       else hPt_n->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc));
+    
+       Float_t phig=TMath::ATan2(mcPy,mcPx);
+       Float_t phi =TMath::ATan2(Py,Px);    
+       //      if(!(rt->PropagateTo(x))) continue;
+       //    if((mcPt > Pt_min) && (mcPt < Pt_max)) {
+       hLambda->Fill((lam - lamg)*1000.);
+       hLambda_n->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2()));
+       if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
+       else hPt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.); 
+       hPhi->Fill((phi - phig)*1000.);
+       hY->Fill(rt->GetY() - y);
+       hZ->Fill(rt->GetZ() - z);
+       hY_n->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2()));
+       hZ_n->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2()));
+      }
+    }  
+  }
+
+  // plot pictures
+
+  if(page[0]) {
+    TCanvas *c0=new TCanvas("c0","",0,0,700,850);
+    gStyle->SetOptStat(111110);
+    c0->Divide(2,2);
+    c0->cd(1); gPad->SetLogy(); hNcl->Draw();
+    c0->cd(2); hNtb->Draw();
+    c0->cd(3); hNep->Draw();
+    c0->cd(4); hNmg->Draw();
+    c0->Print("c0.ps","ps"); 
+  }
+  if(page[1]) {
+    TCanvas *c1=new TCanvas("c1","",0,0,700,850);
+    gStyle->SetOptStat(0);
+    c1->Divide(2,2);
+    c1->cd(1); h2ep->Draw();
+    c1->cd(2); h2ntb->Draw();
+    c1->cd(3); h2mg->Draw();
+    //    c1->cd(4); hNmg->Draw();
+    c1->Print("c1.ps","ps"); 
+  }
+  if(page[2]) {
+    TCanvas *c2=new TCanvas("c2","",0,0,700,850);
+    gStyle->SetOptStat(0);
+    c2->Divide(2,2);
+    c2->cd(1); hPt_all->Draw(); hPt_long->Draw("same"); hPt_short->Draw("same"); 
+    c2->cd(2); hEta_all->Draw(); hEta_long->Draw("same"); hEta_short->Draw("same"); 
+    //    c2->cd(3); h2mg->Draw();
+    //    c2->cd(4); hNmg->Draw();
+    c2->Print("c2.ps","ps"); 
+  }
+  if(page[3]) {
+    TCanvas *c3=new TCanvas("c3","",0,0,700,850);
+    gStyle->SetOptStat(0);
+    c3->Divide(1,2);
+    c3->cd(1);
+    hrtPt_long->Sumw2(); hPt_long->Sumw2(); 
+    hEff_long->Divide(hrtPt_long,hPt_long,1,1.,"b");
+    hEff_long->SetMaximum(1.4);
+    hEff_long->SetYTitle("Matching Efficiency");
+    hEff_long->SetXTitle("Pt (GeV/c)");
+    hEff_long->Draw();          
+    hrtPt_short->Sumw2(); hPt_short->Sumw2(); 
+    hEff_short->Divide(hrtPt_short,hPt_short,1,1.,"b");
+    hEff_short->SetMaximum(1.4);
+    hEff_short->SetYTitle("Matching Efficiency");
+    hEff_short->SetXTitle("Pt (GeV/c)");
+    hEff_short->Draw("same");          
+
+    c3->cd(2);
+    hrtEta_long->Sumw2(); hEta_long->Sumw2(); 
+    hEff_long_eta->Divide(hrtEta_long,hEta_long,1,1.,"b");
+    hEff_long_eta->SetMaximum(1.4);
+    hEff_long_eta->SetYTitle("Matching Efficiency");
+    hEff_long_eta->SetXTitle("Eta");
+    hEff_long_eta->Draw();          
+    hrtEta_short->Sumw2(); hEta_short->Sumw2(); 
+    hEff_short_eta->Divide(hrtEta_short,hEta_short,1,1.,"b");
+    hEff_short_eta->SetMaximum(1.4);
+    hEff_short_eta->SetYTitle("Matching Efficiency");
+    hEff_short_eta->SetXTitle("Eta");
+    hEff_short_eta->Draw("same");          
+    c3->Print("c3.ps","ps"); 
+  }
+  if(page[4]) {
+    TCanvas *c4=new TCanvas("c4","",0,0,700,850);
+    gStyle->SetOptStat(111110);
+    c4->Divide(2,3);
+    c4->cd(1); hY->Draw();
+    c4->cd(2); hZ->Draw();
+    c4->cd(3); hPhi->Draw();
+    c4->cd(4); hLambda->Draw();
+    c4->cd(5); hPt->Draw();
+    c4->Print("c4.ps","ps"); 
+  }
+  if(page[5]) {
+    TCanvas *c5=new TCanvas("c5","",0,0,700,850);
+    gStyle->SetOptStat(111110);
+    c5->Divide(2,3);
+    c5->cd(1); hY_n->Draw();
+    c5->cd(2); hZ_n->Draw();
+    c5->cd(3); hLambda_n->Draw();
+    c5->cd(4); hPt_n->Draw();
+    c5->Print("c5.ps","ps"); 
+  }
+  printf("Efficiency (long) = %d/%d = %f\n",nFound_long,nToFind_long,
+        ((Float_t) nFound_long / (Float_t) nToFind_long));
+  printf("Efficiency (shrt) = %d/%d = %f\n",nFound_short,nToFind_short,
+        ((Float_t) nFound_short / (Float_t) nToFind_short));
+}
+
+