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