Transition to NewIO
[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>
88cb7938 33#include <Riostream.h>
dc8af42e 34
35// --- AliRoot header files ---
36
88cb7938 37#include "AliRunLoader.h"
38#include "AliLoader.h"
39
dc8af42e 40#include "AliFMDdigit.h"
1a42d4f2 41#include "AliFMDhit.h"
dc8af42e 42#include "AliFMDReconstParticles.h"
43#include "AliFMD.h"
44#include "AliFMDv1.h"
45#include "AliFMDReconstruction.h"
46#include "AliRun.h"
88cb7938 47#include "AliConfig.h"
1a42d4f2 48#include "AliHeader.h"
49#include "AliGenEventHeader.h"
88cb7938 50
dc8af42e 51ClassImp(AliFMDReconstruction)
52
53
54//____________________________________________________________________________
55
56AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","")
57{
58 fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
88cb7938 59 fRunLoader = 0x0;
60
dc8af42e 61}
62//____________________________________________________________________________
63
88cb7938 64AliFMDReconstruction::AliFMDReconstruction(AliRunLoader* rl):TTask("AliFMDReconstruction","")
dc8af42e 65{
88cb7938 66
67 if (rl == 0x0)
68 {
69 Fatal("AliFMDReconstruction","Argument AliRunLoader* is null!");
70 return;
71 }
72
dc8af42e 73 fNevents = 0 ; // Number of events to rreconnstraction, 0 means all events in current file
88cb7938 74
75 fRunLoader = rl;
76 AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
77 if (gime == 0x0)
78 {
79 Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
80 return;//never reached
81 }
dc8af42e 82 //add Task to //root/Tasks folder
88cb7938 83 gime->PostReconstructioner(this);
dc8af42e 84}
85
86//____________________________________________________________________________
87
88AliFMDReconstruction::~AliFMDReconstruction()
89{
dc8af42e 90}
88cb7938 91
92//____________________________________________________________________________
dc8af42e 93
37c55dc0 94void AliFMDReconstruction::Exec(Option_t *option)
dc8af42e 95{
88cb7938 96 //Collects all digits in the same active volume into number of particles
a39f3926 97 /*
88cb7938 98 Reconstruct number of particles
a39f3926 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
104 */
105
106
1a42d4f2 107 cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
a39f3926 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];
1a42d4f2 112 Float_t eta, etain,etaout,rad,theta;
cb1df35e 113 Int_t ivol, iSector, iRing;
1a42d4f2 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.};
a39f3926 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];
cb1df35e 120 // number of ring for boundary 0.1 eta
cb1df35e 121
88cb7938 122
123 if (fRunLoader == 0x0)
124 {
125 Error("Exec","Run Loader loader is NULL - Session not opened");
126 return;
127 }
128 AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
129 if (gime == 0x0)
130 {
131 Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
132 return;//never reached
133 }
134
135 fRunLoader->LoadgAlice();
136 Int_t retval;
137
138 retval = gime->LoadHits("READ");
139 if (retval)
140 {
141 Error("Exec","Error occured while loading hits. Exiting.");
142 return;
143 }
144
145 retval = gime->LoadDigits("READ");
146 if (retval)
147 {
148 Error("Exec","Error occured while loading digits. Exiting.");
149 return;
150 }
151
a39f3926 152 AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
153 TClonesArray *fReconParticles=fFMD->ReconParticles();
88cb7938 154
155 TTree* treeD = gime->TreeD();
156 if (treeD == 0x0)
157 {
158 Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
159 return;
160 }
161 if(fNevents == 0) fNevents=(Int_t)treeD->GetEntries();
162 //PH Do we use TreeE (kinematics), or TreeD (digits) toaccess the number
163 //PH of events?
164//PH if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
165//PH cout<<" fNevents "<<fNevents<<endl;
166
167 TBranch *brDigits=0;
168
dc8af42e 169 for(Int_t ievent=0;ievent<fNevents;ievent++)
170 {
88cb7938 171 fRunLoader->GetEvent(ievent) ;
172
173 TTree* treeH = gime->TreeH();
174 if (treeH == 0x0)
175 {
176 Error("Exec","Can not get TreeH");
177 return;
178 }
1a42d4f2 179 cout<<" ievent "<<ievent<<endl;
a39f3926 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 ???
88cb7938 184
185 brDigits=treeD->GetBranch("FMD");
1a42d4f2 186 if (!brDigits){
187 cerr<<"EXEC Branch FMD digits not found"<<endl;
188 return;
189 }
190
88cb7938 191 if(gime->TreeR()==0) gime->MakeTree("R");
192
dc8af42e 193 //Make branches
a39f3926 194 fFMD->MakeBranch("R");
1a42d4f2 195
37c55dc0 196
cb1df35e 197 Int_t zeroADC=1;
1a42d4f2 198
199 AliFMDdigit *fmdDigit;
a39f3926 200 if (fFMD)
dc8af42e 201 {
88cb7938 202 gime->TreeD()->GetEvent(0);
a39f3926 203 TClonesArray *fFMDdigits=fFMD->Digits();
204 Int_t nDigits=fFMDdigits->GetEntries();
1a42d4f2 205 cout<<" nDigits "<<nDigits<<endl;
a39f3926 206 Int_t recParticles[6];
dc8af42e 207 Int_t nRecPart=0 ;
a39f3926 208 Int_t zeroPads=0;
209 Int_t numberOfPads=0; //To avoid warning
cb1df35e 210 Int_t pedestal;
211 Float_t channelWidth=(22400*50)/1024;
dc8af42e 212 for (Int_t digit=0;digit<nDigits;digit++)
213 {
a39f3926 214 fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);
cb1df35e 215 ivol=fmdDigit->Volume();
216 iSector=fmdDigit->NumberOfSector();
217 iRing=fmdDigit->NumberOfRing();
218 pedestal=Int_t(gRandom->Gaus(500,250));
a39f3926 219 padADC[ivol-1][iSector-1][iRing-1]=
cb1df35e 220 fmdDigit->ADCsignal()
221 -Int_t(Float_t(pedestal)/channelWidth);
a39f3926 222 if (padADC[ivol-1][iSector-1][iRing-1]<0)
223 padADC[ivol-1][iSector-1][iRing-1]=0;
dc8af42e 224 } //digit loop
a39f3926 225 Int_t rmin=0; Int_t rmax=0; //To avoid warning
226 Int_t smin=0; Int_t smax=0; //To avoid warning
88cb7938 227 AliHeader *header = fRunLoader->GetHeader();
1a42d4f2 228 AliGenEventHeader* genHeader = header->GenEventHeader();
229 TArrayF *o = new TArrayF(3);
230 genHeader->PrimaryVertex(*o);
231 Float_t zVertex=o->At(2);
a39f3926 232 for (ivol=0; ivol<knumVolumes; ivol++)
cb1df35e 233 {
1a42d4f2 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.));
a39f3926 239 numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
1a42d4f2 240 eta=etain;
a39f3926 241 for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
cb1df35e 242 {
1a42d4f2 243 theta = 2.*TMath::ATan (TMath::Exp (-eta));
244 Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
a39f3926 245 rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
1a42d4f2 246 eta=eta-0.1;
247 theta = 2.*TMath::ATan (TMath::Exp (-eta));
248 rad = TMath::Abs(realZ) * (TMath::Tan (theta));
a39f3926 249 rmax=Int_t( (rad-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
250 zeroPads=0;
251 smin=0;
252 smax=numberOfSectors[ivol];
253 for (Int_t iring=rmin; iring<rmax; iring++)
cb1df35e 254 {
a39f3926 255 numberOfPads=(rmax-rmin)*(smax-smin);
cb1df35e 256 for
257 (Int_t isector=1;
a39f3926 258 isector<=numberOfSectors[ivol];
cb1df35e 259 isector++)
a39f3926 260 if(padADC[ivol][isector-1][iring-1]<zeroADC)
261 zeroPads++;
cb1df35e 262 } //ring
cb1df35e 263 Float_t zerosRatio=
a39f3926 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;
270 if (zerosRatio>0.1 )
271 recParticles[5]=
272 DeterminationByPoisson
273 (padADC,ivol+1,rmin,rmax,smin,smax);
cb1df35e 274 else
a39f3926 275 recParticles[5]=
276 DeterminationByThresholds
277 (padADC,ivol+1,rmin,rmax,smin,smax);
cb1df35e 278 new((*fReconParticles)[nRecPart++])
a39f3926 279 AliFMDReconstParticles(recParticles);
1a42d4f2 280 } // eta
281 } // volume
cb1df35e 282
cb1df35e 283 }//if FMD
88cb7938 284 gime->TreeR()->Reset();
285 gime->TreeR()->Fill();
286 gime->WriteRecPoints("OVERWRITE");
dc8af42e 287 } //event loop
1a42d4f2 288 cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
88cb7938 289}
dc8af42e 290
cb1df35e 291//------------------------------------------------------------
a39f3926 292Int_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)
dc8af42e 295{
a39f3926 296 /*
297 reconstruct number of particles by
298 energy deposited threshold method
299 Using if number of empty pads less then 10%
300
301*/
1a42d4f2 302 cout<<"\nStart threshold method\n";
cb1df35e 303
304 Int_t thresholdInner[30]={
305 0, 14, 28, 42, 57,
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};
311
312 Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
313 132, 165, 198, 231,
314 264, 286, 308, 334,
315 352, 374, 418, 440,
316 462, 484, 506, 528,
317 550, 572, 594, 616};
318 Int_t threshold[30];
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];
322 }
a39f3926 323 Int_t thresholdArraySize = 30;
324 Int_t numPart=0;
cb1df35e 325 //Inner Si
a39f3926 326 for (Int_t iring=rmin; iring<rmax; iring++)
cb1df35e 327 {
a39f3926 328 for (Int_t isector=smin; isector<smax; isector++)
cb1df35e 329 {
a39f3926 330 for (int i=0;i<thresholdArraySize-1;i++)
cb1df35e 331 {
a39f3926 332 if(padAdc[volume-1][isector][iring]>threshold[i]
333 &&padAdc[volume-1][isector][iring]<=threshold[i+1])
cb1df35e 334 {
a39f3926 335 numPart+=i;
cb1df35e 336 break;
337 }; // if in threshol interval
338 }; //threshold_array_size
339 }; //iring
340 }; //sector
1a42d4f2 341 cout<<"\nEnd threshol method"<<endl;
a39f3926 342 return numPart;
dc8af42e 343}
cb1df35e 344
88cb7938 345//__________________________________________________________________
cb1df35e 346
a39f3926 347Int_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)
dc8af42e 351{
a39f3926 352 /*
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%
356
357 */
358
1a42d4f2 359 // cout<<"\nStart Poisson method";
a39f3926 360 Int_t thresholdADC=1;
cb1df35e 361 Int_t zeropads=0;
362 for (Int_t i=rmin;i<rmax;i++)
363 {
364 for (Int_t j=secmin;j<secmax;j++)
365 {
a39f3926 366 if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
cb1df35e 367 };
368 };
369 Float_t lambda=-TMath::Log(Float_t(zeropads)/
370 ( Float_t(secmax-secmin)*
371 Float_t(rmax-rmin))); //+1 zdes' ne nado
a39f3926 372 Int_t fRecon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
cb1df35e 373 // cout<<"\nEnd Poisson method"<<endl;
a39f3926 374 return fRecon;
88cb7938 375}
37c55dc0 376
1a42d4f2 377