1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // This is a TTask that constructs ReconstParticles (reconstructed particles)
20 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
21 //////////////////////////////////////////////////////////////////////////////
23 // --- ROOT system ---
31 // --- Standard library ---
33 #include <Riostream.h>
35 // --- AliRoot header files ---
37 #include "AliRunLoader.h"
38 #include "AliLoader.h"
40 #include "AliFMDdigit.h"
41 #include "AliFMDhit.h"
42 #include "AliFMDReconstParticles.h"
45 #include "AliFMDReconstruction.h"
47 #include "AliConfig.h"
48 #include "AliHeader.h"
49 #include "AliGenEventHeader.h"
51 ClassImp(AliFMDReconstruction)
54 //____________________________________________________________________________
56 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","")
58 fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
62 //____________________________________________________________________________
64 AliFMDReconstruction::AliFMDReconstruction(AliRunLoader* rl):TTask("AliFMDReconstruction","")
69 Fatal("AliFMDReconstruction","Argument AliRunLoader* is null!");
73 fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
76 AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
79 Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
80 return;//never reached
82 //add Task to //root/Tasks folder
83 gime->PostReconstructioner(this);
86 //____________________________________________________________________________
88 AliFMDReconstruction::~AliFMDReconstruction()
92 //____________________________________________________________________________
94 void AliFMDReconstruction::Exec(Option_t *option)
96 //Collects all digits in the same active volume into number of particles
98 Reconstruct number of particles
99 in given group of pads for given FMDvolume
100 determine by numberOfVolume ,
101 numberOfMinSector,numberOfMaxSector,
102 numberOfMinRing, numberOgMaxRing
103 Reconstruction method choose dependence on number of empty pads
107 cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
108 Int_t const knumVolumes=5;
109 Int_t const knumSectors=40;
110 Int_t const knumRings=768;
111 Int_t padADC[10][50][800];
112 Float_t eta, etain,etaout,rad,theta;
113 Int_t ivol, iSector, iRing;
114 Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
115 Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
116 Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
117 Int_t numberOfRings[5]={768,384,768,384,768};
118 Int_t numberOfSectors[5]= {20,40,20,40,20};
119 Int_t numberOfEtaIntervals[5];
120 // number of ring for boundary 0.1 eta
123 if (fRunLoader == 0x0)
125 Error("Exec","Run Loader loader is NULL - Session not opened");
128 AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
131 Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
132 return;//never reached
135 fRunLoader->LoadgAlice();
138 retval = gime->LoadHits("READ");
141 Error("Exec","Error occured while loading hits. Exiting.");
145 retval = gime->LoadDigits("READ");
148 Error("Exec","Error occured while loading digits. Exiting.");
152 AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
153 TClonesArray *fReconParticles=fFMD->ReconParticles();
155 TTree* treeD = gime->TreeD();
158 Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
161 if(fNevents == 0) fNevents=(Int_t)treeD->GetEntries();
162 //PH Do we use TreeE (kinematics), or TreeD (digits) toaccess the number
164 //PH if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
165 //PH cout<<" fNevents "<<fNevents<<endl;
169 for(Int_t ievent=0;ievent<fNevents;ievent++)
171 fRunLoader->GetEvent(ievent) ;
173 TTree* treeH = gime->TreeH();
176 Error("Exec","Can not get TreeH");
179 cout<<" ievent "<<ievent<<endl;
180 for (Int_t i=0; i<knumVolumes; i++)
181 for(Int_t j=0; j<knumSectors; j++)
182 for(Int_t ij=0; ij<knumRings; ij++)
183 padADC[i][j][ij]=0; //zhachem ???
185 brDigits=treeD->GetBranch("FMD");
187 cerr<<"EXEC Branch FMD digits not found"<<endl;
191 if(gime->TreeR()==0) gime->MakeTree("R");
194 fFMD->MakeBranch("R");
199 AliFMDdigit *fmdDigit;
202 gime->TreeD()->GetEvent(0);
203 TClonesArray *fFMDdigits=fFMD->Digits();
204 Int_t nDigits=fFMDdigits->GetEntries();
205 cout<<" nDigits "<<nDigits<<endl;
206 Int_t recParticles[6];
209 Int_t numberOfPads=0; //To avoid warning
211 Float_t channelWidth=(22400*50)/1024;
212 for (Int_t digit=0;digit<nDigits;digit++)
214 fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);
215 ivol=fmdDigit->Volume();
216 iSector=fmdDigit->NumberOfSector();
217 iRing=fmdDigit->NumberOfRing();
218 pedestal=Int_t(gRandom->Gaus(500,250));
219 padADC[ivol-1][iSector-1][iRing-1]=
220 fmdDigit->ADCsignal()
221 -Int_t(Float_t(pedestal)/channelWidth);
222 if (padADC[ivol-1][iSector-1][iRing-1]<0)
223 padADC[ivol-1][iSector-1][iRing-1]=0;
225 Int_t rmin=0; Int_t rmax=0; //To avoid warning
226 Int_t smin=0; Int_t smax=0; //To avoid warning
227 AliHeader *header = fRunLoader->GetHeader();
228 AliGenEventHeader* genHeader = header->GenEventHeader();
229 TArrayF *o = new TArrayF(3);
230 genHeader->PrimaryVertex(*o);
231 Float_t zVertex=o->At(2);
232 for (ivol=0; ivol<knumVolumes; ivol++)
234 Float_t realZ=zVertex+z[ivol];
235 theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
236 etain = - TMath::Log( TMath::Tan(theta/2.));
237 theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
238 etaout=- TMath::Log( TMath::Tan(theta/2.));
239 numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
241 for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
243 theta = 2.*TMath::ATan (TMath::Exp (-eta));
244 Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
245 rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
247 theta = 2.*TMath::ATan (TMath::Exp (-eta));
248 rad = TMath::Abs(realZ) * (TMath::Tan (theta));
249 rmax=Int_t( (rad-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
252 smax=numberOfSectors[ivol];
253 for (Int_t iring=rmin; iring<rmax; iring++)
255 numberOfPads=(rmax-rmin)*(smax-smin);
258 isector<=numberOfSectors[ivol];
260 if(padADC[ivol][isector-1][iring-1]<zeroADC)
264 (Float_t)zeroPads/(Float_t)numberOfPads;
265 recParticles[0]=ivol;
266 recParticles[1]=smin;
267 recParticles[2]=smax;
268 recParticles[3]=rmin;
269 recParticles[4]=rmax;
272 DeterminationByPoisson
273 (padADC,ivol+1,rmin,rmax,smin,smax);
276 DeterminationByThresholds
277 (padADC,ivol+1,rmin,rmax,smin,smax);
278 new((*fReconParticles)[nRecPart++])
279 AliFMDReconstParticles(recParticles);
284 gime->TreeR()->Reset();
285 gime->TreeR()->Fill();
286 gime->WriteRecPoints("OVERWRITE");
288 cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
291 //------------------------------------------------------------
292 Int_t AliFMDReconstruction::DeterminationByThresholds
293 (Int_t padAdc[][50][800],Int_t volume, Int_t rmin, Int_t rmax,
294 Int_t smin, Int_t smax)
297 reconstruct number of particles by
298 energy deposited threshold method
299 Using if number of empty pads less then 10%
302 cout<<"\nStart threshold method\n";
304 Int_t thresholdInner[30]={
306 72, 89, 104, 124, 129,
307 164, 174, 179, 208, 228,
308 248, 268, 317, 337, 357,
309 392, 407, 416, 426, 436,
310 461, 468, 493, 506, 515};
312 Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
319 for (Int_t it=0; it<30; it++) {
320 if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
321 if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
323 Int_t thresholdArraySize = 30;
326 for (Int_t iring=rmin; iring<rmax; iring++)
328 for (Int_t isector=smin; isector<smax; isector++)
330 for (int i=0;i<thresholdArraySize-1;i++)
332 if(padAdc[volume-1][isector][iring]>threshold[i]
333 &&padAdc[volume-1][isector][iring]<=threshold[i+1])
337 }; // if in threshol interval
338 }; //threshold_array_size
341 cout<<"\nEnd threshol method"<<endl;
345 //__________________________________________________________________
347 Int_t AliFMDReconstruction::DeterminationByPoisson
348 (Int_t padAdc[][50][800],
349 Int_t vol, Int_t rmin, Int_t rmax,
350 Int_t secmin, Int_t secmax)
353 reconstruct number of particles by Poisson statistic method
354 according number of empty pad in chosen group of pads
355 using if number of empty pads more then 10%
359 // cout<<"\nStart Poisson method";
360 Int_t thresholdADC=1;
362 for (Int_t i=rmin;i<rmax;i++)
364 for (Int_t j=secmin;j<secmax;j++)
366 if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
369 Float_t lambda=-TMath::Log(Float_t(zeropads)/
370 ( Float_t(secmax-secmin)*
371 Float_t(rmax-rmin))); //+1 zdes' ne nado
372 Int_t fRecon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
373 // cout<<"\nEnd Poisson method"<<endl;