]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCFast.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity fix
[u/mrichter/AliRoot.git] / TPC / AliTPCFast.cxx
CommitLineData
6d75e4b6 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/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// fast TPC cluster simulation //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
24#include <TParticle.h>
e8d02863 25#include <TVector.h>
e939a978 26#include <TRandom.h>
6d75e4b6 27
28#include "AliRunLoader.h"
29#include "AliRun.h"
30#include "AliMC.h"
31#include "AliTPC.h"
32#include "AliTPCParam.h"
33#include "AliTPCClustersArray.h"
34#include "AliTPCClustersRow.h"
35#include "AliTPCcluster.h"
36#include "AliComplexCluster.h"
37#include "AliTPCFast.h"
8c2b3fd7 38#include "AliLog.h"
6d75e4b6 39
40ClassImp(AliTPCFast)
179c6296 41 //____________________________________________________________________
42AliTPCFast::AliTPCFast(const AliTPCFast &param)
43 :TObject(param),fParam(0)
6d75e4b6 44
179c6296 45{
46 //
47 // copy constructor - dummy
48 //
49 fParam = param.fParam;
50}
51AliTPCFast & AliTPCFast::operator =(const AliTPCFast & param)
52{
53 //
54 // assignment operator - dummy
55 //
56 fParam=param.fParam;
57 return (*this);
58}
6d75e4b6 59
60//_____________________________________________________________________________
61void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
62{
63 //--------------------------------------------------------
64 // TPC simple cluster generator from hits
65 // obtained from the TPC Fast Simulator
66 // The point errors are taken from the parametrization
67 //--------------------------------------------------------
68
69 //-----------------------------------------------------------------
70 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
71 //-----------------------------------------------------------------
72 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
73 // Jouri.Belikov@cern.ch
74 //----------------------------------------------------------------
75
76 /////////////////////////////////////////////////////////////////////////////
77 //
78 //---------------------------------------------------------------------
79 // ALICE TPC Cluster Parameters
80 //--------------------------------------------------------------------
81
82
83
84 // Cluster width in rphi
85 const Float_t kACrphi=0.18322;
86 const Float_t kBCrphi=0.59551e-3;
87 const Float_t kCCrphi=0.60952e-1;
88 // Cluster width in z
89 const Float_t kACz=0.19081;
90 const Float_t kBCz=0.55938e-3;
91 const Float_t kCCz=0.30428;
92
93
94 AliLoader* loader = runLoader->GetLoader("TPCLoader");
95 if (!loader) {
8c2b3fd7 96 AliError("No TPC loader found");
6d75e4b6 97 return;
98 }
99 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
100 AliRun* aliRun = runLoader->GetAliRun();
101 if (!aliRun) {
8c2b3fd7 102 AliError("Couldn't get AliRun object");
6d75e4b6 103 return;
104 }
105 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
106 if (!tpc) {
8c2b3fd7 107 AliError("Couldn't get TPC detector");
6d75e4b6 108 return;
109 }
110 AliTPCParam* param = tpc->GetParam();
111 if (!param) {
8c2b3fd7 112 AliError("No TPC parameters available");
6d75e4b6 113 return;
114 }
115
116 //if(fDefaults == 0) SetDefaults();
117
118 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
119 //
120 TParticle *particle; // pointer to a given particle
121 AliTPChit *tpcHit; // pointer to a sigle TPC hit
122 Int_t sector;
123 Int_t ipart;
124 Float_t xyz[5];
125 Float_t pl,pt,tanth,rpad,ratio;
126 Float_t cph,sph;
127
128 //---------------------------------------------------------------
129 // Get the access to the tracks
130 //---------------------------------------------------------------
131
132 TTree *tH = loader->TreeH();
75b27bdb 133 if (tH == 0x0){
8c2b3fd7 134 AliFatal("Can not find TreeH in folder");
75b27bdb 135 return;
136 }
6d75e4b6 137 tpc->SetTreeAddress();
138
139 Stat_t ntracks = tH->GetEntries();
140
141 //Switch to the output file
142
143 if (loader->TreeR() == 0x0) loader->MakeTree("R");
144
8c2b3fd7 145 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
6d75e4b6 146
147 runLoader->CdGAFile();
148 //param->Write(param->GetTitle());
149
150 AliTPCClustersArray carray;
151 carray.Setup(param);
152 carray.SetClusterType("AliTPCcluster");
153 carray.MakeTree(loader->TreeR());
154
155 Int_t nclusters=0; //cluster counter
156
157 //------------------------------------------------------------
158 // Loop over all sectors (72 sectors for 20 deg
159 // segmentation for both lower and upper sectors)
160 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
161 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
162 //
163 // First cluster for sector 0 starts at "0"
164 //------------------------------------------------------------
165
166 for(Int_t isec=0;isec<param->GetNSector();isec++){
167 //MI change
168 param->AdjustCosSin(isec,cph,sph);
169
170 //------------------------------------------------------------
171 // Loop over tracks
172 //------------------------------------------------------------
173
174 for(Int_t track=0;track<ntracks;track++){
175 tpc->ResetHits();
176 tH->GetEvent(track);
177 //
178 // Get number of the TPC hits
179 //
180 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
181
182 // Loop over hits
183 //
184 while(tpcHit){
8ebf1d48 185
6d75e4b6 186 if (tpcHit->fQ == 0.) {
187 tpcHit = (AliTPChit*) tpc->NextHit();
188 continue; //information about track (I.Belikov)
189 }
8ebf1d48 190 sector=tpcHit->fSector; // sector number
191
192 if(sector != isec){
193 tpcHit = (AliTPChit*) tpc->NextHit();
194 continue;
195 }
196 ipart=tpcHit->Track();
197 particle=aliRun->GetMCApp()->Particle(ipart);
198 pl=particle->Pz();
199 pt=particle->Pt();
200 if(pt < 1.e-9) pt=1.e-9;
201 tanth=pl/pt;
202 tanth = TMath::Abs(tanth);
203 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
204 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
205
206 // space-point resolutions
207
6d75e4b6 208 sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
209 sigmaZ =AliTPCcluster::SigmaZ2(rpad,tanth );
210
211 // cluster widths
212
213 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
214 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
215
216 // temporary protection
217
218 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
219 if(sigmaZ < 0.) sigmaZ=0.4e-3;
220 if(clRphi < 0.) clRphi=2.5e-3;
221 if(clZ < 0.) clZ=2.5e-5;
222
223 //
224
225 //
226 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
227 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
228 //
229 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
230 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
231 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
232 Float_t alpha=(isec < param->GetNInnerSector()) ?
233 param->GetInnerAngle() : param->GetOuterAngle();
234 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
235 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
236 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
a1ec4d07 237 if (TMath::Abs(xyz[1])>param->GetZLength(isec)) xyz[1]=tpcHit->Z();
6d75e4b6 238 xyz[2]=sigmaRphi; // fSigmaY2
239 xyz[3]=sigmaZ; // fSigmaZ2
240 xyz[4]=tpcHit->fQ; // q
241
242 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
243 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
244
245 Int_t tracks[3]={tpcHit->Track(), -1, -1};
246 AliTPCcluster cluster(tracks,xyz);
247
248 clrow->InsertCluster(&cluster); nclusters++;
249
250 tpcHit = (AliTPChit*)tpc->NextHit();
251
252
253 } // end of loop over hits
254
255 } // end of loop over tracks
256
257 Int_t nrows=param->GetNRow(isec);
258 for (Int_t irow=0; irow<nrows; irow++) {
259 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
260 if (!clrow) continue;
261 carray.StoreRow(isec,irow);
262 carray.ClearRow(isec,irow);
263 }
264
265 } // end of loop over sectors
266
267 // cerr<<"Number of made clusters : "<<nclusters<<" \n";
268 loader->WriteRecPoints("OVERWRITE");
269
270
271} // end of function
272
8ebf1d48 273
274
275//_________________________________________________________________
276
277
278
279void AliTPCFast::Hits2ExactClusters(AliRunLoader* runLoader) const{
280
281
282
283 AliLoader* loader = runLoader->GetLoader("TPCLoader");
284 if (!loader) {
285 AliError("No TPC loader found");
286 return;
287 }
288 AliTPC* tpc = (AliTPC*) gAlice->GetDetector("TPC");
289 if (!tpc) {
290 AliError("Couldn't get TPC detector");
291 return;
292 }
293 AliTPCParam* param = tpc->GetParam();
294 if (!param) {
295 AliError("No TPC parameters available");
296 return;
297 }
298 //---------------------------------------------------------------
299 // Get the access to the tracks
300 //---------------------------------------------------------------
301
302 TTree *tH = loader->TreeH();
303 if (tH == 0x0) { runLoader->LoadHits("TPC"); loader->TreeH();}
304
305 tpc->SetTreeAddress();
306
307
308 //Switch to the output file
309
310 if (loader->TreeR() == 0x0) loader->MakeTree("R");
311
312 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
313
314 runLoader->CdGAFile();
315 //param->Write(param->GetTitle());
316
317 AliTPCClustersArray carray;
318 carray.Setup(param);
319 carray.SetClusterType("AliTPCclusterMI");
320 carray.MakeTree(loader->TreeR());
321
322 //------------------------------------------------------------
323 // Loop over all sectors (72 sectors for 20 deg)
324 //------------------------------------------------------------
325
326 for(Int_t isec=0;isec<param->GetNSector();isec++){
327 Hits2ExactClustersSector(runLoader, &carray, isec);
328 } // end of loop over sectors
329 loader->WriteRecPoints("OVERWRITE");
330}
331
332
333
334
6d75e4b6 335//_________________________________________________________________
336void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
337 AliTPCClustersArray* clustersArray,
338 Int_t isec) const
339{
340 //--------------------------------------------------------
341 //calculate exact cross point of track and given pad row
8ebf1d48 342 //resulting values are expressed in "local" coordinata
343 //the sigmaY2 and sigmaZ2 of cluster are the shape of cluster parameters
344 // - thay are used later on for error parameterization
6d75e4b6 345 //--------------------------------------------------------
346
347 //-----------------------------------------------------------------
348 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
349 //-----------------------------------------------------------------
350 //
351 if (clustersArray==0){
352 return;
353 }
354 AliLoader* loader = runLoader->GetLoader("TPCLoader");
355 if (!loader) {
8c2b3fd7 356 AliError("No TPC loader found");
6d75e4b6 357 return;
358 }
359 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
360 AliRun* aliRun = runLoader->GetAliRun();
361 if (!aliRun) {
8c2b3fd7 362 AliError("Couldn't get AliRun object");
6d75e4b6 363 return;
364 }
365 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
366 if (!tpc) {
8c2b3fd7 367 AliError("Couldn't get TPC detector");
6d75e4b6 368 return;
369 }
370 AliTPCParam* param = tpc->GetParam();
371 if (!param) {
8c2b3fd7 372 AliError("No TPC parameters available");
6d75e4b6 373 return;
374 }
375 //
8ebf1d48 376 //
6d75e4b6 377 AliTPChit *tpcHit; // pointer to a sigle TPC hit
378 // Int_t sector,nhits;
379 Int_t ipart;
380 const Int_t kcmaxhits=30000;
8c2b3fd7 381 TVector * xxxx = new TVector(kcmaxhits*4);
382 TVector & xxx = *xxxx;
6d75e4b6 383 Int_t maxhits = kcmaxhits;
384 //construct array for each padrow
385 for (Int_t i=0; i<param->GetNRow(isec);i++)
386 clustersArray->CreateRow(isec,i);
387
388 //---------------------------------------------------------------
389 // Get the access to the tracks
390 //---------------------------------------------------------------
391
392 TTree *tH = loader->TreeH();
75b27bdb 393 if (tH == 0x0){
8c2b3fd7 394 AliFatal("Can not find TreeH in folder");
75b27bdb 395 return;
396 }
6d75e4b6 397 tpc->SetTreeAddress();
398
399 Stat_t ntracks = tH->GetEntries();
6d75e4b6 400 //MI change
401 TBranch * branch=0;
402 if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
403 else branch = tH->GetBranch("TPC");
404
405 //------------------------------------------------------------
406 // Loop over tracks
407 //------------------------------------------------------------
408
409 for(Int_t track=0;track<ntracks;track++){
410 Bool_t isInSector=kTRUE;
411 tpc->ResetHits();
8ebf1d48 412 isInSector = tpc->TrackInVolume(isec,track);
6d75e4b6 413 if (!isInSector) continue;
414 //MI change
415 branch->GetEntry(track); // get next track
416 //
6d75e4b6 417 // Loop over hits
418 //
419 Int_t currentIndex=0;
420 Int_t lastrow=-1; //last writen row
421
422 //M.I. changes
423
424 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
425 while(tpcHit){
426
427 Int_t sector=tpcHit->fSector; // sector number
428 if(sector != isec){
429 tpcHit = (AliTPChit*) tpc->NextHit();
430 continue;
431 }
432
433 ipart=tpcHit->Track();
6d75e4b6 434
435 //find row number
436
437 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
438 Int_t index[3]={1,isec,0};
8ebf1d48 439 Int_t currentrow = param->GetPadRow(x,index);
440
6d75e4b6 441 if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
442 if (lastrow<0) lastrow=currentrow;
8ebf1d48 443 //
444 //
445 if (currentrow!=lastrow){
6d75e4b6 446 if (currentIndex>2){
447 Float_t sumx=0;
448 Float_t sumx2=0;
449 Float_t sumx3=0;
450 Float_t sumx4=0;
451 Float_t sumy=0;
452 Float_t sumxy=0;
453 Float_t sumx2y=0;
454 Float_t sumz=0;
455 Float_t sumxz=0;
456 Float_t sumx2z=0;
457 Float_t sumq=0;
d23650ed 458 for (Int_t cindex=0;cindex<currentIndex;cindex++){
459 Float_t x1,x2,x3,x4;
460 x1=x2=x3=x4=xxx(cindex*4);
461 x2*=x1;
6d75e4b6 462 x3*=x2;
463 x4*=x3;
d23650ed 464 sumx+=x1;
6d75e4b6 465 sumx2+=x2;
466 sumx3+=x3;
467 sumx4+=x4;
d23650ed 468 sumy+=xxx(cindex*4+1);
469 sumxy+=xxx(cindex*4+1)*x1;
470 sumx2y+=xxx(cindex*4+1)*x2;
471 sumz+=xxx(cindex*4+2);
472 sumxz+=xxx(cindex*4+2)*x1;
473 sumx2z+=xxx(cindex*4+2)*x2;
474 sumq+=xxx(cindex*4+3);
6d75e4b6 475 }
6d75e4b6 476 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
477 sumx2*(sumx*sumx3-sumx2*sumx2);
478
479 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
480 sumx2*(sumxy*sumx3-sumx2y*sumx2);
481 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
482 sumx2*(sumxz*sumx3-sumx2z*sumx2);
483
484 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
485 sumx2*(sumx*sumx2y-sumx2*sumxy);
486 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
487 sumx2*(sumx*sumx2z-sumx2*sumxz);
488
8ebf1d48 489 if (TMath::Abs(det)<0.00000000000000001){
6d75e4b6 490 tpcHit = (AliTPChit*)tpc->NextHit();
491 continue;
492 }
493
8ebf1d48 494 Float_t y=detay/det;
495 Float_t z=detaz/det;
6d75e4b6 496 Float_t by=detby/det; //y angle
497 Float_t bz=detbz/det; //z angle
498 sumy/=Float_t(currentIndex);
499 sumz/=Float_t(currentIndex);
500
6d75e4b6 501
8ebf1d48 502 //
503 //
504 Float_t sign = (tpcHit->Z()>0)? 1.:-1.;
505 y = (y+0.5)*param->GetPadPitchWidth(isec);
506 z = z*param->GetZWidth();
507 //
508 // y expected shape
509 Double_t sigmay2 = z*fParam->GetDiffL();
510 sigmay2 *= sigmay2;
511 Float_t angulary = by*param->GetPadPitchLength(isec,lastrow);
512 angulary*=angulary;
513 angulary/=12;
514 sigmay2 +=angulary+0.25*param->GetPadPitchWidth(isec)*param->GetPadPitchWidth(isec);
515 //
516 // z expected shape
517 Double_t sigmaz2 = z*fParam->GetDiffT();
518 sigmaz2 *= sigmaz2;
519 Float_t angularz = bz*param->GetPadPitchLength(isec,lastrow);
520 angularz*=angularz;
521 angularz/=12;
522 sigmaz2 +=angularz+0.25;
523 //
524 sigmaz2 = TMath::Min(sigmaz2,1.);
525 sigmay2 = TMath::Min(sigmay2,1.);
526 //
527 //
a1ec4d07 528 z = sign*(param->GetZLength(isec) - z);
529 if (TMath::Abs(z)< param->GetZLength(isec)-1){
8ebf1d48 530 AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
531 if (row!=0) {
532 AliTPCclusterMI* cl = new((AliTPCclusterMI*)row->Append()) AliTPCclusterMI ;
533 // AliTPCclusterMI cl;
534 cl->SetZ(z);
535 cl->SetY(y);
536 cl->SetQ(sumq);
537 if (TMath::Abs(sumq)<1000){
538 cl->SetMax(Short_t(sumq));
539 }
540 else{
541 cl->SetMax(0);
542 }
543 cl->SetSigmaZ2(sigmaz2);
544 cl->SetSigmaY2(sigmay2);
545 cl->SetLabel(ipart,0);
546 cl->SetLabel(-1,1);
547 cl->SetLabel(-1,2);
548 cl->SetType(0);
549 }
550 }
551 } //end of calculating cluster for given row
552 currentIndex=0;
553 lastrow=currentrow;
554 } // end of crossing rows
555 //
6d75e4b6 556 if ( currentIndex>=maxhits){
557 maxhits+=kcmaxhits;
558 xxx.ResizeTo(4*maxhits);
559 }
560 xxx(currentIndex*4)=x[0];
561 xxx(currentIndex*4+1)=x[1];
562 xxx(currentIndex*4+2)=x[2];
563 xxx(currentIndex*4+3)=tpcHit->fQ;
564 currentIndex++;
8ebf1d48 565 tpcHit = (AliTPChit*)tpc->NextHit();
566 } // end of loop over hits
567 } // end of loop over tracks
568 //write padrows to tree
569 for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
570 clustersArray->StoreRow(isec,ii);
571 clustersArray->ClearRow(isec,ii);
572 }
6d75e4b6 573 xxxx->Delete();
8ebf1d48 574
6d75e4b6 575}
576
8ebf1d48 577