Adaption to new fluka common blocks (E. Futo)
[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{
a39f3926 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
1a42d4f2 89 cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
a39f3926 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];
1a42d4f2 94 Float_t eta, etain,etaout,rad,theta;
cb1df35e 95 Int_t ivol, iSector, iRing;
1a42d4f2 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.};
a39f3926 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];
cb1df35e 102 // number of ring for boundary 0.1 eta
1a42d4f2 103 TBranch *brDigits=0;
cb1df35e 104
a39f3926 105 AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
106 TClonesArray *fReconParticles=fFMD->ReconParticles();
cb1df35e 107 if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
1a42d4f2 108 cout<<" fNevents "<<fNevents<<endl;
dc8af42e 109 for(Int_t ievent=0;ievent<fNevents;ievent++)
110 {
1a42d4f2 111 cout<<" ievent "<<ievent<<endl;
a39f3926 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 ???
dc8af42e 116 gAlice->GetEvent(ievent) ;
117 if(gAlice->TreeH()==0) return;
1a42d4f2 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
dc8af42e 125 if(gAlice->TreeR()==0) gAlice->MakeTree("R");
126 //Make branches
a39f3926 127 fFMD->MakeBranch("R");
1a42d4f2 128
37c55dc0 129
cb1df35e 130 Int_t zeroADC=1;
1a42d4f2 131
132 AliFMDdigit *fmdDigit;
a39f3926 133 if (fFMD)
dc8af42e 134 {
135 gAlice->TreeD()->GetEvent(0);
a39f3926 136 TClonesArray *fFMDdigits=fFMD->Digits();
137 Int_t nDigits=fFMDdigits->GetEntries();
1a42d4f2 138 cout<<" nDigits "<<nDigits<<endl;
a39f3926 139 Int_t recParticles[6];
dc8af42e 140 Int_t nRecPart=0 ;
a39f3926 141 Int_t zeroPads=0;
142 Int_t numberOfPads=0; //To avoid warning
cb1df35e 143 Int_t pedestal;
144 Float_t channelWidth=(22400*50)/1024;
dc8af42e 145 for (Int_t digit=0;digit<nDigits;digit++)
146 {
a39f3926 147 fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);
cb1df35e 148 ivol=fmdDigit->Volume();
149 iSector=fmdDigit->NumberOfSector();
150 iRing=fmdDigit->NumberOfRing();
151 pedestal=Int_t(gRandom->Gaus(500,250));
a39f3926 152 padADC[ivol-1][iSector-1][iRing-1]=
cb1df35e 153 fmdDigit->ADCsignal()
154 -Int_t(Float_t(pedestal)/channelWidth);
a39f3926 155 if (padADC[ivol-1][iSector-1][iRing-1]<0)
156 padADC[ivol-1][iSector-1][iRing-1]=0;
dc8af42e 157 } //digit loop
a39f3926 158 Int_t rmin=0; Int_t rmax=0; //To avoid warning
159 Int_t smin=0; Int_t smax=0; //To avoid warning
1a42d4f2 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);
a39f3926 165 for (ivol=0; ivol<knumVolumes; ivol++)
cb1df35e 166 {
1a42d4f2 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.));
a39f3926 172 numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
1a42d4f2 173 eta=etain;
a39f3926 174 for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
cb1df35e 175 {
1a42d4f2 176 theta = 2.*TMath::ATan (TMath::Exp (-eta));
177 Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
a39f3926 178 rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
1a42d4f2 179 eta=eta-0.1;
180 theta = 2.*TMath::ATan (TMath::Exp (-eta));
181 rad = TMath::Abs(realZ) * (TMath::Tan (theta));
a39f3926 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++)
cb1df35e 187 {
a39f3926 188 numberOfPads=(rmax-rmin)*(smax-smin);
cb1df35e 189 for
190 (Int_t isector=1;
a39f3926 191 isector<=numberOfSectors[ivol];
cb1df35e 192 isector++)
a39f3926 193 if(padADC[ivol][isector-1][iring-1]<zeroADC)
194 zeroPads++;
cb1df35e 195 } //ring
cb1df35e 196 Float_t zerosRatio=
a39f3926 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);
cb1df35e 207 else
a39f3926 208 recParticles[5]=
209 DeterminationByThresholds
210 (padADC,ivol+1,rmin,rmax,smin,smax);
cb1df35e 211 new((*fReconParticles)[nRecPart++])
a39f3926 212 AliFMDReconstParticles(recParticles);
1a42d4f2 213 } // eta
214 } // volume
cb1df35e 215
cb1df35e 216 }//if FMD
1a42d4f2 217 gAlice->TreeR()->Reset();
218 gAlice->TreeR()->Fill();
219 gAlice->TreeR()->Write(0,TObject::kOverwrite);
dc8af42e 220 } //event loop
1a42d4f2 221 cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
cb1df35e 222};
dc8af42e 223
cb1df35e 224//------------------------------------------------------------
a39f3926 225Int_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)
dc8af42e 228{
a39f3926 229 /*
230 reconstruct number of particles by
231 energy deposited threshold method
232 Using if number of empty pads less then 10%
233
234*/
1a42d4f2 235 cout<<"\nStart threshold method\n";
cb1df35e 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 }
a39f3926 256 Int_t thresholdArraySize = 30;
257 Int_t numPart=0;
cb1df35e 258 //Inner Si
a39f3926 259 for (Int_t iring=rmin; iring<rmax; iring++)
cb1df35e 260 {
a39f3926 261 for (Int_t isector=smin; isector<smax; isector++)
cb1df35e 262 {
a39f3926 263 for (int i=0;i<thresholdArraySize-1;i++)
cb1df35e 264 {
a39f3926 265 if(padAdc[volume-1][isector][iring]>threshold[i]
266 &&padAdc[volume-1][isector][iring]<=threshold[i+1])
cb1df35e 267 {
a39f3926 268 numPart+=i;
cb1df35e 269 break;
270 }; // if in threshol interval
271 }; //threshold_array_size
272 }; //iring
273 }; //sector
1a42d4f2 274 cout<<"\nEnd threshol method"<<endl;
a39f3926 275 return numPart;
dc8af42e 276}
cb1df35e 277
278
279 //__________________________________________________________________
280
a39f3926 281Int_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)
dc8af42e 285{
a39f3926 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
1a42d4f2 293 // cout<<"\nStart Poisson method";
a39f3926 294 Int_t thresholdADC=1;
cb1df35e 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 {
a39f3926 300 if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
cb1df35e 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
a39f3926 306 Int_t fRecon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
cb1df35e 307 // cout<<"\nEnd Poisson method"<<endl;
a39f3926 308 return fRecon;
cb1df35e 309};
310
311//------------------------------------------------------------------
312
313
314
315
316
dc8af42e 317
318
319
37c55dc0 320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
1a42d4f2 338
dc8af42e 339
340