git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New code from Piergiorgio added
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
CommitLineData
4c039060 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$Log$
73042f01 18Revision 1.19.2.4 2000/06/26 07:39:42 kowal2
19Changes to obey the coding rules
20
21Revision 1.19.2.3 2000/06/25 08:38:41 kowal2
22Splitted from AliTPCtracking
23
24Revision 1.19.2.2 2000/06/16 12:59:28 kowal2
25Changed parameter settings
26
27Revision 1.19.2.1 2000/06/09 07:15:07 kowal2
28
29Defaults loaded automatically (hard-wired)
30Optional parameters can be set via macro called in the constructor
31
32Revision 1.19 2000/04/18 19:00:59 fca
33Small bug fixes to TPC files
34
4d68a14a 35Revision 1.18 2000/04/17 09:37:33 kowal2
36removed obsolete AliTPCDigitsDisplay.C
37
cc80f89e 38Revision 1.17.2.2 2000/04/10 08:15:12 kowal2
39
40New, experimental data structure from M. Ivanov
41New tracking algorithm
42Different pad geometry for different sectors
43Digitization rewritten
44
45Revision 1.17.2.1 2000/04/10 07:56:53 kowal2
46Not used anymore - removed
47
48Revision 1.17 2000/01/19 17:17:30 fca
49Introducing a list of lists of hits -- more hits allowed for detector now
50
1cedd08a 51Revision 1.16 1999/11/05 09:29:23 fca
52Accept only signals > 0
53
921bf71a 54Revision 1.15 1999/10/08 06:26:53 fca
55Removed ClustersIndex - not used anymore
56
e674b993 57Revision 1.14 1999/09/29 09:24:33 fca
58Introduction of the Copyright and cvs Log
59
4c039060 60*/
61
fe4da5cc 62///////////////////////////////////////////////////////////////////////////////
63// //
64// Time Projection Chamber //
65// This class contains the basic functions for the Time Projection Chamber //
66// detector. Functions specific to one particular geometry are //
67// contained in the derived classes //
68// //
69//Begin_Html
70/*
1439f98e 71<img src="picts/AliTPCClass.gif">
fe4da5cc 72*/
73//End_Html
74// //
8c555625 75// //
fe4da5cc 76///////////////////////////////////////////////////////////////////////////////
77
73042f01 78//
79
fe4da5cc 80#include <TMath.h>
81#include <TRandom.h>
82#include <TVector.h>
8c555625 83#include <TMatrix.h>
fe4da5cc 84#include <TGeometry.h>
85#include <TNode.h>
86#include <TTUBS.h>
87#include <TObjectTable.h>
1578254f 88#include "TParticle.h"
fe4da5cc 89#include "AliTPC.h"
73042f01 90#include <TFile.h>
fe4da5cc 91#include "AliRun.h"
92#include <iostream.h>
93#include <fstream.h>
94#include "AliMC.h"
95
cc80f89e 96
73042f01 97#include "AliTPCParamSR.h"
8c555625 98#include "AliTPCPRF2D.h"
99#include "AliTPCRF1D.h"
cc80f89e 100#include "AliDigits.h"
101#include "AliSimDigits.h"
102
103#include "AliTPCDigitsArray.h"
104#include "AliCluster.h"
105#include "AliClusters.h"
106#include "AliTPCClustersRow.h"
107#include "AliTPCClustersArray.h"
8c555625 108
73042f01 109#include "AliTPCcluster.h"
110#include "AliTPCclusterer.h"
111#include "AliTPCtracker.h"
8c555625 112
73042f01 113#include <TInterpreter.h>
1283eee5 114
fe4da5cc 115ClassImp(AliTPC)
116
117//_____________________________________________________________________________
118AliTPC::AliTPC()
119{
120 //
121 // Default constructor
122 //
123 fIshunt = 0;
fe4da5cc 124 fHits = 0;
125 fDigits = 0;
fe4da5cc 126 fNsectors = 0;
cc80f89e 127 //MI changes
128 fDigitsArray = 0;
129 fClustersArray = 0;
73042f01 130 fTPCParam=0;
fe4da5cc 131}
132
133//_____________________________________________________________________________
134AliTPC::AliTPC(const char *name, const char *title)
135 : AliDetector(name,title)
136{
137 //
138 // Standard constructor
139 //
140
141 //
142 // Initialise arrays of hits and digits
143 fHits = new TClonesArray("AliTPChit", 176);
4d68a14a 144 gAlice->AddHitList(fHits);
cc80f89e 145 //MI change
146 fDigitsArray = 0;
147 fClustersArray= 0;
fe4da5cc 148 //
149 // Initialise counters
cc80f89e 150 fNsectors = 0;
cc80f89e 151
fe4da5cc 152 //
153 fIshunt = 0;
154 //
155 // Initialise color attributes
156 SetMarkerColor(kYellow);
73042f01 157
158 //
159 // Set TPC parameters
160 //
161
162 if (!strcmp(title,"Default")) {
163 fTPCParam = new AliTPCParamSR;
164 } else {
165 cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
166 fTPCParam=0;
167 }
168
fe4da5cc 169}
170
171//_____________________________________________________________________________
172AliTPC::~AliTPC()
173{
174 //
175 // TPC destructor
176 //
73042f01 177
fe4da5cc 178 fIshunt = 0;
179 delete fHits;
180 delete fDigits;
73042f01 181 delete fTPCParam;
fe4da5cc 182}
183
fe4da5cc 184//_____________________________________________________________________________
185void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
186{
187 //
188 // Add a hit to the list
189 //
190 TClonesArray &lhits = *fHits;
191 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
192}
193
fe4da5cc 194//_____________________________________________________________________________
195void AliTPC::BuildGeometry()
196{
cc80f89e 197
fe4da5cc 198 //
199 // Build TPC ROOT TNode geometry for the event display
200 //
73042f01 201 TNode *nNode, *nTop;
fe4da5cc 202 TTUBS *tubs;
203 Int_t i;
204 const int kColorTPC=19;
1283eee5 205 char name[5], title[25];
fe4da5cc 206 const Double_t kDegrad=TMath::Pi()/180;
1283eee5 207 const Double_t kRaddeg=180./TMath::Pi();
208
1283eee5 209
73042f01 210 Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
211 Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
1283eee5 212
73042f01 213 Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
214 Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
1283eee5 215
216 Int_t nLo = fTPCParam->GetNInnerSector()/2;
217 Int_t nHi = fTPCParam->GetNOuterSector()/2;
218
73042f01 219 const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
220 const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
221 const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
222 const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
1283eee5 223
224
73042f01 225 const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
226 const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
1283eee5 227
228 Double_t rl,ru;
229
230
fe4da5cc 231 //
232 // Get ALICE top node
fe4da5cc 233 //
1283eee5 234
73042f01 235 nTop=gAlice->GetGeometry()->GetNode("alice");
1283eee5 236
237 // inner sectors
238
cc80f89e 239 rl = fTPCParam->GetInnerRadiusLow();
240 ru = fTPCParam->GetInnerRadiusUp();
1283eee5 241
242
fe4da5cc 243 for(i=0;i<nLo;i++) {
244 sprintf(name,"LS%2.2d",i);
1283eee5 245 name[4]='\0';
246 sprintf(title,"TPC low sector %3d",i);
247 title[24]='\0';
248
73042f01 249 tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
250 kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
fe4da5cc 251 tubs->SetNumberOfDivisions(1);
73042f01 252 nTop->cd();
253 nNode = new TNode(name,title,name,0,0,0,"");
254 nNode->SetLineColor(kColorTPC);
255 fNodes->Add(nNode);
fe4da5cc 256 }
1283eee5 257
fe4da5cc 258 // Outer sectors
1283eee5 259
cc80f89e 260 rl = fTPCParam->GetOuterRadiusLow();
261 ru = fTPCParam->GetOuterRadiusUp();
1283eee5 262
fe4da5cc 263 for(i=0;i<nHi;i++) {
264 sprintf(name,"US%2.2d",i);
1283eee5 265 name[4]='\0';
fe4da5cc 266 sprintf(title,"TPC upper sector %d",i);
1283eee5 267 title[24]='\0';
73042f01 268 tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
269 khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
fe4da5cc 270 tubs->SetNumberOfDivisions(1);
73042f01 271 nTop->cd();
272 nNode = new TNode(name,title,name,0,0,0,"");
273 nNode->SetLineColor(kColorTPC);
274 fNodes->Add(nNode);
fe4da5cc 275 }
cc80f89e 276
73042f01 277}
1283eee5 278
fe4da5cc 279//_____________________________________________________________________________
280Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
281{
282 //
283 // Calculate distance from TPC to mouse on the display
284 // Dummy procedure
285 //
286 return 9999;
287}
288
73042f01 289void AliTPC::Clusters2Tracks(TFile *of) {
3c0f9266 290 //-----------------------------------------------------------------
291 // This is a track finder.
3c0f9266 292 //-----------------------------------------------------------------
73042f01 293 AliTPCtracker::Clusters2Tracks(fTPCParam,of);
fe4da5cc 294}
295
296//_____________________________________________________________________________
297void AliTPC::CreateMaterials()
298{
8c555625 299 //-----------------------------------------------
37831078 300 // Create Materials for for TPC simulations
8c555625 301 //-----------------------------------------------
302
303 //-----------------------------------------------------------------
304 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
305 //-----------------------------------------------------------------
1283eee5 306
73042f01 307 Int_t iSXFLD=gAlice->Field()->Integ();
308 Float_t sXMGMX=gAlice->Field()->Max();
1283eee5 309
310 Float_t amat[5]; // atomic numbers
311 Float_t zmat[5]; // z
312 Float_t wmat[5]; // proportions
313
314 Float_t density;
37831078 315 Float_t apure[2];
1283eee5 316
1283eee5 317
37831078 318 //***************** Gases *************************
319
320 //-------------------------------------------------
1283eee5 321 // pure gases
37831078 322 //-------------------------------------------------
1283eee5 323
37831078 324 // Neon
1283eee5 325
326
37831078 327 amat[0]= 20.18;
328 zmat[0]= 10.;
1283eee5 329 density = 0.0009;
37831078 330
331 apure[0]=amat[0];
1283eee5 332
37831078 333 AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
1283eee5 334
37831078 335 // Argon
1283eee5 336
37831078 337 amat[0]= 39.948;
338 zmat[0]= 18.;
339 density = 0.001782;
1283eee5 340
37831078 341 apure[1]=amat[0];
1283eee5 342
37831078 343 AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
344
1283eee5 345
346 //--------------------------------------------------------------
347 // gases - compounds
348 //--------------------------------------------------------------
349
350 Float_t amol[3];
351
37831078 352 // CO2
1283eee5 353
354 amat[0]=12.011;
355 amat[1]=15.9994;
356
357 zmat[0]=6.;
358 zmat[1]=8.;
359
360 wmat[0]=1.;
361 wmat[1]=2.;
362
363 density=0.001977;
364
365 amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
366
367 AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
37831078 368
1283eee5 369 // CF4
370
371 amat[0]=12.011;
372 amat[1]=18.998;
373
374 zmat[0]=6.;
375 zmat[1]=9.;
376
377 wmat[0]=1.;
378 wmat[1]=4.;
379
380 density=0.003034;
381
382 amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
383
384 AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
385
37831078 386
1283eee5 387 // CH4
388
389 amat[0]=12.011;
390 amat[1]=1.;
391
392 zmat[0]=6.;
393 zmat[1]=1.;
394
395 wmat[0]=1.;
396 wmat[1]=4.;
397
398 density=0.000717;
399
400 amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
401
402 AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
403
404 //----------------------------------------------------------------
405 // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
406 //----------------------------------------------------------------
37831078 407
408 char namate[21];
1283eee5 409 density = 0.;
410 Float_t am=0;
411 Int_t nc;
37831078 412 Float_t rho,absl,X0,buf[1];
1283eee5 413 Int_t nbuf;
37831078 414 Float_t a,z;
1283eee5 415
416 for(nc = 0;nc<fNoComp;nc++)
417 {
418
419 // retrive material constants
420
37831078 421 gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
1283eee5 422
423 amat[nc] = a;
424 zmat[nc] = z;
425
426 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
427
37831078 428 am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
1283eee5 429 density += fMixtProp[nc]*rho; // density of the mixture
430
431 }
432
433 // mixture proportions by weight!
434
435 for(nc = 0;nc<fNoComp;nc++)
436 {
437
438 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
439
37831078 440 wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
441 apure[nnc] : amol[nnc])/am;
442
443 }
444
445 // Drift gases 1 - nonsensitive, 2 - sensitive
1283eee5 446
1283eee5 447 AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
448 AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
1283eee5 449
1283eee5 450
37831078 451 // Air
452
453 amat[0] = 14.61;
454 zmat[0] = 7.3;
455 density = 0.001205;
1283eee5 456
37831078 457 AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
1283eee5 458
1283eee5 459
460 //----------------------------------------------------------------------
461 // solid materials
462 //----------------------------------------------------------------------
463
1283eee5 464
37831078 465 // Kevlar C14H22O2N2
1283eee5 466
37831078 467 amat[0] = 12.011;
468 amat[1] = 1.;
469 amat[2] = 15.999;
470 amat[3] = 14.006;
1283eee5 471
37831078 472 zmat[0] = 6.;
473 zmat[1] = 1.;
474 zmat[2] = 8.;
475 zmat[3] = 7.;
476
477 wmat[0] = 14.;
478 wmat[1] = 22.;
479 wmat[2] = 2.;
480 wmat[3] = 2.;
481
482 density = 1.45;
483
484 AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
485
486 // NOMEX
1283eee5 487
37831078 488 amat[0] = 12.011;
489 amat[1] = 1.;
490 amat[2] = 15.999;
491 amat[3] = 14.006;
492
493 zmat[0] = 6.;
494 zmat[1] = 1.;
495 zmat[2] = 8.;
496 zmat[3] = 7.;
497
498 wmat[0] = 14.;
499 wmat[1] = 22.;
500 wmat[2] = 2.;
501 wmat[3] = 2.;
502
503 density = 0.03;
1283eee5 504
fe4da5cc 505
37831078 506 AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
507
508 // Makrolon C16H18O3
509
510 amat[0] = 12.011;
511 amat[1] = 1.;
512 amat[2] = 15.999;
1283eee5 513
37831078 514 zmat[0] = 6.;
515 zmat[1] = 1.;
516 zmat[2] = 8.;
517
518 wmat[0] = 16.;
519 wmat[1] = 18.;
520 wmat[2] = 3.;
521
522 density = 1.2;
523
524 AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
525
1283eee5 526 // Mylar C5H4O2
527
528 amat[0]=12.011;
529 amat[1]=1.;
530 amat[2]=15.9994;
531
532 zmat[0]=6.;
533 zmat[1]=1.;
534 zmat[2]=8.;
535
536 wmat[0]=5.;
537 wmat[1]=4.;
538 wmat[2]=2.;
539
540 density = 1.39;
fe4da5cc 541
37831078 542 AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
1283eee5 543
37831078 544 // G10 60% SiO2 + 40% epoxy, I use A and Z for SiO2
1283eee5 545
37831078 546 amat[0]=28.086;
547 amat[1]=15.9994;
1283eee5 548
37831078 549 zmat[0]=14.;
550 zmat[1]=8.;
1283eee5 551
37831078 552 wmat[0]=1.;
553 wmat[1]=2.;
1283eee5 554
37831078 555 density = 1.7;
1283eee5 556
37831078 557 AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz
558
559 gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf);
1283eee5 560
37831078 561 AliMaterial(39,"G10",amat[0],zmat[0],density,999.,999.);
1283eee5 562
37831078 563 // Al
1283eee5 564
37831078 565 amat[0] = 26.98;
566 zmat[0] = 13.;
1283eee5 567
37831078 568 density = 2.7;
1283eee5 569
37831078 570 AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
1283eee5 571
37831078 572 // Si
1283eee5 573
37831078 574 amat[0] = 28.086;
575 zmat[0] = 14.,
1283eee5 576
37831078 577 density = 2.33;
1283eee5 578
37831078 579 AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
1283eee5 580
37831078 581 // Cu
1283eee5 582
37831078 583 amat[0] = 63.546;
584 zmat[0] = 29.;
1283eee5 585
37831078 586 density = 8.96;
1283eee5 587
37831078 588 AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
1283eee5 589
37831078 590 // Tedlar C2H3F
1283eee5 591
37831078 592 amat[0] = 12.011;
593 amat[1] = 1.;
594 amat[2] = 18.998;
1283eee5 595
37831078 596 zmat[0] = 6.;
597 zmat[1] = 1.;
598 zmat[2] = 9.;
1283eee5 599
37831078 600 wmat[0] = 2.;
601 wmat[1] = 3.;
602 wmat[2] = 1.;
1283eee5 603
37831078 604 density = 1.71;
1283eee5 605
37831078 606 AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
1283eee5 607
1283eee5 608
37831078 609 // Plexiglas C5H8O2
1283eee5 610
37831078 611 amat[0]=12.011;
612 amat[1]=1.;
613 amat[2]=15.9994;
1283eee5 614
37831078 615 zmat[0]=6.;
616 zmat[1]=1.;
617 zmat[2]=8.;
1283eee5 618
37831078 619 wmat[0]=5.;
620 wmat[1]=8.;
621 wmat[2]=2.;
1283eee5 622
37831078 623 density=1.18;
1283eee5 624
37831078 625 AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
1283eee5 626
1283eee5 627
37831078 628
629 //----------------------------------------------------------
630 // tracking media for gases
631 //----------------------------------------------------------
632
633 AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
634 AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
635 AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
636 AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
637
638 //-----------------------------------------------------------
639 // tracking media for solids
640 //-----------------------------------------------------------
641
642 AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
643 AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
644 AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
645 AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
646 AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
647 AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
648 AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
649 AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
650 AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
651 AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
1283eee5 652
fe4da5cc 653}
654
fe4da5cc 655
73042f01 656void AliTPC::Digits2Clusters(TFile *of)
fe4da5cc 657{
3c0f9266 658 //-----------------------------------------------------------------
659 // This is a simple cluster finder.
3c0f9266 660 //-----------------------------------------------------------------
73042f01 661 AliTPCclusterer::Digits2Clusters(fTPCParam,of);
fe4da5cc 662}
663
73042f01 664extern Double_t SigmaY2(Double_t, Double_t, Double_t);
665extern Double_t SigmaZ2(Double_t, Double_t);
fe4da5cc 666//_____________________________________________________________________________
73042f01 667void AliTPC::Hits2Clusters(TFile *of)
fe4da5cc 668{
8c555625 669 //--------------------------------------------------------
fe4da5cc 670 // TPC simple cluster generator from hits
671 // obtained from the TPC Fast Simulator
8c555625 672 // The point errors are taken from the parametrization
673 //--------------------------------------------------------
674
675 //-----------------------------------------------------------------
676 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
677 //-----------------------------------------------------------------
73042f01 678 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
679 // Jouri.Belikov@cern.ch
680 //----------------------------------------------------------------
681
682 /////////////////////////////////////////////////////////////////////////////
683 //
684 //---------------------------------------------------------------------
685 // ALICE TPC Cluster Parameters
686 //--------------------------------------------------------------------
687
688
689
690 // Cluster width in rphi
691 const Float_t kACrphi=0.18322;
692 const Float_t kBCrphi=0.59551e-3;
693 const Float_t kCCrphi=0.60952e-1;
694 // Cluster width in z
695 const Float_t kACz=0.19081;
696 const Float_t kBCz=0.55938e-3;
697 const Float_t kCCz=0.30428;
698
699 TDirectory *savedir=gDirectory;
700
701 if (!of->IsOpen()) {
702 cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
703 return;
704 }
3c0f9266 705
cc80f89e 706 if(fTPCParam == 0){
707 printf("AliTPCParam MUST be created firstly\n");
708 return;
709 }
710
73042f01 711 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
fe4da5cc 712 //
1578254f 713 TParticle *particle; // pointer to a given particle
fe4da5cc 714 AliTPChit *tpcHit; // pointer to a sigle TPC hit
73042f01 715 TClonesArray *particles; //pointer to the particle list
fe4da5cc 716 Int_t sector,nhits;
717 Int_t ipart;
718 Float_t xyz[5];
719 Float_t pl,pt,tanth,rpad,ratio;
fe4da5cc 720 Float_t cph,sph;
721
722 //---------------------------------------------------------------
723 // Get the access to the tracks
724 //---------------------------------------------------------------
725
73042f01 726 TTree *tH = gAlice->TreeH();
727 Stat_t ntracks = tH->GetEntries();
728 particles=gAlice->Particles();
729
730 //Switch to the output file
731 of->cd();
732
733 fTPCParam->Write(fTPCParam->GetTitle());
734 AliTPCClustersArray carray;
735 carray.Setup(fTPCParam);
736 carray.SetClusterType("AliTPCcluster");
737 carray.MakeTree();
738
739 Int_t nclusters=0; //cluster counter
fe4da5cc 740
741 //------------------------------------------------------------
cc80f89e 742 // Loop over all sectors (72 sectors for 20 deg
743 // segmentation for both lower and upper sectors)
744 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
3c0f9266 745 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
fe4da5cc 746 //
3c0f9266 747 // First cluster for sector 0 starts at "0"
fe4da5cc 748 //------------------------------------------------------------
cc80f89e 749
750 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
8c555625 751 //MI change
cc80f89e 752 fTPCParam->AdjustCosSin(isec,cph,sph);
fe4da5cc 753
754 //------------------------------------------------------------
755 // Loop over tracks
756 //------------------------------------------------------------
757
758 for(Int_t track=0;track<ntracks;track++){
759 ResetHits();
73042f01 760 tH->GetEvent(track);
fe4da5cc 761 //
73042f01 762 // Get number of the TPC hits
fe4da5cc 763 //
764 nhits=fHits->GetEntriesFast();
fe4da5cc 765 //
766 // Loop over hits
767 //
768 for(Int_t hit=0;hit<nhits;hit++){
769 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
cc80f89e 770 if (tpcHit->fQ == 0.) continue; //information about track (I.Belikov)
fe4da5cc 771 sector=tpcHit->fSector; // sector number
772 if(sector != isec) continue; //terminate iteration
773 ipart=tpcHit->fTrack;
73042f01 774 particle=(TParticle*)particles->UncheckedAt(ipart);
1578254f 775 pl=particle->Pz();
776 pt=particle->Pt();
fe4da5cc 777 if(pt < 1.e-9) pt=1.e-9;
778 tanth=pl/pt;
779 tanth = TMath::Abs(tanth);
780 rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
781 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
73042f01 782
fe4da5cc 783 // space-point resolutions
784
73042f01 785 sigmaRphi=SigmaY2(rpad,tanth,pt);
786 sigmaZ =SigmaZ2(rpad,tanth );
fe4da5cc 787
788 // cluster widths
789
73042f01 790 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
791 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
fe4da5cc 792
793 // temporary protection
794
73042f01 795 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
796 if(sigmaZ < 0.) sigmaZ=0.4e-3;
797 if(clRphi < 0.) clRphi=2.5e-3;
798 if(clZ < 0.) clZ=2.5e-5;
fe4da5cc 799
800 //
cc80f89e 801
802 //
803 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
fe4da5cc 804 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
805 //
73042f01 806 Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
fe4da5cc 807 Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
73042f01 808 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
809 Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
810 fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
811 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
812 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
813 xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigmaZ)); // z
814 if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->fZ;
cc80f89e 815 xyz[2]=tpcHit->fQ; // q
73042f01 816 xyz[3]=sigmaRphi; // fSigmaY2
817 xyz[4]=sigmaZ; // fSigmaZ2
818
819 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
820 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
821
822 Int_t tracks[3]={tpcHit->fTrack, -1, -1};
823 AliTPCcluster cluster(xyz,tracks);
824
825 clrow->InsertCluster(&cluster); nclusters++;
826
fe4da5cc 827 } // end of loop over hits
73042f01 828
829 } // end of loop over tracks
830
831 Int_t nrows=fTPCParam->GetNRow(isec);
832 for (Int_t irow=0; irow<nrows; irow++) {
833 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
834 if (!clrow) continue;
835 carray.StoreRow(isec,irow);
836 carray.ClearRow(isec,irow);
837 }
838
cc80f89e 839 } // end of loop over sectors
73042f01 840
841 cerr<<"Number of made clusters : "<<nclusters<<" \n";
842
843 carray.GetTree()->Write();
844
845 savedir->cd(); //switch back to the input file
cc80f89e 846
847} // end of function
848
849//_________________________________________________________________
850void AliTPC::Hits2ExactClustersSector(Int_t isec)
851{
852 //--------------------------------------------------------
853 //calculate exact cross point of track and given pad row
854 //resulting values are expressed in "digit" coordinata
855 //--------------------------------------------------------
856
857 //-----------------------------------------------------------------
858 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
859 //-----------------------------------------------------------------
860 //
861 if (fClustersArray==0){
862 return;
863 }
864 //
865 TParticle *particle; // pointer to a given particle
866 AliTPChit *tpcHit; // pointer to a sigle TPC hit
73042f01 867 TClonesArray *particles; //pointer to the particle list
cc80f89e 868 Int_t sector,nhits;
869 Int_t ipart;
73042f01 870 const Int_t kcmaxhits=30000;
871 TVector * xxxx = new TVector(kcmaxhits*4);
cc80f89e 872 TVector & xxx = *xxxx;
73042f01 873 Int_t maxhits = kcmaxhits;
cc80f89e 874 //construct array for each padrow
875 for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++)
876 fClustersArray->CreateRow(isec,i);
fe4da5cc 877
cc80f89e 878 //---------------------------------------------------------------
879 // Get the access to the tracks
880 //---------------------------------------------------------------
881
73042f01 882 TTree *tH = gAlice->TreeH();
883 Stat_t ntracks = tH->GetEntries();
884 particles=gAlice->Particles();
885 Int_t npart = particles->GetEntriesFast();
cc80f89e 886
887 //------------------------------------------------------------
888 // Loop over tracks
889 //------------------------------------------------------------
fe4da5cc 890
cc80f89e 891 for(Int_t track=0;track<ntracks;track++){
892 ResetHits();
73042f01 893 tH->GetEvent(track);
cc80f89e 894 //
895 // Get number of the TPC hits and a pointer
896 // to the particles
897 //
898 nhits=fHits->GetEntriesFast();
899 //
900 // Loop over hits
901 //
902 Int_t currentIndex=0;
903 Int_t lastrow=-1; //last writen row
904 for(Int_t hit=0;hit<nhits;hit++){
905 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
906 if (tpcHit==0) continue;
907 sector=tpcHit->fSector; // sector number
908 if(sector != isec) continue;
909 ipart=tpcHit->fTrack;
73042f01 910 if (ipart<npart) particle=(TParticle*)particles->UncheckedAt(ipart);
cc80f89e 911
912 //find row number
913
914 Float_t x[3]={tpcHit->fX,tpcHit->fY,tpcHit->fZ};
915 Int_t index[3]={1,isec,0};
916 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
917 if (currentrow<0) continue;
918 if (lastrow<0) lastrow=currentrow;
919 if (currentrow==lastrow){
920 if ( currentIndex>=maxhits){
73042f01 921 maxhits+=kcmaxhits;
cc80f89e 922 xxx.ResizeTo(4*maxhits);
923 }
924 xxx(currentIndex*4)=x[0];
925 xxx(currentIndex*4+1)=x[1];
926 xxx(currentIndex*4+2)=x[2];
927 xxx(currentIndex*4+3)=tpcHit->fQ;
928 currentIndex++;
929 }
930 else
931 if (currentIndex>2){
932 Float_t sumx=0;
933 Float_t sumx2=0;
934 Float_t sumx3=0;
935 Float_t sumx4=0;
936 Float_t sumy=0;
937 Float_t sumxy=0;
938 Float_t sumx2y=0;
939 Float_t sumz=0;
940 Float_t sumxz=0;
941 Float_t sumx2z=0;
942 Float_t sumq=0;
943 for (Int_t index=0;index<currentIndex;index++){
944 Float_t x,x2,x3,x4;
945 x=x2=x3=x4=xxx(index*4);
946 x2*=x;
947 x3*=x2;
948 x4*=x3;
949 sumx+=x;
950 sumx2+=x2;
951 sumx3+=x3;
952 sumx4+=x4;
953 sumy+=xxx(index*4+1);
954 sumxy+=xxx(index*4+1)*x;
955 sumx2y+=xxx(index*4+1)*x2;
956 sumz+=xxx(index*4+2);
957 sumxz+=xxx(index*4+2)*x;
958 sumx2z+=xxx(index*4+2)*x2;
959 sumq+=xxx(index*4+3);
960 }
73042f01 961 Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
cc80f89e 962 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
963 sumx2*(sumx*sumx3-sumx2*sumx2);
964
965 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
966 sumx2*(sumxy*sumx3-sumx2y*sumx2);
967 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
968 sumx2*(sumxz*sumx3-sumx2z*sumx2);
969
970 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
971 sumx2*(sumx*sumx2y-sumx2*sumxy);
972 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
973 sumx2*(sumx*sumx2z-sumx2*sumxz);
974
73042f01 975 Float_t y=detay/det+centralPad;
cc80f89e 976 Float_t z=detaz/det;
977 Float_t by=detby/det; //y angle
978 Float_t bz=detbz/det; //z angle
979 sumy/=Float_t(currentIndex);
980 sumz/=Float_t(currentIndex);
981 AliCluster cl;
982 cl.fX=z;
983 cl.fY=y;
984 cl.fQ=sumq;
985 cl.fSigmaX2=bz;
986 cl.fSigmaY2=by;
987 cl.fTracks[0]=ipart;
988
989 AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
990 if (row!=0) row->InsertCluster(&cl);
991 currentIndex=0;
992 lastrow=currentrow;
993 } //end of calculating cluster for given row
994
995
996
997 } // end of loop over hits
998 } // end of loop over tracks
999 //write padrows to tree
1000 for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1001 fClustersArray->StoreRow(isec,ii);
1002 fClustersArray->ClearRow(isec,ii);
1003 }
1004 xxxx->Delete();
1005
fe4da5cc 1006}
1007
cc80f89e 1008//__________________________________________________________________
8c555625 1009void AliTPC::Hits2Digits()
1010{
8c555625 1011 //----------------------------------------------------
cc80f89e 1012 // Loop over all sectors
1013 //----------------------------------------------------
1014
1015 if(fTPCParam == 0){
1016 printf("AliTPCParam MUST be created firstly\n");
1017 return;
1018 }
1019
1020 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1021
8c555625 1022}
1023
1024
fe4da5cc 1025//_____________________________________________________________________________
8c555625 1026void AliTPC::Hits2DigitsSector(Int_t isec)
fe4da5cc 1027{
8c555625 1028 //-------------------------------------------------------------------
fe4da5cc 1029 // TPC conversion from hits to digits.
8c555625 1030 //-------------------------------------------------------------------
1031
1032 //-----------------------------------------------------------------
1033 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1034 //-----------------------------------------------------------------
1035
fe4da5cc 1036 //-------------------------------------------------------
8c555625 1037 // Get the access to the track hits
fe4da5cc 1038 //-------------------------------------------------------
8c555625 1039
cc80f89e 1040
73042f01 1041 TTree *tH = gAlice->TreeH(); // pointer to the hits tree
1042 Stat_t ntracks = tH->GetEntries();
8c555625 1043
1044 if( ntracks > 0){
1045
1046 //-------------------------------------------
1047 // Only if there are any tracks...
1048 //-------------------------------------------
1049
8c555625 1050 TObjArray **row;
fe4da5cc 1051
8c555625 1052 printf("*** Processing sector number %d ***\n",isec);
1053
1054 Int_t nrows =fTPCParam->GetNRow(isec);
1055
1056 row= new TObjArray* [nrows];
fe4da5cc 1057
73042f01 1058 MakeSector(isec,nrows,tH,ntracks,row);
8c555625 1059
1060 //--------------------------------------------------------
1061 // Digitize this sector, row by row
1062 // row[i] is the pointer to the TObjArray of TVectors,
1063 // each one containing electrons accepted on this
1064 // row, assigned into tracks
1065 //--------------------------------------------------------
1066
1067 Int_t i;
1068
cc80f89e 1069 if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree();
8c555625 1070
cc80f89e 1071 for (i=0;i<nrows;i++){
8c555625 1072
cc80f89e 1073 AliDigits * dig = fDigitsArray->CreateRow(isec,i);
8c555625 1074
cc80f89e 1075 DigitizeRow(i,isec,row);
8c555625 1076
cc80f89e 1077 fDigitsArray->StoreRow(isec,i);
8c555625 1078
73042f01 1079 Int_t ndig = dig->GetDigitSize();
cc80f89e 1080
1081 printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1082
1083 fDigitsArray->ClearRow(isec,i);
8c555625 1084
cc80f89e 1085
8c555625 1086 } // end of the sector digitization
8c555625 1087
cc80f89e 1088 for(i=0;i<nrows;i++){
1089 row[i]->Delete();
1090 }
1091
8c555625 1092 delete [] row; // delete the array of pointers to TObjArray-s
1093
1094 } // ntracks >0
8c555625 1095
cc80f89e 1096} // end of Hits2DigitsSector
8c555625 1097
8c555625 1098
8c555625 1099//_____________________________________________________________________________
cc80f89e 1100void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
8c555625 1101{
1102 //-----------------------------------------------------------
1103 // Single row digitization, coupling from the neighbouring
1104 // rows taken into account
1105 //-----------------------------------------------------------
1106
1107 //-----------------------------------------------------------------
1108 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
cc80f89e 1109 // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
8c555625 1110 //-----------------------------------------------------------------
1111
1112
8c555625 1113 Float_t zerosup = fTPCParam->GetZeroSup();
1114 Int_t nrows =fTPCParam->GetNRow(isec);
cc80f89e 1115 fCurrentIndex[1]= isec;
8c555625 1116
8c555625 1117
73042f01 1118 Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1119 Int_t nofTbins = fTPCParam->GetMaxTBin();
1120 Int_t indexRange[4];
8c555625 1121 //
1122 // Integrated signal for this row
1123 // and a single track signal
cc80f89e 1124 //
73042f01 1125 TMatrix *m1 = new TMatrix(0,nofPads,0,nofTbins); // integrated
1126 TMatrix *m2 = new TMatrix(0,nofPads,0,nofTbins); // single
8c555625 1127 //
73042f01 1128 TMatrix &total = *m1;
8c555625 1129
1130 // Array of pointers to the label-signal list
1131
73042f01 1132 Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1133 Float_t **pList = new Float_t* [nofDigits];
8c555625 1134
1135 Int_t lp;
cc80f89e 1136 Int_t i1;
73042f01 1137 for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
8c555625 1138 //
cc80f89e 1139 //calculate signal
1140 //
1141 Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1142 Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1143 for (Int_t row= row1;row<=row2;row++){
1144 Int_t nTracks= rows[row]->GetEntries();
1145 for (i1=0;i1<nTracks;i1++){
1146 fCurrentIndex[2]= row;
1147 fCurrentIndex[3]=irow;
1148 if (row==irow){
1149 m2->Zero(); // clear single track signal matrix
73042f01 1150 Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1151 GetList(trackLabel,nofPads,m2,indexRange,pList);
cc80f89e 1152 }
73042f01 1153 else GetSignal(rows[row],i1,0,m1,indexRange);
cc80f89e 1154 }
8c555625 1155 }
cc80f89e 1156
8c555625 1157 Int_t tracks[3];
8c555625 1158
cc80f89e 1159 AliDigits *dig = fDigitsArray->GetRow(isec,irow);
73042f01 1160 for(Int_t ip=0;ip<nofPads;ip++){
1161 for(Int_t it=0;it<nofTbins;it++){
8c555625 1162
73042f01 1163 Float_t q = total(ip,it);
8c555625 1164
73042f01 1165 Int_t gi =it*nofPads+ip; // global index
8c555625 1166
cc80f89e 1167 q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
8c555625 1168
cc80f89e 1169 q = (Int_t)q;
1170
1171 if(q <=zerosup) continue; // do not fill zeros
73042f01 1172 if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation
8c555625 1173
1174 //
1175 // "real" signal or electronic noise (list = -1)?
1176 //
1177
1178 for(Int_t j1=0;j1<3;j1++){
1179 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1180 }
1181
cc80f89e 1182//Begin_Html
1183/*
1184 <A NAME="AliDigits"></A>
1185 using of AliDigits object
1186*/
1187//End_Html
1188 dig->SetDigitFast((Short_t)q,it,ip);
1189 if (fDigitsArray->IsSimulated())
1190 {
1191 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1192 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1193 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1194 }
1195
8c555625 1196
1197 } // end of loop over time buckets
1198 } // end of lop over pads
1199
1200 //
1201 // This row has been digitized, delete nonused stuff
1202 //
1203
73042f01 1204 for(lp=0;lp<nofDigits;lp++){
8c555625 1205 if(pList[lp]) delete [] pList[lp];
1206 }
1207
1208 delete [] pList;
1209
1210 delete m1;
1211 delete m2;
cc80f89e 1212 // delete m3;
8c555625 1213
1214} // end of DigitizeRow
cc80f89e 1215
8c555625 1216//_____________________________________________________________________________
cc80f89e 1217
1218Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
73042f01 1219 Int_t *indexRange)
8c555625 1220{
1221
1222 //---------------------------------------------------------------
1223 // Calculates 2-D signal (pad,time) for a single track,
1224 // returns a pointer to the signal matrix and the track label
1225 // No digitization is performed at this level!!!
1226 //---------------------------------------------------------------
1227
1228 //-----------------------------------------------------------------
1229 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
cc80f89e 1230 // Modified: Marian Ivanov
8c555625 1231 //-----------------------------------------------------------------
1232
1233 TVector *tv;
8c555625 1234
8c555625 1235 tv = (TVector*)p1->At(ntr); // pointer to a track
1236 TVector &v = *tv;
1237
1238 Float_t label = v(0);
73042f01 1239 Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
8c555625 1240
8c555625 1241 Int_t nElectrons = (tv->GetNrows()-1)/4;
73042f01 1242 indexRange[0]=9999; // min pad
1243 indexRange[1]=-1; // max pad
1244 indexRange[2]=9999; //min time
1245 indexRange[3]=-1; // max time
8c555625 1246
cc80f89e 1247 // Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1248
1249 TMatrix &signal = *m1;
1250 TMatrix &total = *m2;
8c555625 1251 //
1252 // Loop over all electrons
1253 //
8c555625 1254 for(Int_t nel=0; nel<nElectrons; nel++){
cc80f89e 1255 Int_t idx=nel*4;
1256 Float_t aval = v(idx+4);
1257 Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1258 Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1259 Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
8c555625 1260
cc80f89e 1261 if (n>0) for (Int_t i =0; i<n; i++){
1262 Int_t *index = fTPCParam->GetResBin(i);
73042f01 1263 Int_t pad=index[1]+centralPad; //in digit coordinates central pad has coordinate 0
cc80f89e 1264 if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
1265 Int_t time=index[2];
1266 Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
1267 weight *= eltoadcfac;
1268
1269 if (m1!=0) signal(pad,time)+=weight;
1270 total(pad,time)+=weight;
73042f01 1271 indexRange[0]=TMath::Min(indexRange[0],pad);
1272 indexRange[1]=TMath::Max(indexRange[1],pad);
1273 indexRange[2]=TMath::Min(indexRange[2],time);
1274 indexRange[3]=TMath::Max(indexRange[3],time);
cc80f89e 1275 }
1276 }
8c555625 1277 } // end of loop over electrons
cc80f89e 1278
8c555625 1279 return label; // returns track label when finished
1280}
1281
1282//_____________________________________________________________________________
73042f01 1283void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
8c555625 1284 Float_t **pList)
1285{
1286 //----------------------------------------------------------------------
1287 // Updates the list of tracks contributing to digits for a given row
1288 //----------------------------------------------------------------------
1289
1290 //-----------------------------------------------------------------
1291 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1292 //-----------------------------------------------------------------
1293
1294 TMatrix &signal = *m;
1295
1296 // lop over nonzero digits
1297
73042f01 1298 for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1299 for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
8c555625 1300
1301
921bf71a 1302 // accept only the contribution larger than 500 electrons (1/2 s_noise)
1303
cc80f89e 1304 if(signal(ip,it)<0.5) continue;
921bf71a 1305
1306
73042f01 1307 Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
8c555625 1308
73042f01 1309 if(!pList[globalIndex]){
8c555625 1310
1311 //
1312 // Create new list (6 elements - 3 signals and 3 labels),
8c555625 1313 //
1314
73042f01 1315 pList[globalIndex] = new Float_t [6];
8c555625 1316
1317 // set list to -1
1318
73042f01 1319 *pList[globalIndex] = -1.;
1320 *(pList[globalIndex]+1) = -1.;
1321 *(pList[globalIndex]+2) = -1.;
1322 *(pList[globalIndex]+3) = -1.;
1323 *(pList[globalIndex]+4) = -1.;
1324 *(pList[globalIndex]+5) = -1.;
8c555625 1325
1326
73042f01 1327 *pList[globalIndex] = label;
1328 *(pList[globalIndex]+3) = signal(ip,it);
8c555625 1329 }
1330 else{
1331
1332 // check the signal magnitude
1333
73042f01 1334 Float_t highest = *(pList[globalIndex]+3);
1335 Float_t middle = *(pList[globalIndex]+4);
1336 Float_t lowest = *(pList[globalIndex]+5);
8c555625 1337
1338 //
1339 // compare the new signal with already existing list
1340 //
1341
1342 if(signal(ip,it)<lowest) continue; // neglect this track
1343
1344 //
1345
1346 if (signal(ip,it)>highest){
73042f01 1347 *(pList[globalIndex]+5) = middle;
1348 *(pList[globalIndex]+4) = highest;
1349 *(pList[globalIndex]+3) = signal(ip,it);
8c555625 1350
73042f01 1351 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1352 *(pList[globalIndex]+1) = *pList[globalIndex];
1353 *pList[globalIndex] = label;
8c555625 1354 }
1355 else if (signal(ip,it)>middle){
73042f01 1356 *(pList[globalIndex]+5) = middle;
1357 *(pList[globalIndex]+4) = signal(ip,it);
8c555625 1358
73042f01 1359 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1360 *(pList[globalIndex]+1) = label;
8c555625 1361 }
1362 else{
73042f01 1363 *(pList[globalIndex]+5) = signal(ip,it);
1364 *(pList[globalIndex]+2) = label;
8c555625 1365 }
1366 }
1367
1368 } // end of loop over pads
1369 } // end of loop over time bins
1370
1371
1372
8c555625 1373}//end of GetList
1374//___________________________________________________________________
1375void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1376 Stat_t ntracks,TObjArray **row)
1377{
1378
1379 //-----------------------------------------------------------------
1380 // Prepares the sector digitization, creates the vectors of
1381 // tracks for each row of this sector. The track vector
1382 // contains the track label and the position of electrons.
1383 //-----------------------------------------------------------------
1384
1385 //-----------------------------------------------------------------
1386 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1387 //-----------------------------------------------------------------
1388
cc80f89e 1389 Float_t gasgain = fTPCParam->GetGasGain();
8c555625 1390 Int_t i;
cc80f89e 1391 Float_t xyz[4];
8c555625 1392
1393 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1394
1395 //----------------------------------------------
1396 // Create TObjArray-s, one for each row,
1397 // each TObjArray will store the TVectors
1398 // of electrons, one TVector per each track.
1399 //----------------------------------------------
1400
f74bb6f5 1401 Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
1402 TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
8c555625 1403 for(i=0; i<nrows; i++){
1404 row[i] = new TObjArray;
f74bb6f5 1405 nofElectrons[i]=0;
1406 tracks[i]=0;
8c555625 1407 }
8c555625 1408
37831078 1409
1410
8c555625 1411 //--------------------------------------------------------------------
1412 // Loop over tracks, the "track" contains the full history
1413 //--------------------------------------------------------------------
1414
1415 Int_t previousTrack,currentTrack;
1416 previousTrack = -1; // nothing to store so far!
1417
1418 for(Int_t track=0;track<ntracks;track++){
1419
1420 ResetHits();
1421
1422 TH->GetEvent(track); // get next track
1423 Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track
1424
1425 if(nhits == 0) continue; // no hits in the TPC for this track
1426
1427 //--------------------------------------------------------------
1428 // Loop over hits
1429 //--------------------------------------------------------------
1430
1431 for(Int_t hit=0;hit<nhits;hit++){
1432
1433 tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit
1434
1435 Int_t sector=tpcHit->fSector; // sector number
cc80f89e 1436 if(sector != isec) continue;
8c555625 1437
1438 currentTrack = tpcHit->fTrack; // track number
1439 if(currentTrack != previousTrack){
1440
1441 // store already filled fTrack
1442
1443 for(i=0;i<nrows;i++){
1444 if(previousTrack != -1){
73042f01 1445 if(nofElectrons[i]>0){
cc80f89e 1446 TVector &v = *tracks[i];
8c555625 1447 v(0) = previousTrack;
73042f01 1448 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
cc80f89e 1449 row[i]->Add(tracks[i]);
8c555625 1450 }
1451 else{
cc80f89e 1452 delete tracks[i]; // delete empty TVector
1453 tracks[i]=0;
8c555625 1454 }
1455 }
1456
73042f01 1457 nofElectrons[i]=0;
cc80f89e 1458 tracks[i] = new TVector(481); // TVectors for the next fTrack
8c555625 1459
1460 } // end of loop over rows
1461
1462 previousTrack=currentTrack; // update track label
1463 }
1464
73042f01 1465 Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
8c555625 1466
1467 //---------------------------------------------------
1468 // Calculate the electron attachment probability
1469 //---------------------------------------------------
1470
cc80f89e 1471
1472 Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->fZ))
1473 /fTPCParam->GetDriftV();
8c555625 1474 // in microseconds!
73042f01 1475 Float_t attProb = fTPCParam->GetAttCoef()*
8c555625 1476 fTPCParam->GetOxyCont()*time; // fraction!
1477
1478 //-----------------------------------------------
1479 // Loop over electrons
1480 //-----------------------------------------------
cc80f89e 1481 Int_t index[3];
1482 index[1]=isec;
73042f01 1483 for(Int_t nel=0;nel<qI;nel++){
8c555625 1484 // skip if electron lost due to the attachment
73042f01 1485 if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
8c555625 1486 xyz[0]=tpcHit->fX;
1487 xyz[1]=tpcHit->fY;
cc80f89e 1488 xyz[2]=tpcHit->fZ;
1489 xyz[3]= (Float_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1490 index[0]=1;
1491
1492 TransportElectron(xyz,index); //MI change -august
73042f01 1493 Int_t rowNumber;
cc80f89e 1494 fTPCParam->GetPadRow(xyz,index); //MI change august
73042f01 1495 rowNumber = index[2];
cc80f89e 1496 //transform position to local digit coordinates
1497 //relative to nearest pad row
73042f01 1498 if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue;
1499 nofElectrons[rowNumber]++;
8c555625 1500 //----------------------------------
1501 // Expand vector if necessary
1502 //----------------------------------
73042f01 1503 if(nofElectrons[rowNumber]>120){
1504 Int_t range = tracks[rowNumber]->GetNrows();
1505 if((nofElectrons[rowNumber])>(range-1)/4){
cc80f89e 1506
73042f01 1507 tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
fe4da5cc 1508 }
1509 }
1510
73042f01 1511 TVector &v = *tracks[rowNumber];
1512 Int_t idx = 4*nofElectrons[rowNumber]-3;
8c555625 1513
cc80f89e 1514 v(idx)= xyz[0]; // X - pad row coordinate
1515 v(idx+1)=xyz[1]; // Y - pad coordinate (along the pad-row)
1516 v(idx+2)=xyz[2]; // Z - time bin coordinate
1517 v(idx+3)=xyz[3]; // avalanche size
8c555625 1518 } // end of loop over electrons
1519
1520 } // end of loop over hits
1521 } // end of loop over tracks
1522
1523 //
1524 // store remaining track (the last one) if not empty
1525 //
1526
1527 for(i=0;i<nrows;i++){
73042f01 1528 if(nofElectrons[i]>0){
cc80f89e 1529 TVector &v = *tracks[i];
8c555625 1530 v(0) = previousTrack;
73042f01 1531 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
cc80f89e 1532 row[i]->Add(tracks[i]);
fe4da5cc 1533 }
1534 else{
cc80f89e 1535 delete tracks[i];
1536 tracks[i]=0;
8c555625 1537 }
1538 }
1539
cc80f89e 1540 delete [] tracks;
73042f01 1541 delete [] nofElectrons;
8c555625 1542
8c555625 1543
cc80f89e 1544} // end of MakeSector
8c555625 1545
fe4da5cc 1546
1547//_____________________________________________________________________________
1548void AliTPC::Init()
1549{
1550 //
1551 // Initialise TPC detector after definition of geometry
1552 //
1553 Int_t i;
1554 //
1555 printf("\n");
1556 for(i=0;i<35;i++) printf("*");
1557 printf(" TPC_INIT ");
1558 for(i=0;i<35;i++) printf("*");
1559 printf("\n");
1560 //
1561 for(i=0;i<80;i++) printf("*");
1562 printf("\n");
1563}
1564
1565//_____________________________________________________________________________
1566void AliTPC::MakeBranch(Option_t* option)
1567{
1568 //
1569 // Create Tree branches for the TPC.
1570 //
1571 Int_t buffersize = 4000;
1572 char branchname[10];
1573 sprintf(branchname,"%s",GetName());
1574
1575 AliDetector::MakeBranch(option);
1576
73042f01 1577 char *d = strstr(option,"D");
fe4da5cc 1578
73042f01 1579 if (fDigits && gAlice->TreeD() && d) {
fe4da5cc 1580 gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1581 printf("Making Branch %s for digits\n",branchname);
1582 }
fe4da5cc 1583}
1584
1585//_____________________________________________________________________________
1586void AliTPC::ResetDigits()
1587{
1588 //
1589 // Reset number of digits and the digits array for this detector
fe4da5cc 1590 //
1591 fNdigits = 0;
cc80f89e 1592 if (fDigits) fDigits->Clear();
fe4da5cc 1593}
1594
1595//_____________________________________________________________________________
1596void AliTPC::SetSecAL(Int_t sec)
1597{
8c555625 1598 //---------------------------------------------------
fe4da5cc 1599 // Activate/deactivate selection for lower sectors
8c555625 1600 //---------------------------------------------------
1601
1602 //-----------------------------------------------------------------
1603 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1604 //-----------------------------------------------------------------
1605
fe4da5cc 1606 fSecAL = sec;
1607}
1608
1609//_____________________________________________________________________________
1610void AliTPC::SetSecAU(Int_t sec)
1611{
8c555625 1612 //----------------------------------------------------
fe4da5cc 1613 // Activate/deactivate selection for upper sectors
8c555625 1614 //---------------------------------------------------
1615
1616 //-----------------------------------------------------------------
1617 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1618 //-----------------------------------------------------------------
1619
fe4da5cc 1620 fSecAU = sec;
1621}
1622
1623//_____________________________________________________________________________
1624void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
1625{
8c555625 1626 //----------------------------------------
fe4da5cc 1627 // Select active lower sectors
8c555625 1628 //----------------------------------------
1629
1630 //-----------------------------------------------------------------
1631 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1632 //-----------------------------------------------------------------
1633
fe4da5cc 1634 fSecLows[0] = s1;
1635 fSecLows[1] = s2;
1636 fSecLows[2] = s3;
1637 fSecLows[3] = s4;
1638 fSecLows[4] = s5;
1639 fSecLows[5] = s6;
1640}
1641
1642//_____________________________________________________________________________
1643void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
1644 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
1645 Int_t s11 , Int_t s12)
1646{
8c555625 1647 //--------------------------------
fe4da5cc 1648 // Select active upper sectors
8c555625 1649 //--------------------------------
1650
1651 //-----------------------------------------------------------------
1652 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1653 //-----------------------------------------------------------------
1654
fe4da5cc 1655 fSecUps[0] = s1;
1656 fSecUps[1] = s2;
1657 fSecUps[2] = s3;
1658 fSecUps[3] = s4;
1659 fSecUps[4] = s5;
1660 fSecUps[5] = s6;
1661 fSecUps[6] = s7;
1662 fSecUps[7] = s8;
1663 fSecUps[8] = s9;
1664 fSecUps[9] = s10;
1665 fSecUps[10] = s11;
1666 fSecUps[11] = s12;
1667}
1668
1669//_____________________________________________________________________________
1670void AliTPC::SetSens(Int_t sens)
1671{
8c555625 1672
1673 //-------------------------------------------------------------
1674 // Activates/deactivates the sensitive strips at the center of
1675 // the pad row -- this is for the space-point resolution calculations
1676 //-------------------------------------------------------------
1677
1678 //-----------------------------------------------------------------
1679 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1680 //-----------------------------------------------------------------
1681
fe4da5cc 1682 fSens = sens;
1683}
4b0fdcad 1684
73042f01 1685void AliTPC::SetSide(Float_t side=0.)
4b0fdcad 1686{
73042f01 1687 // choice of the TPC side
1688
4b0fdcad 1689 fSide = side;
1690
1691}
1283eee5 1692//____________________________________________________________________________
1693void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
1694 Float_t p2,Float_t p3)
1695{
fe4da5cc 1696
73042f01 1697 // gax mixture definition
1698
1283eee5 1699 fNoComp = nc;
1700
1701 fMixtComp[0]=c1;
1702 fMixtComp[1]=c2;
1703 fMixtComp[2]=c3;
1704
1705 fMixtProp[0]=p1;
1706 fMixtProp[1]=p2;
1707 fMixtProp[2]=p3;
1708
1709
cc80f89e 1710}
1711//_____________________________________________________________________________
1712
1713void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
1714{
1715 //
1716 // electron transport taking into account:
1717 // 1. diffusion,
1718 // 2.ExB at the wires
1719 // 3. nonisochronity
1720 //
1721 // xyz and index must be already transformed to system 1
1722 //
1723
1724 fTPCParam->Transform1to2(xyz,index);
1725
1726 //add diffusion
1727 Float_t driftl=xyz[2];
1728 if(driftl<0.01) driftl=0.01;
1729 driftl=TMath::Sqrt(driftl);
73042f01 1730 Float_t sigT = driftl*(fTPCParam->GetDiffT());
1731 Float_t sigL = driftl*(fTPCParam->GetDiffL());
1732 xyz[0]=gRandom->Gaus(xyz[0],sigT);
1733 xyz[1]=gRandom->Gaus(xyz[1],sigT);
1734 xyz[2]=gRandom->Gaus(xyz[2],sigL);
cc80f89e 1735
1736 // ExB
1737
1738 if (fTPCParam->GetMWPCReadout()==kTRUE){
1739 Float_t x1=xyz[0];
1740 fTPCParam->Transform2to2NearestWire(xyz,index);
1741 Float_t dx=xyz[0]-x1;
1742 xyz[1]+=dx*(fTPCParam->GetOmegaTau());
1743 }
1744 //add nonisochronity (not implemented yet)
1745
1283eee5 1746}
fe4da5cc 1747//_____________________________________________________________________________
1748void AliTPC::Streamer(TBuffer &R__b)
1749{
1750 //
1751 // Stream an object of class AliTPC.
1752 //
1753 if (R__b.IsReading()) {
1754 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1755 AliDetector::Streamer(R__b);
1756 if (R__v < 2) return;
1757 R__b >> fNsectors;
fe4da5cc 1758 } else {
1759 R__b.WriteVersion(AliTPC::IsA());
1760 AliDetector::Streamer(R__b);
1761 R__b << fNsectors;
fe4da5cc 1762 }
1763}
1764
fe4da5cc 1765ClassImp(AliTPCdigit)
1766
1767//_____________________________________________________________________________
1768AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
1769 AliDigit(tracks)
1770{
1771 //
1772 // Creates a TPC digit object
1773 //
1774 fSector = digits[0];
1775 fPadRow = digits[1];
1776 fPad = digits[2];
1777 fTime = digits[3];
1778 fSignal = digits[4];
1779}
1780
1781
1782ClassImp(AliTPChit)
1783
1784//_____________________________________________________________________________
1785AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1786AliHit(shunt,track)
1787{
1788 //
1789 // Creates a TPC hit object
1790 //
1791 fSector = vol[0];
1792 fPadRow = vol[1];
1793 fX = hits[0];
1794 fY = hits[1];
1795 fZ = hits[2];
1796 fQ = hits[3];
1797}
1798
8c555625 1799