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