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