Reconstruction in 0.1 eta over all sectors
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstruction.cxx
CommitLineData
37c55dc0 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
dc8af42e 23// --- ROOT system ---
24#include "TTask.h"
25#include "TTree.h"
26#include "TSystem.h"
27#include "TFile.h"
37c55dc0 28#include "TROOT.h"
29#include "TFolder.h"
cb1df35e 30#include <Riostream.h>
37c55dc0 31
dc8af42e 32// --- Standard library ---
37c55dc0 33#include <stdlib.h>
dc8af42e 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"
dc8af42e 43ClassImp(AliFMDReconstruction)
44
45
46//____________________________________________________________________________
47
48AliFMDReconstruction::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
57AliFMDReconstruction::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
69AliFMDReconstruction::~AliFMDReconstruction()
70{
dc8af42e 71}
cb1df35e 72//----------------------------------------------------------------------------
dc8af42e 73
37c55dc0 74void AliFMDReconstruction::Exec(Option_t *option)
dc8af42e 75{
cb1df35e 76 printf (" AliFMDReconstruction starting \n");
dc8af42e 77 //Collects all digits in the same active volume into number of particles
cb1df35e 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");
37c55dc0 105 TClonesArray *fReconParticles=FMD->ReconParticles();
cb1df35e 106 if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
dc8af42e 107 for(Int_t ievent=0;ievent<fNevents;ievent++)
108 {
cb1df35e 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 ???
dc8af42e 113 gAlice->GetEvent(ievent) ;
114 if(gAlice->TreeH()==0) return;
115 if(gAlice->TreeR()==0) gAlice->MakeTree("R");
116 //Make branches
cb1df35e 117 AliFMDdigit *fmdDigit;
dc8af42e 118 FMD->MakeBranch("R");
37c55dc0 119
cb1df35e 120 Int_t zeroADC=1;
121 // Int_t threshold_array_size=30;
37c55dc0 122
cb1df35e 123 // cout<<" AliFMDdigit "<<AliFMDdigit<<endl;
dc8af42e 124 if (FMD)
125 {
126 gAlice->TreeD()->GetEvent(0);
127 TClonesArray *FMDdigits=FMD->Digits();
37c55dc0 128 Int_t nDigits=FMDdigits->GetEntries();
cb1df35e 129 Int_t RecParticles[6];
dc8af42e 130 Int_t nRecPart=0 ;
cb1df35e 131 Int_t ZeroPads=0;
132 Int_t NumberOfPads=0; //To avoid warning
133 Int_t pedestal;
134 Float_t channelWidth=(22400*50)/1024;
dc8af42e 135 for (Int_t digit=0;digit<nDigits;digit++)
136 {
cb1df35e 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;
dc8af42e 147 } //digit loop
cb1df35e 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);
dc8af42e 220 } //event loop
cb1df35e 221 // cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
222};
dc8af42e 223
cb1df35e 224//------------------------------------------------------------
225Int_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)
dc8af42e 228{
cb1df35e 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;
dc8af42e 269}
cb1df35e 270
271
272 //__________________________________________________________________
273
274Int_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)
dc8af42e 277{
cb1df35e 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
dc8af42e 301
302
303
37c55dc0 304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
dc8af42e 322
323