Additional protection
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstructor.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 class 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 "AliFMDReconstructor.h"
47 #include "AliRun.h"
48 #include "AliConfig.h"
49 #include "AliHeader.h"
50 #include "AliGenEventHeader.h"
51
52 ClassImp(AliFMDReconstructor)
53
54         
55 //____________________________________________________________________________
56
57 void AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
58
59  //Collects all digits in the same active volume into number of particles
60   /*
61     Reconstruct number of particles 
62     in given group of pads for given FMDvolume
63     determine by numberOfVolume , 
64     numberOfMinSector,numberOfMaxSector,
65     numberOfMinRing, numberOgMaxRing
66     Reconstructor method choose dependence on number of empty pads  
67   */
68
69
70 #ifdef DEBUG
71   Info("Exec ","Start");
72 #endif
73
74   Int_t const knumVolumes=5;
75   Int_t padADC;
76   Float_t eta, etain,etaout,rad,theta;
77   Int_t ivol, iSector, iRing;
78   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
79   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
80   Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
81   Int_t numberOfRings[5]={512,256,512,256,512};
82   Int_t numberOfSectors[5]=  {20,40,20,40,20};
83   Int_t numberOfEtaIntervals[5];
84   // number of ring for boundary 0.1 eta
85
86   
87   if (runLoader == 0x0)
88    {
89     Error("Exec","Run Loader loader is NULL - Session not opened");
90     return;
91    }
92
93   AliLoader* plFMD = runLoader->GetLoader("FMDLoader");
94   if (plFMD == 0x0)
95    {
96      Fatal("AliFMDReconstructor","Can not find FMD (loader) in specified event");
97      return;//never reached
98    }
99    
100   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
101   if (!runLoader->TreeE()) runLoader->LoadHeader();
102
103   TDirectory* cwd = gDirectory;
104   gDirectory = 0x0;
105   Text_t buf1[20];
106   TH2F* hTotal[10];
107   for (Int_t j=1; j<=5; j++){
108     sprintf(buf1,"hTotal%d",j);
109     
110     hTotal[j] = new TH2F(buf1," Number of primary particles ",
111                          numberOfSectors[j-1],1,numberOfSectors[j-1],
112                          numberOfRings[j-1],1,numberOfRings[j-1]);
113   }   
114   gDirectory = cwd;
115  
116   plFMD->LoadRecPoints("RECREATE");
117   Int_t retval=0;     
118   Int_t nevents=Int_t (runLoader->TreeE()->GetEntries()); 
119 #ifdef DEBUG
120   cout<<" nevents "<<nevents<<endl;
121 #endif
122    for(Int_t ievent=0;ievent<nevents;ievent++)
123     { 
124 #ifdef DEBUG
125       cout<<" *** event "<<ievent<<endl; 
126 #endif
127       runLoader->GetEvent(ievent) ;
128       //event z-vertex for correction eta-rad dependence      
129       AliHeader *header = runLoader->GetHeader();
130       AliGenEventHeader* genHeader = header->GenEventHeader();
131       TArrayF *o = new TArrayF(3); 
132       if (genHeader) genHeader->PrimaryVertex(*o);
133       Float_t zVertex=o->At(2);
134
135       for (Int_t i=0; i<5; i++)
136         hTotal[i+1]->Reset();
137       
138       retval = plFMD->LoadDigits("READ"); 
139       if (retval)
140         {
141           Error("Exec","Error occured while loading digits. Exiting.");
142           return;
143         }
144       
145       AliFMD * fFMD = (AliFMD *)gAlice->GetDetector("FMD");
146       TClonesArray *fReconParticles=fFMD->ReconParticles();
147       TClonesArray *fDigits=fFMD->Digits();
148  
149       TTree* treeD = plFMD->TreeD();
150       if (treeD == 0x0)
151         {
152           Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
153           return;
154         }
155       
156       TBranch *brDigits=0;
157             
158       brDigits=treeD->GetBranch("FMD");
159
160       if (brDigits) {
161         brDigits->SetAddress(&fDigits);
162       }else{
163         cerr<<"EXEC Branch FMD digits not found"<<endl;
164         return;
165       } 
166           
167       if(plFMD->TreeR()==0) plFMD->MakeTree("R");
168
169       //Make branches
170       fFMD->MakeBranch("R");
171
172       
173       Int_t zeroADC=6;
174       // read Digits 
175       AliFMDdigit  *fmdDigit;
176        if (fFMD)
177         {
178           plFMD->TreeD()->GetEvent(0); 
179                   
180           Int_t nDigits=fDigits->GetEntries();
181           Int_t recParticles[6];
182           Int_t nRecPart=0 ;
183           Int_t zeroPads=0;
184           Int_t numberOfPads=0;
185           Int_t pedestal;
186           Float_t channelWidth=(22400*50)/1024;
187           for (Int_t digit=0;digit<nDigits;digit++) 
188             {
189               fmdDigit=(AliFMDdigit*)fDigits->UncheckedAt(digit);    
190               ivol=fmdDigit->Volume();
191               iSector=fmdDigit->NumberOfSector();
192               iRing=fmdDigit->NumberOfRing();
193               pedestal=Int_t(gRandom->Gaus(500,250));
194               padADC= fmdDigit->ADCsignal()
195                 -Int_t(Float_t(pedestal)/channelWidth);
196               if (padADC<0) padADC=0;
197               hTotal[ivol]->Fill(iSector,iRing,padADC);
198             } //digit loop
199
200           //reconstruct multiplicity in 0.1 eta according Poisson distribution 
201
202           Int_t rmin=0; Int_t rmax=0; 
203           Int_t smin=0; Int_t smax=0; 
204           for (ivol=0; ivol<knumVolumes; ivol++)
205             {
206               Float_t ring2number=Float_t (numberOfRings[ivol])/
207                 (rout[ivol]-rin[ivol]);
208               Float_t realZ=zVertex+z[ivol];
209               theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
210               etain = - TMath::Log( TMath::Tan(theta/2.));
211               theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
212               etaout=- TMath::Log( TMath::Tan(theta/2.));
213               numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
214               eta=etain;
215               for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++) 
216                 {
217                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
218                   Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
219                   rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
220                   eta=eta+0.1;
221                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
222                   rad = TMath::Abs(realZ) * (TMath::Tan (theta));
223                   rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
224                   
225                   zeroPads=0;
226                   smin=0;
227                   smax=numberOfSectors[ivol]; 
228                   numberOfPads=(rmax-rmin)*(smax-smin);
229                   for (Int_t iring=rmin; iring<rmax; iring++) 
230                     {
231                       for 
232                         (Int_t isector=0;
233                          isector<numberOfSectors[ivol]; 
234                          isector++) 
235                         {
236                           if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
237                              <zeroADC) zeroPads++;}
238                       
239                     } //ring
240                   Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
241
242                   Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
243                   Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
244                   recParticles[0]=ivol+1;
245                   recParticles[1]=smin;
246                   recParticles[2]=smax;
247                   recParticles[3]=rmin;
248                   recParticles[4]=rmax;
249                   recParticles[5]= fRecon;
250                   new((*fReconParticles)[nRecPart++])   AliFMDReconstParticles(recParticles);             
251                   
252
253                  } // eta
254              } // volume
255            
256         }//if FMD
257        plFMD->TreeR()->Reset();
258        plFMD->TreeR()->Fill(); 
259        plFMD->WriteRecPoints("OVERWRITE");
260        plFMD->UnloadDigits();
261     } //event loop
262   plFMD->UnloadRecPoints();
263 #ifdef DEBUG
264   Info(" Exec"," finished");
265 #endif
266   //  delete hTotal[10]; 
267
268 }
269
270
271 //_____________________________________________________________________________
272 void AliFMDReconstructor::FillESD(AliRunLoader* /*runLoader*/, 
273                                   AliESD* /*esd*/) const
274 {
275 // nothing to be done
276
277 }
278