code clean
[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
31 // --- Standard library ---
32 #include <stdlib.h>
33 #include "Riostream.h"
34
35 // --- AliRoot header files ---
36
37 #include "AliFMDdigit.h"
38 #include "AliFMDhit.h"
39 #include "AliFMDReconstParticles.h"
40 #include "AliFMD.h"
41 #include "AliFMDv1.h"
42 #include "AliFMDReconstruction.h"
43 #include "AliRun.h"
44 #include "AliHeader.h"
45 #include "AliGenEventHeader.h"
46 ClassImp(AliFMDReconstruction)
47
48         
49 //____________________________________________________________________________ 
50
51 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","") 
52 {
53   fNevents = 0 ;  // Number of events to rreconnstraction, 0 means all events in current file
54   // add Task to //root/Tasks folder
55   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
56   roottasks->Add(this) ; 
57 }
58 //____________________________________________________________________________ 
59
60 AliFMDReconstruction::AliFMDReconstruction(char* HeaderFile, char *ReconstParticlesFile):TTask("AliFMDReconstruction","")
61 {
62   fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file
63   fReconstParticlesFile=ReconstParticlesFile ;
64   fHeadersFile=HeaderFile ;
65   //add Task to //root/Tasks folder
66   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
67   roottasks->Add(this) ;     
68 }
69
70 //____________________________________________________________________________ 
71
72 AliFMDReconstruction::~AliFMDReconstruction()
73 {
74 }
75 //----------------------------------------------------------------------------
76
77 void AliFMDReconstruction::Exec(Option_t *option) 
78
79   /*
80     Reconstruct nember of particles 
81     in given group of pads for given FMDvolume
82     determine by numberOfVolume , 
83     numberOfMinSector,numberOfMaxSector,
84     numberOfMinRing, numberOgMaxRing
85     Reconstruction method choose dependence on number of empty pads  
86   */
87
88
89   cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
90   Int_t const knumVolumes=5;
91   Int_t const knumSectors=40;
92   Int_t const knumRings=768;
93   Int_t padADC[10][50][800];
94   Float_t eta, etain,etaout,rad,theta;
95   Int_t ivol, iSector, iRing;
96   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
97   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
98   Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
99   Int_t numberOfRings[5]={768,384,768,384,768};
100   Int_t numberOfSectors[5]=  {20,40,20,40,20};
101   Int_t numberOfEtaIntervals[5];
102   // number of ring for boundary 0.1 eta
103   TBranch *brDigits=0;
104
105   AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
106   TClonesArray *fReconParticles=fFMD->ReconParticles();
107   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries(); 
108   cout<<" fNevents "<<fNevents<<endl;
109   for(Int_t ievent=0;ievent<fNevents;ievent++)
110     { 
111       cout<<" ievent "<<ievent<<endl;
112       for (Int_t i=0; i<knumVolumes; i++)
113         for(Int_t j=0; j<knumSectors; j++)
114           for(Int_t ij=0; ij<knumRings; ij++)
115             padADC[i][j][ij]=0;                    //zhachem ???
116       gAlice->GetEvent(ievent) ;
117       if(gAlice->TreeH()==0) return; 
118       if(gAlice->TreeD()==0) return;
119         brDigits=gAlice->TreeD()->GetBranch("FMD");
120         if (!brDigits){
121           cerr<<"EXEC Branch FMD digits not found"<<endl;
122         return;
123       } 
124  
125       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
126       //Make branches
127       fFMD->MakeBranch("R");
128
129       
130       Int_t zeroADC=1;
131  
132       AliFMDdigit  *fmdDigit;
133        if (fFMD)
134         {
135           gAlice->TreeD()->GetEvent(0); 
136           TClonesArray *fFMDdigits=fFMD->Digits();
137           Int_t nDigits=fFMDdigits->GetEntries();
138           cout<<" nDigits "<<nDigits<<endl;
139            Int_t recParticles[6];
140            Int_t nRecPart=0 ;
141            Int_t zeroPads=0;
142            Int_t numberOfPads=0; //To avoid warning
143            Int_t pedestal;
144            Float_t channelWidth=(22400*50)/1024;
145            for (Int_t digit=0;digit<nDigits;digit++) 
146              {
147                fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);    
148                ivol=fmdDigit->Volume();
149                iSector=fmdDigit->NumberOfSector();
150                iRing=fmdDigit->NumberOfRing();
151                pedestal=Int_t(gRandom->Gaus(500,250));
152                padADC[ivol-1][iSector-1][iRing-1]=
153                  fmdDigit->ADCsignal()
154                  -Int_t(Float_t(pedestal)/channelWidth);
155                if (padADC[ivol-1][iSector-1][iRing-1]<0) 
156                  padADC[ivol-1][iSector-1][iRing-1]=0;
157              } //digit loop
158            Int_t rmin=0; Int_t rmax=0; //To avoid warning
159            Int_t smin=0; Int_t smax=0; //To avoid warning
160            AliHeader *header = gAlice->GetHeader();
161            AliGenEventHeader* genHeader = header->GenEventHeader();
162            TArrayF *o = new TArrayF(3); 
163            genHeader->PrimaryVertex(*o);
164            Float_t zVertex=o->At(2);
165            for (ivol=0; ivol<knumVolumes; ivol++)
166              {
167                Float_t realZ=zVertex+z[ivol];
168                theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
169                etain = - TMath::Log( TMath::Tan(theta/2.));
170                theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
171                etaout=- TMath::Log( TMath::Tan(theta/2.));
172                numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
173                eta=etain;
174                for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++) 
175                  {
176                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
177                    Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
178                    rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
179                    eta=eta-0.1;
180                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
181                    rad = TMath::Abs(realZ) * (TMath::Tan (theta));
182                    rmax=Int_t( (rad-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
183                    zeroPads=0;
184                    smin=0;
185                    smax=numberOfSectors[ivol]; 
186                    for (Int_t iring=rmin; iring<rmax; iring++) 
187                      {
188                        numberOfPads=(rmax-rmin)*(smax-smin);
189                        for 
190                          (Int_t isector=1;
191                           isector<=numberOfSectors[ivol]; 
192                           isector++) 
193                          if(padADC[ivol][isector-1][iring-1]<zeroADC)
194                            zeroPads++;
195                      } //ring
196                    Float_t zerosRatio= 
197                      (Float_t)zeroPads/(Float_t)numberOfPads;
198                    recParticles[0]=ivol;
199                    recParticles[1]=smin;
200                    recParticles[2]=smax;
201                    recParticles[3]=rmin;
202                    recParticles[4]=rmax;
203                    if (zerosRatio>0.1 )
204                      recParticles[5]=
205                        DeterminationByPoisson
206                        (padADC,ivol+1,rmin,rmax,smin,smax);
207                    else
208                      recParticles[5]=
209                        DeterminationByThresholds
210                        (padADC,ivol+1,rmin,rmax,smin,smax);
211                    new((*fReconParticles)[nRecPart++]) 
212                  AliFMDReconstParticles(recParticles);                 
213                  } // eta
214              } // volume
215            
216         }//if FMD
217        gAlice->TreeR()->Reset();
218        gAlice->TreeR()->Fill(); 
219        gAlice->TreeR()->Write(0,TObject::kOverwrite);
220     } //event loop
221   cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
222 };
223
224 //------------------------------------------------------------
225 Int_t AliFMDReconstruction::DeterminationByThresholds
226 (Int_t padAdc[][50][800],Int_t volume, Int_t rmin, Int_t rmax, 
227  Int_t smin, Int_t smax)
228 {
229   /*
230     reconstruct number of particles by 
231     energy deposited threshold method 
232     Using if number of empty pads less then 10%
233
234 */
235   cout<<"\nStart threshold method\n";
236   
237   Int_t thresholdInner[30]={
238     0,     14,  28,    42,   57,     
239     72,    89,  104,  124,  129, 
240     164,  174,  179,  208,  228, 
241     248,  268,   317,  337,  357, 
242     392,  407,  416,  426,  436, 
243     461,  468,  493,  506,  515}; 
244
245   Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
246                             132, 165, 198, 231, 
247                             264, 286, 308, 334, 
248                             352, 374, 418, 440,
249                             462, 484, 506, 528,
250                             550, 572, 594, 616}; 
251   Int_t threshold[30];
252   for (Int_t it=0; it<30; it++) {
253     if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
254     if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
255   }
256   Int_t thresholdArraySize = 30;     
257   Int_t numPart=0;
258   //Inner Si
259   for (Int_t iring=rmin; iring<rmax; iring++) 
260     {
261       for (Int_t isector=smin; isector<smax; isector++) 
262         {
263           for (int i=0;i<thresholdArraySize-1;i++)
264             {         
265               if(padAdc[volume-1][isector][iring]>threshold[i]
266                  &&padAdc[volume-1][isector][iring]<=threshold[i+1])
267                 {
268                   numPart+=i;
269                   break;
270                 }; // if in threshol interval
271             }; //threshold_array_size
272         }; //iring
273     }; //sector
274   cout<<"\nEnd threshol method"<<endl;
275   return numPart;
276 }
277
278
279  //__________________________________________________________________
280
281 Int_t AliFMDReconstruction::DeterminationByPoisson 
282 (Int_t padAdc[][50][800],
283  Int_t vol, Int_t rmin, Int_t rmax, 
284  Int_t secmin, Int_t secmax)
285 {
286   /* 
287      reconstruct number of particles by Poisson statistic method 
288      according number of empty pad in chosen group of pads 
289      using if number of empty pads more then 10% 
290
291   */
292
293   //  cout<<"\nStart Poisson method";
294   Int_t thresholdADC=1;   
295   Int_t zeropads=0;
296   for (Int_t i=rmin;i<rmax;i++)
297     {
298       for (Int_t j=secmin;j<secmax;j++)
299         {
300           if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
301         };
302     };
303   Float_t lambda=-TMath::Log(Float_t(zeropads)/
304                              ( Float_t(secmax-secmin)*
305                                Float_t(rmax-rmin))); //+1 zdes' ne nado
306   Int_t fRecon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
307   //  cout<<"\nEnd Poisson method"<<endl;
308   return fRecon;
309 };
310
311 //------------------------------------------------------------------
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340