AliFMDMerger.h
[u/mrichter/AliRoot.git] / 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