iostream.h replaced by Riostream.h
[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  //Collects all digits in the same active volume into number of particles
80   cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
81   Int_t const NumVolums=5;
82   Int_t const NumSectors=40;
83   Int_t const NumRings=768;
84   Int_t PadADC[10][50][800];
85   Float_t eta, etain,etaout,rad,theta;
86   Int_t ivol, iSector, iRing;
87   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
88   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
89   Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
90   Int_t NumberOfRings[5]={768,384,768,384,768};
91   Int_t NumberOfSectors[5]=  {20,40,20,40,20};
92   Int_t NumberOfEtaIntervals[5];
93   // number of ring for boundary 0.1 eta
94   TBranch *brDigits=0;
95
96   AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");
97   TClonesArray *fReconParticles=FMD->ReconParticles();
98   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries(); 
99   cout<<" fNevents "<<fNevents<<endl;
100   for(Int_t ievent=0;ievent<fNevents;ievent++)
101     { 
102       cout<<" ievent "<<ievent<<endl;
103       for (Int_t i=0; i<NumVolums; i++)
104         for(Int_t j=0; j<NumSectors; j++)
105           for(Int_t ij=0; ij<NumRings; ij++)
106             PadADC[i][j][ij]=0;                    //zhachem ???
107       gAlice->GetEvent(ievent) ;
108       if(gAlice->TreeH()==0) return; 
109       if(gAlice->TreeD()==0) return;
110         brDigits=gAlice->TreeD()->GetBranch("FMD");
111         if (!brDigits){
112           cerr<<"EXEC Branch FMD digits not found"<<endl;
113         return;
114       } 
115  
116       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
117       //Make branches
118       FMD->MakeBranch("R");
119
120       
121       Int_t zeroADC=1;
122  
123       AliFMDdigit  *fmdDigit;
124        if (FMD)
125         {
126           gAlice->TreeD()->GetEvent(0); 
127           TClonesArray *FMDdigits=FMD->Digits();
128           Int_t nDigits=FMDdigits->GetEntries();
129           cout<<" nDigits "<<nDigits<<endl;
130            Int_t RecParticles[6];
131            Int_t nRecPart=0 ;
132            Int_t ZeroPads=0;
133            Int_t NumberOfPads=0; //To avoid warning
134            Int_t pedestal;
135            Float_t channelWidth=(22400*50)/1024;
136            for (Int_t digit=0;digit<nDigits;digit++) 
137              {
138                fmdDigit=(AliFMDdigit*)FMDdigits->UncheckedAt(digit);    
139                ivol=fmdDigit->Volume();
140                iSector=fmdDigit->NumberOfSector();
141                iRing=fmdDigit->NumberOfRing();
142                pedestal=Int_t(gRandom->Gaus(500,250));
143                PadADC[ivol-1][iSector-1][iRing-1]=
144                  fmdDigit->ADCsignal()
145                  -Int_t(Float_t(pedestal)/channelWidth);
146                if (PadADC[ivol-1][iSector-1][iRing-1]<0) 
147                  PadADC[ivol-1][iSector-1][iRing-1]=0;
148              } //digit loop
149            Int_t Rmin=0; Int_t Rmax=0; //To avoid warning
150            Int_t Smin=0; Int_t Smax=0; //To avoid warning
151            AliHeader *header = gAlice->GetHeader();
152            AliGenEventHeader* genHeader = header->GenEventHeader();
153            TArrayF *o = new TArrayF(3); 
154            genHeader->PrimaryVertex(*o);
155            Float_t zVertex=o->At(2);
156            for (ivol=0; ivol<NumVolums; ivol++)
157              {
158                Float_t realZ=zVertex+z[ivol];
159                theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
160                etain = - TMath::Log( TMath::Tan(theta/2.));
161                theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
162                etaout=- TMath::Log( TMath::Tan(theta/2.));
163                NumberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
164                eta=etain;
165                for (Int_t e1=0;e1<=NumberOfEtaIntervals[ivol];e1++) 
166                  {
167                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
168                    Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
169                    Rmin= Int_t ( (radmin-rin[ivol])*NumberOfRings[ivol]/(rout[ivol]-rin[ivol]));
170                    eta=eta-0.1;
171                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
172                    rad = TMath::Abs(realZ) * (TMath::Tan (theta));
173                    Rmax=Int_t( (rad-rin[ivol])*NumberOfRings[ivol]/(rout[ivol]-rin[ivol]));
174                    ZeroPads=0;
175                    Smin=0;
176                    Smax=NumberOfSectors[ivol]; 
177                    for (Int_t iring=Rmin; iring<Rmax; iring++) 
178                      {
179                        NumberOfPads=(Rmax-Rmin)*(Smax-Smin);
180                        for 
181                          (Int_t isector=1;
182                           isector<=NumberOfSectors[ivol]; 
183                           isector++) 
184                          if(PadADC[ivol][isector-1][iring-1]<zeroADC)
185                            ZeroPads++;
186                      } //ring
187                    Float_t zerosRatio= 
188                      (Float_t)ZeroPads/(Float_t)NumberOfPads;
189                    RecParticles[0]=ivol;
190                    RecParticles[1]=Smin;
191                    RecParticles[2]=Smax;
192                    RecParticles[3]=Rmin;
193                    RecParticles[4]=Rmax;
194                    if (zerosRatio>0.1 ||(ivol==2||ivol==4))
195                      RecParticles[5]=
196                        Determination_by_Poisson
197                        (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
198                    else
199                      RecParticles[5]=
200                        Determination_by_thresholds
201                        (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
202                    new((*fReconParticles)[nRecPart++]) 
203                  AliFMDReconstParticles(RecParticles);                 
204                  } // eta
205              } // volume
206            
207         }//if FMD
208        gAlice->TreeR()->Reset();
209        gAlice->TreeR()->Fill(); 
210        gAlice->TreeR()->Write(0,TObject::kOverwrite);
211     } //event loop
212   cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
213 };
214
215 //------------------------------------------------------------
216 Int_t AliFMDReconstruction::Determination_by_thresholds
217 (Int_t PadADC[][50][800],Int_t volume, Int_t Rmin, Int_t Rmax, 
218  Int_t Smin, Int_t Smax)
219 {
220   cout<<"\nStart threshold method\n";
221   
222   Int_t thresholdInner[30]={
223     0,     14,  28,    42,   57,     
224     72,    89,  104,  124,  129, 
225     164,  174,  179,  208,  228, 
226     248,  268,   317,  337,  357, 
227     392,  407,  416,  426,  436, 
228     461,  468,  493,  506,  515}; 
229
230   Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
231                             132, 165, 198, 231, 
232                             264, 286, 308, 334, 
233                             352, 374, 418, 440,
234                             462, 484, 506, 528,
235                             550, 572, 594, 616}; 
236   Int_t threshold[30];
237   for (Int_t it=0; it<30; it++) {
238     if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
239     if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
240   }
241   Int_t threshold_array_size = 30;     
242   Int_t NumPart=0;
243   //Inner Si
244   for (Int_t iring=Rmin; iring<Rmax; iring++) 
245     {
246       for (Int_t isector=Smin; isector<Smax; isector++) 
247         {
248           for (int i=0;i<threshold_array_size-1;i++)
249             {         
250               if(PadADC[volume-1][isector][iring]>threshold[i]
251                  &&PadADC[volume-1][isector][iring]<=threshold[i+1])
252                 {
253                   NumPart+=i;
254                   break;
255                 }; // if in threshol interval
256             }; //threshold_array_size
257         }; //iring
258     }; //sector
259   cout<<"\nEnd threshol method"<<endl;
260   return NumPart;
261 }
262
263
264  //__________________________________________________________________
265
266 Int_t AliFMDReconstruction::Determination_by_Poisson (Int_t PadADC[][50][800],
267                                            Int_t vol, Int_t rmin, Int_t rmax, 
268                                            Int_t secmin, Int_t secmax)
269 {
270   //  cout<<"\nStart Poisson method";
271   Int_t threshold_adc=1;   
272   Int_t zeropads=0;
273   for (Int_t i=rmin;i<rmax;i++)
274     {
275       for (Int_t j=secmin;j<secmax;j++)
276         {
277           if (PadADC[vol-1][j][i]<threshold_adc) zeropads++;
278         };
279     };
280   Float_t lambda=-TMath::Log(Float_t(zeropads)/
281                              ( Float_t(secmax-secmin)*
282                                Float_t(rmax-rmin))); //+1 zdes' ne nado
283   Int_t Recon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
284   //  cout<<"\nEnd Poisson method"<<endl;
285   return Recon;
286 };
287
288 //------------------------------------------------------------------
289
290
291
292
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