49480abbf4dd816382798580f2546bc5760a5541
[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 <Riostream.h>
31
32 // --- Standard library ---
33 #include <stdlib.h>
34
35 // --- AliRoot header files ---
36
37 #include "AliFMDdigit.h"
38 #include "AliFMDReconstParticles.h"
39 #include "AliFMD.h"
40 #include "AliFMDv1.h"
41 #include "AliFMDReconstruction.h"
42 #include "AliRun.h"
43 ClassImp(AliFMDReconstruction)
44
45         
46 //____________________________________________________________________________ 
47
48 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","") 
49 {
50   fNevents = 0 ;  // Number of events to rreconnstraction, 0 means all events in current file
51   // add Task to //root/Tasks folder
52   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
53   roottasks->Add(this) ; 
54 }
55 //____________________________________________________________________________ 
56
57 AliFMDReconstruction::AliFMDReconstruction(char* HeaderFile, char *ReconstParticlesFile):TTask("AliFMDReconstruction","")
58 {
59   fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file
60   fReconstParticlesFile=ReconstParticlesFile ;
61   fHeadersFile=HeaderFile ;
62   //add Task to //root/Tasks folder
63   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
64   roottasks->Add(this) ;     
65 }
66
67 //____________________________________________________________________________ 
68
69 AliFMDReconstruction::~AliFMDReconstruction()
70 {
71 }
72 //----------------------------------------------------------------------------
73
74 void AliFMDReconstruction::Exec(Option_t *option) 
75
76   printf (" AliFMDReconstruction starting \n");
77  //Collects all digits in the same active volume into number of particles
78   
79   Int_t const NumVolums=5;
80   Int_t const NumSectors=40;
81   Int_t const NumRings=768;
82   Int_t PadADC[10][50][800];
83   Int_t ivol, iSector, iRing;
84   Int_t Ne1; 
85   //Int_t NumberOfRings[5]=  {256,128,256,128,256};
86   Int_t NumberOfSectors[5]=  {20,40,20,40,20};
87   // number of ring for boundary 0.1 eta
88    Int_t EtaIntervalInner []=
89       {0, 55, 110, 165, 221, 276, 331, 386, 442,
90         497, 552, 607, 663, 718, 767 };
91    /*
92      {0,  18, 36,   55,  73,
93       92, 110, 128, 147, 165,
94       184, 202, 221, 239, 255};*/
95    Int_t  EtaIntervalOuter []=  //{0, 21, 43, 65, 86, 108, 127};
96      {0, 65, 130, 195, 260, 325, 383};
97  Int_t EtaInterval[20];
98   
99
100   Int_t size_EtaIntervalInner=sizeof (EtaIntervalInner)/sizeof(EtaIntervalInner[0]);
101
102   Int_t size_EtaIntervalOuter=sizeof (EtaIntervalOuter)/sizeof(EtaIntervalOuter[0]);
103
104   AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");  
105   TClonesArray *fReconParticles=FMD->ReconParticles();
106   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries(); 
107   for(Int_t ievent=0;ievent<fNevents;ievent++)
108     { 
109       for (Int_t i=0; i<NumVolums; i++)
110         for(Int_t j=0; j<NumSectors; j++)
111           for(Int_t ij=0; ij<NumRings; ij++)
112             PadADC[i][j][ij]=0;                    //zhachem ???
113       gAlice->GetEvent(ievent) ;
114       if(gAlice->TreeH()==0) return; 
115       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
116       //Make branches
117       AliFMDdigit  *fmdDigit;
118       FMD->MakeBranch("R");
119       
120       Int_t zeroADC=1;
121       // Int_t  threshold_array_size=30;
122
123       // cout<<" AliFMDdigit "<<AliFMDdigit<<endl;
124       if (FMD)
125         {
126           gAlice->TreeD()->GetEvent(0); 
127           TClonesArray *FMDdigits=FMD->Digits();
128           Int_t nDigits=FMDdigits->GetEntries();
129            Int_t RecParticles[6];
130            Int_t nRecPart=0 ;
131            Int_t ZeroPads=0;
132            Int_t NumberOfPads=0; //To avoid warning
133            Int_t pedestal;
134            Float_t channelWidth=(22400*50)/1024;
135            for (Int_t digit=0;digit<nDigits;digit++) 
136              {
137                fmdDigit=(AliFMDdigit*)FMDdigits->UncheckedAt(digit);    
138                ivol=fmdDigit->Volume();
139                iSector=fmdDigit->NumberOfSector();
140                iRing=fmdDigit->NumberOfRing();
141                pedestal=Int_t(gRandom->Gaus(500,250));
142                PadADC[ivol-1][iSector-1][iRing-1]=
143                  fmdDigit->ADCsignal()
144                  -Int_t(Float_t(pedestal)/channelWidth);
145                if (PadADC[ivol-1][iSector-1][iRing-1]<0) 
146                  PadADC[ivol-1][iSector-1][iRing-1]=0;
147              } //digit loop
148            Int_t Rmin=0; Int_t Rmax=0; //To avoid warning
149            Int_t Smin=0; Int_t Smax=0; //To avoid warning
150                    
151            for (ivol=1; ivol<=NumVolums; ivol++)
152              {
153                if (ivol==1||ivol==3||ivol==5)
154                  {
155                    Ne1=size_EtaIntervalInner;
156                    for(Int_t ieta=0; ieta<Ne1; ieta++)
157                      EtaInterval[ieta]=EtaIntervalInner[ieta];
158                  }
159                if (ivol==2||ivol==4)
160                  {
161                    Ne1=size_EtaIntervalOuter;
162                    for( Int_t ieta=0; ieta<Ne1; ieta++)
163                      EtaInterval[ieta]=EtaIntervalOuter[ieta];
164                  }
165
166
167                    for (Int_t e1=0;e1<Ne1-1;e1++) // vol<=NumVolums
168                      {
169                        Rmin=EtaInterval[e1];
170                    Rmax=EtaInterval[e1+1];
171                    ZeroPads=0;
172                    Smin=0;
173                    Smax=NumberOfSectors[ivol-1]; 
174                    for (Int_t iring=Rmin; iring<Rmax; iring++) 
175                      {
176                        NumberOfPads=(Rmax-Rmin)*(Smax-Smin);
177                        for 
178                          (Int_t isector=1;
179                           isector<=NumberOfSectors[ivol-1]; 
180                           isector++) 
181                          if(PadADC[ivol-1][isector-1][iring-1]<zeroADC)
182                            ZeroPads++;
183                      } //ring
184                    /*
185                    cout<<"\nRmin="<<Rmin;
186                    cout<<"\nRmax="<<Rmax;
187                    cout<<"\nSmin="<<Smin;
188                    cout<<"\nSmax="<<Smax;
189                
190                    cout<<"\nvolume "<<ivol<<" zero "<<ZeroPads<<
191                      " NumberOfPads "<<NumberOfPads<<endl;
192                    */
193                    Float_t zerosRatio= 
194                      (Float_t)ZeroPads/(Float_t)NumberOfPads;
195                    //  cout<<"\nzerosRatio="<<zerosRatio;                     
196                    RecParticles[0]=ivol;
197                    RecParticles[1]=Smin;
198                    RecParticles[2]=Smax;
199                    RecParticles[3]=Rmin;
200                    RecParticles[4]=Rmax;
201                    if (zerosRatio>0.1 ||(ivol==2||ivol==4))
202                      RecParticles[5]=
203                        Determination_by_Poisson
204                        (PadADC,ivol,Rmin,Rmax,Smin,Smax);
205                    else
206                      RecParticles[5]=
207                        Determination_by_thresholds
208                        (PadADC,ivol,Rmin,Rmax,Smin,Smax);
209                    //  cout<<" nDeterm "<<RecParticles[5]<<endl;
210                    new((*fReconParticles)[nRecPart++]) 
211                  AliFMDReconstParticles(RecParticles);                 
212         } // eta
213     } // volume
214            
215            // if(zerosRatio<0.01) Determination_by_counting (ZeroPads) ;           
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::Determination_by_thresholds
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   Int_t thresholdInner[30]={
231     0,     14,  28,    42,   57,     
232     72,    89,  104,  124,  129, 
233     164,  174,  179,  208,  228, 
234     248,  268,   317,  337,  357, 
235     392,  407,  416,  426,  436, 
236     461,  468,  493,  506,  515}; 
237
238   Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
239                             132, 165, 198, 231, 
240                             264, 286, 308, 334, 
241                             352, 374, 418, 440,
242                             462, 484, 506, 528,
243                             550, 572, 594, 616}; 
244   Int_t threshold[30];
245   for (Int_t it=0; it<30; it++) {
246     if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
247     if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
248   }
249   Int_t threshold_array_size = 30;     
250   Int_t NumPart=0;
251   //Inner Si
252   for (Int_t iring=Rmin; iring<Rmax; iring++) 
253     {
254       for (Int_t isector=Smin; isector<Smax; isector++) 
255         {
256           for (int i=0;i<threshold_array_size-1;i++)
257             {         
258               if(PadADC[volume-1][isector][iring]>threshold[i]
259                  &&PadADC[volume-1][isector][iring]<=threshold[i+1])
260                 {
261                   NumPart+=i;
262                   break;
263                 }; // if in threshol interval
264             }; //threshold_array_size
265         }; //iring
266     }; //sector
267   // cout<<"\nEnd threshol method"<<endl;
268   return NumPart;
269 }
270
271
272  //__________________________________________________________________
273
274 Int_t AliFMDReconstruction::Determination_by_Poisson (Int_t PadADC[][50][800],
275                                            Int_t vol, Int_t rmin, Int_t rmax, 
276                                            Int_t secmin, Int_t secmax)
277 {
278   Int_t threshold_adc=1;   
279   Int_t zeropads=0;
280   for (Int_t i=rmin;i<rmax;i++)
281     {
282       for (Int_t j=secmin;j<secmax;j++)
283         {
284           if (PadADC[vol-1][j][i]<threshold_adc) zeropads++;
285         };
286     };
287   Float_t lambda=-TMath::Log(Float_t(zeropads)/
288                              ( Float_t(secmax-secmin)*
289                                Float_t(rmax-rmin))); //+1 zdes' ne nado
290   Int_t Recon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
291   //  cout<<"\nEnd Poisson method"<<endl;
292   return Recon;
293 };
294
295 //------------------------------------------------------------------
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323