iostream.h replaced by Riostream.h
[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"
30
dc8af42e 31// --- Standard library ---
37c55dc0 32#include <stdlib.h>
e75fab8f 33#include "Riostream.h"
dc8af42e 34
35// --- AliRoot header files ---
36
37#include "AliFMDdigit.h"
1a42d4f2 38#include "AliFMDhit.h"
dc8af42e 39#include "AliFMDReconstParticles.h"
40#include "AliFMD.h"
41#include "AliFMDv1.h"
42#include "AliFMDReconstruction.h"
43#include "AliRun.h"
1a42d4f2 44#include "AliHeader.h"
45#include "AliGenEventHeader.h"
dc8af42e 46ClassImp(AliFMDReconstruction)
47
48
49//____________________________________________________________________________
50
51AliFMDReconstruction::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
60AliFMDReconstruction::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
72AliFMDReconstruction::~AliFMDReconstruction()
73{
dc8af42e 74}
cb1df35e 75//----------------------------------------------------------------------------
dc8af42e 76
37c55dc0 77void AliFMDReconstruction::Exec(Option_t *option)
dc8af42e 78{
79 //Collects all digits in the same active volume into number of particles
1a42d4f2 80 cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
cb1df35e 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];
1a42d4f2 85 Float_t eta, etain,etaout,rad,theta;
cb1df35e 86 Int_t ivol, iSector, iRing;
1a42d4f2 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};
cb1df35e 91 Int_t NumberOfSectors[5]= {20,40,20,40,20};
1a42d4f2 92 Int_t NumberOfEtaIntervals[5];
cb1df35e 93 // number of ring for boundary 0.1 eta
1a42d4f2 94 TBranch *brDigits=0;
cb1df35e 95
1a42d4f2 96 AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");
37c55dc0 97 TClonesArray *fReconParticles=FMD->ReconParticles();
cb1df35e 98 if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
1a42d4f2 99 cout<<" fNevents "<<fNevents<<endl;
dc8af42e 100 for(Int_t ievent=0;ievent<fNevents;ievent++)
101 {
1a42d4f2 102 cout<<" ievent "<<ievent<<endl;
cb1df35e 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 ???
dc8af42e 107 gAlice->GetEvent(ievent) ;
108 if(gAlice->TreeH()==0) return;
1a42d4f2 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
dc8af42e 116 if(gAlice->TreeR()==0) gAlice->MakeTree("R");
117 //Make branches
118 FMD->MakeBranch("R");
1a42d4f2 119
37c55dc0 120
cb1df35e 121 Int_t zeroADC=1;
1a42d4f2 122
123 AliFMDdigit *fmdDigit;
124 if (FMD)
dc8af42e 125 {
126 gAlice->TreeD()->GetEvent(0);
127 TClonesArray *FMDdigits=FMD->Digits();
37c55dc0 128 Int_t nDigits=FMDdigits->GetEntries();
1a42d4f2 129 cout<<" nDigits "<<nDigits<<endl;
cb1df35e 130 Int_t RecParticles[6];
dc8af42e 131 Int_t nRecPart=0 ;
cb1df35e 132 Int_t ZeroPads=0;
133 Int_t NumberOfPads=0; //To avoid warning
134 Int_t pedestal;
135 Float_t channelWidth=(22400*50)/1024;
dc8af42e 136 for (Int_t digit=0;digit<nDigits;digit++)
137 {
cb1df35e 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;
dc8af42e 148 } //digit loop
cb1df35e 149 Int_t Rmin=0; Int_t Rmax=0; //To avoid warning
150 Int_t Smin=0; Int_t Smax=0; //To avoid warning
1a42d4f2 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++)
cb1df35e 157 {
1a42d4f2 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++)
cb1df35e 166 {
1a42d4f2 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]));
cb1df35e 174 ZeroPads=0;
175 Smin=0;
1a42d4f2 176 Smax=NumberOfSectors[ivol];
cb1df35e 177 for (Int_t iring=Rmin; iring<Rmax; iring++)
178 {
179 NumberOfPads=(Rmax-Rmin)*(Smax-Smin);
180 for
181 (Int_t isector=1;
1a42d4f2 182 isector<=NumberOfSectors[ivol];
cb1df35e 183 isector++)
1a42d4f2 184 if(PadADC[ivol][isector-1][iring-1]<zeroADC)
cb1df35e 185 ZeroPads++;
186 } //ring
cb1df35e 187 Float_t zerosRatio=
188 (Float_t)ZeroPads/(Float_t)NumberOfPads;
cb1df35e 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
1a42d4f2 197 (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
cb1df35e 198 else
199 RecParticles[5]=
200 Determination_by_thresholds
1a42d4f2 201 (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
cb1df35e 202 new((*fReconParticles)[nRecPart++])
203 AliFMDReconstParticles(RecParticles);
1a42d4f2 204 } // eta
205 } // volume
cb1df35e 206
cb1df35e 207 }//if FMD
1a42d4f2 208 gAlice->TreeR()->Reset();
209 gAlice->TreeR()->Fill();
210 gAlice->TreeR()->Write(0,TObject::kOverwrite);
dc8af42e 211 } //event loop
1a42d4f2 212 cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
cb1df35e 213};
dc8af42e 214
cb1df35e 215//------------------------------------------------------------
216Int_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)
dc8af42e 219{
1a42d4f2 220 cout<<"\nStart threshold method\n";
cb1df35e 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
1a42d4f2 259 cout<<"\nEnd threshol method"<<endl;
cb1df35e 260 return NumPart;
dc8af42e 261}
cb1df35e 262
263
264 //__________________________________________________________________
265
266Int_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)
dc8af42e 269{
1a42d4f2 270 // cout<<"\nStart Poisson method";
cb1df35e 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
dc8af42e 294
295
296
37c55dc0 297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
1a42d4f2 315
dc8af42e 316
317