reconstruction working now under AliReconstruction
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstruction.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 // This is a TTask that constructs ReconstParticles (reconstructed particles) 
18 // out of Digits
19 // 
20 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
21 //////////////////////////////////////////////////////////////////////////////
22
23 // --- ROOT system ---
24 #include "TTask.h"
25 #include "TTree.h"
26 #include "TSystem.h"
27 #include "TFile.h"
28 #include "TROOT.h"
29 #include "TFolder.h"
30 #include "TH2F.h"
31
32 // --- Standard library ---
33 #include <stdlib.h>
34 #include <Riostream.h>
35
36 // --- AliRoot header files ---
37
38 #include "AliRunLoader.h"
39 #include "AliLoader.h"
40
41 #include "AliFMDdigit.h"
42 #include "AliFMDhit.h"
43 #include "AliFMDReconstParticles.h"
44 #include "AliFMD.h"
45 #include "AliFMDv1.h"
46 #include "AliFMDReconstruction.h"
47 #include "AliRun.h"
48 #include "AliConfig.h"
49 #include "AliHeader.h"
50 #include "AliGenEventHeader.h"
51
52 ClassImp(AliFMDReconstruction)
53
54         
55 //____________________________________________________________________________ 
56
57 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","") 
58 {
59   fNevents = 0 ;  // Number of events to rreconnstraction, 0 means all events in current file
60   fRunLoader = 0x0;
61   
62 }
63 //____________________________________________________________________________ 
64
65 AliFMDReconstruction::AliFMDReconstruction(AliRunLoader* rl):TTask("AliFMDReconstruction","")
66 {
67
68   if (rl == 0x0)
69    {
70      Fatal("AliFMDReconstruction","Argument AliRunLoader* is null!");
71      return;
72    }
73    
74   fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file
75
76   fRunLoader = rl;
77   AliLoader* plFMD = fRunLoader->GetLoader("FMDLoader");
78   if (plFMD == 0x0)
79    {
80      Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
81      return;//never reached
82    }
83   //add Task to //root/Tasks folder
84 //  gime->PostReconstructioner(this);
85 }
86
87 //____________________________________________________________________________ 
88
89 AliFMDReconstruction::~AliFMDReconstruction()
90 {
91   AliLoader* plFMD = fRunLoader->GetLoader("FMDLoader");
92   plFMD->CleanReconstructioner();
93 }
94
95 //____________________________________________________________________________
96
97 void AliFMDReconstruction::Exec() 
98
99  //Collects all digits in the same active volume into number of particles
100   /*
101     Reconstruct number of particles 
102     in given group of pads for given FMDvolume
103     determine by numberOfVolume , 
104     numberOfMinSector,numberOfMaxSector,
105     numberOfMinRing, numberOgMaxRing
106     Reconstruction method choose dependence on number of empty pads  
107   */
108
109
110 #ifdef DEBUG
111   Info("Exec ","Start");
112 #endif
113
114   Int_t const knumVolumes=5;
115   Int_t padADC;
116   Float_t eta, etain,etaout,rad,theta;
117   Int_t ivol, iSector, iRing;
118   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
119   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
120   Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
121   Int_t numberOfRings[5]={512,256,512,256,512};
122   Int_t numberOfSectors[5]=  {20,40,20,40,20};
123   Int_t numberOfEtaIntervals[5];
124   // number of ring for boundary 0.1 eta
125
126   
127   if (fRunLoader == 0x0)
128    {
129     Error("Exec","Run Loader loader is NULL - Session not opened");
130     return;
131    }
132
133   AliLoader* plFMD = fRunLoader->GetLoader("FMDLoader");
134   if (plFMD == 0x0)
135    {
136      Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
137      return;//never reached
138    }
139    
140   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
141   if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
142
143   TDirectory* cwd = gDirectory;
144   gDirectory = 0x0;
145   Text_t buf1[20];
146   TH2F* hTotal[10];
147   for (Int_t j=1; j<=5; j++){
148     sprintf(buf1,"hTotal%d",j);
149     
150     hTotal[j] = new TH2F(buf1," Number of primary particles ",
151                          numberOfSectors[j-1],1,numberOfSectors[j-1],
152                          numberOfRings[j-1],1,numberOfRings[j-1]);
153   }   
154   gDirectory = cwd;
155  
156   Int_t retval=0;     
157   if(fNevents == 0) fNevents=Int_t (fRunLoader->TreeE()->GetEntries()); 
158 #ifdef DEBUG
159   cout<<" fNevents "<<fNevents<<endl;
160 #endif
161    for(Int_t ievent=0;ievent<fNevents;ievent++)
162     { 
163 #ifdef DEBUG
164       cout<<" *** event "<<ievent<<endl; 
165 #endif
166       fRunLoader->GetEvent(ievent) ;
167       //event z-vertex for correction eta-rad dependence      
168       AliHeader *header = fRunLoader->GetHeader();
169       AliGenEventHeader* genHeader = header->GenEventHeader();
170       TArrayF *o = new TArrayF(3); 
171       genHeader->PrimaryVertex(*o);
172       Float_t zVertex=o->At(2);
173
174       for (Int_t i=0; i<5; i++)
175         hTotal[i+1]->Reset();
176       
177       retval = plFMD->LoadDigits("READ"); 
178       if (retval)
179         {
180           Error("Exec","Error occured while loading digits. Exiting.");
181           return;
182         }
183       
184       AliFMD * fFMD = (AliFMD *)gAlice->GetDetector("FMD");
185       TClonesArray *fReconParticles=fFMD->ReconParticles();
186       TClonesArray *fDigits=fFMD->Digits();
187  
188       TTree* treeD = plFMD->TreeD();
189       if (treeD == 0x0)
190         {
191           Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
192           return;
193         }
194       
195       TBranch *brDigits=0;
196             
197       brDigits=treeD->GetBranch("FMD");
198
199       if (brDigits) {
200         brDigits->SetAddress(&fDigits);
201       }else{
202         cerr<<"EXEC Branch FMD digits not found"<<endl;
203         return;
204       } 
205           
206       if(plFMD->TreeR()==0) plFMD->MakeTree("R");
207
208       //Make branches
209       fFMD->MakeBranch("R");
210
211       
212       Int_t zeroADC=6;
213       // read Digits 
214       AliFMDdigit  *fmdDigit;
215        if (fFMD)
216         {
217           plFMD->TreeD()->GetEvent(0); 
218                   
219           Int_t nDigits=fDigits->GetEntries();
220           Int_t recParticles[6];
221           Int_t nRecPart=0 ;
222           Int_t zeroPads=0;
223           Int_t numberOfPads=0;
224           Int_t pedestal;
225           Float_t channelWidth=(22400*50)/1024;
226           for (Int_t digit=0;digit<nDigits;digit++) 
227             {
228               fmdDigit=(AliFMDdigit*)fDigits->UncheckedAt(digit);    
229               ivol=fmdDigit->Volume();
230               iSector=fmdDigit->NumberOfSector();
231               iRing=fmdDigit->NumberOfRing();
232               pedestal=Int_t(gRandom->Gaus(500,250));
233               padADC= fmdDigit->ADCsignal()
234                 -Int_t(Float_t(pedestal)/channelWidth);
235               if (padADC<0) padADC=0;
236               hTotal[ivol]->Fill(iSector,iRing,padADC);
237             } //digit loop
238
239           //reconstruct multiplicity in 0.1 eta according Poisson distribution 
240
241           Int_t rmin=0; Int_t rmax=0; 
242           Int_t smin=0; Int_t smax=0; 
243           for (ivol=0; ivol<knumVolumes; ivol++)
244             {
245               Float_t ring2number=Float_t (numberOfRings[ivol])/
246                 (rout[ivol]-rin[ivol]);
247               Float_t realZ=zVertex+z[ivol];
248               theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
249               etain = - TMath::Log( TMath::Tan(theta/2.));
250               theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
251               etaout=- TMath::Log( TMath::Tan(theta/2.));
252               numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
253               eta=etain;
254               for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++) 
255                 {
256                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
257                   Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
258                   rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
259                   eta=eta+0.1;
260                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
261                   rad = TMath::Abs(realZ) * (TMath::Tan (theta));
262                   rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
263                   
264                   zeroPads=0;
265                   smin=0;
266                   smax=numberOfSectors[ivol]; 
267                   numberOfPads=(rmax-rmin)*(smax-smin);
268                   for (Int_t iring=rmin; iring<rmax; iring++) 
269                     {
270                       for 
271                         (Int_t isector=0;
272                          isector<numberOfSectors[ivol]; 
273                          isector++) 
274                         {
275                           if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
276                              <zeroADC) zeroPads++;}
277                       
278                     } //ring
279                   Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
280
281                   Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
282                   Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
283                   recParticles[0]=ivol+1;
284                   recParticles[1]=smin;
285                   recParticles[2]=smax;
286                   recParticles[3]=rmin;
287                   recParticles[4]=rmax;
288                   recParticles[5]= fRecon;
289                   new((*fReconParticles)[nRecPart++])   AliFMDReconstParticles(recParticles);             
290                   
291
292                  } // eta
293              } // volume
294            
295         }//if FMD
296        plFMD->TreeR()->Reset();
297        plFMD->TreeR()->Fill(); 
298        plFMD->WriteRecPoints("OVERWRITE");
299        plFMD->UnloadDigits();
300     } //event loop
301 #ifdef DEBUG
302   Info(" Exec"," finished");
303 #endif
304   //  delete hTotal[10]; 
305
306 }
307