FMD/TestReconstruction.C

   1 void  TestReconstruction (Int_t vol=1, const Int_t nRings=128, const Int_t nSectors=20)
   2 {
   3   // Dynamically link some shared libs
   4   if (gClassTable->GetID("AliRun") < 0)
   5     {
   6       gROOT->LoadMacro("loadlibs.C");
   7       loadlibs();
   8     }
   9   // Connect the Root Galice file containing Geometry, Kine and Hits
  10   char filename[]="galice.root";
  11   TFile *file =  (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
  12   if (!file) file = new TFile(filename);
  13   // Get AliRun object from file or create it if not on file
  14   if (!gAlice)
  15     {
  16       gAlice = (AliRun*)file->Get("gAlice");
  17       if (gAlice) printf("\nAliRun object found on file\n");
  18       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  19     }
  20   Int_t nbytes = 0;
  21   Int_t ipart;
  22   Int_t volume;
  23   Int_t np[nRings][nSectors];
  24
  25
  26   for (int i=0;i<nRings;i++)
  27     for(int j=0;j<nSectors;j++)
  28       np[i][j]=0;
  29
  30   TH1F *hNReal = new TH1F("hNReal","Real number of particles",50,0,50);
  31   TH1F *hNRec  = new TH1F("hNRec ","Reconst. number of particles",50,0,50);
  32
  33   Int_t nparticles = gAlice->GetEvent(0);
  34   if (nparticles <= 0) return;
  35   printf("\nnparticles=%d\n",nparticles);
  36
  37   gAlice->TreeR()->GetEvent(0);
  38   AliFMD *FMD  = (AliFMD*)gAlice->GetDetector("FMD");
  39   TClonesArray *Particles = gAlice->Particles();
  40   TParticle *particle;
  41   AliFMDhit *fmdHit;
  42   AliFMDReconstParticles *fmdRP;
  43   if (FMD)
  44     {
  45       TClonesArray *FMDhits   = FMD->Hits();
  46       TClonesArray *FMDrec    = FMD->ReconParticles();
  47       TTree *TH = gAlice->TreeH();
  48       Int_t ntracks    = TH->GetEntries();
  49       if (ntracks<=0) return;
  50
  51       Int_t  nPads=FMDrec->GetEntries();
  52
  53
  54 #ifdef DEBUG
  55       cout<<"\n(AliFMDReconstParticles*)FMDrec->UncheckedAt(0)="<<(AliFMDReconstParticles*)FMDrec->UncheckedAt(0);
  56       cout<<"\nFMDrec->UncheckedAt(0)="<<FMDrec->UncheckedAt(0);
  57 #endif
  58
  59       for (Int_t track=0; track<ntracks;track++)
  60         {
  61           gAlice->ResetHits();
  62           nbytes += TH->GetEvent(track);//?
  63           particle=(TParticle*)Particles->UncheckedAt(track);
  64           //      Int_t numpart=particle->GetKF();
  65           //Float_t eta=particle->GetEta();
  66
  67           Int_t  nhits=FMDhits->GetEntriesFast();
  68           for (Int_t hit=0;hit<nhits;hit++)
  69             {
  70               fmdHit  = (AliFMDhit*)FMDhits->UncheckedAt(hit);
  71               volume=fmdHit->Volume();
  72               if(volume==vol)
  73                 {
  74                   np[fmdHit->NumberOfRing()-1][fmdHit->NumberOfSector()-1]++;
  75                 }
  76             }
  77         }
  78       //Int_t  nRecPart=FMDrec->GetEntriesFast();
  79       Int_t nDeterm=0; Int_t nReal=0;
  80       for (Int_t pad=0;pad<nPads;pad++)
  81         {
  82           fmdRP  = (AliFMDReconstParticles*)FMDrec->UncheckedAt(pad);
  83           volume=fmdRP->GetVolume();
  84           if(volume==vol)
  85             {
  86 ;
  87 #ifdef DEBUG
  88               fmdDigit  = (AliFMDdigit*)FMDdig->UncheckedAt(pad);
  89               cout<<"\nfmdDigit->ADCsignal()="<<fmdDigit->ADCsignal();
  90               cout<<"\nfmdDigit->NumberOfRing()="<<fmdDigit->NumberOfRing();
  91               cout<<"\nfmdDigit->NumberOfSector()="<<fmdDigit->NumberOfSector();
  92 #endif
  93              nDeterm+=fmdRP->GetNumberOfReconstParticles();
  94              nReal+=np[fmdRP->GetNumberOfRing()-1][fmdRP->GetNumberOfSector()-1]; //-1=?
  95               Int_t RecRing=fmdRP->GetNumberOfRing()-1;
  96               Int_t RecSector=fmdRP->GetNumberOfSector()-1;
  97               hNReal->Fill(np[RecRing][RecSector]);
  98               hNRec->Fill(fmdRP->GetNumberOfReconstParticles());
  99             }
 100         }
 101     }
 102   cout<<"\nReal="<<nReal<<
 103     " nDeterm="<<nDeterm<<
 104     "\nerror="<<float(nDeterm-nReal)/float(nReal)<<endl;
 105
 106
 107   TCanvas *c1 = new TCanvas("c1","Alice FMD ",400,10,800,800);
 108   hNReal->SetFillColor(2);
 109   hNReal->Draw();
 110   hNRec->SetFillColor(4);
 111   hNRec->Draw("same");
 112
 113 }
 114
 115
 116
 117
 118
 119