]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPC.cxx
Bug fix after sect renamed to sector
[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
88cb7938 16/* $Id$ */
4c039060 17
fe4da5cc 18///////////////////////////////////////////////////////////////////////////////
19// //
20// Time Projection Chamber //
21// This class contains the basic functions for the Time Projection Chamber //
22// detector. Functions specific to one particular geometry are //
23// contained in the derived classes //
24// //
25//Begin_Html
26/*
1439f98e 27<img src="picts/AliTPCClass.gif">
fe4da5cc 28*/
29//End_Html
30// //
8c555625 31// //
fe4da5cc 32///////////////////////////////////////////////////////////////////////////////
33
73042f01 34//
35
88cb7938 36#include <Riostream.h>
37#include <stdlib.h>
38
39#include <TFile.h>
40#include <TGeometry.h>
41#include <TInterpreter.h>
fe4da5cc 42#include <TMath.h>
e8d02863 43#include <TMatrixF.h>
44#include <TVector.h>
fe4da5cc 45#include <TNode.h>
fe4da5cc 46#include <TObjectTable.h>
88cb7938 47#include <TParticle.h>
afc42102 48#include <TROOT.h>
88cb7938 49#include <TRandom.h>
afc42102 50#include <TSystem.h>
88cb7938 51#include <TTUBS.h>
52#include <TTree.h>
88cb7938 53#include <TVirtualMC.h>
2a336e15 54#include <TString.h>
55#include <TF2.h>
4c57c771 56#include <TStopwatch.h>
cc80f89e 57
cc80f89e 58#include "AliDigits.h"
88cb7938 59#include "AliMagF.h"
60#include "AliPoints.h"
61#include "AliRun.h"
62#include "AliRunLoader.h"
cc80f89e 63#include "AliSimDigits.h"
88cb7938 64#include "AliTPC.h"
65#include "AliTPC.h"
88cb7938 66#include "AliTPCDigitsArray.h"
67#include "AliTPCLoader.h"
68#include "AliTPCPRF2D.h"
69#include "AliTPCParamSR.h"
70#include "AliTPCRF1D.h"
be5ffbfe 71//#include "AliTPCTrackHits.h"
792bb11c 72#include "AliTPCTrackHitsV2.h"
88cb7938 73#include "AliTrackReference.h"
5d12ce38 74#include "AliMC.h"
85a5290f 75#include "AliTPCDigitizer.h"
0421c3d1 76#include "AliTPCBuffer.h"
77#include "AliTPCDDLRawData.h"
8c2b3fd7 78#include "AliLog.h"
39c8eb58 79
80
fe4da5cc 81ClassImp(AliTPC)
fe4da5cc 82//_____________________________________________________________________________
83AliTPC::AliTPC()
84{
85 //
86 // Default constructor
87 //
88 fIshunt = 0;
fe4da5cc 89 fHits = 0;
90 fDigits = 0;
fe4da5cc 91 fNsectors = 0;
cc80f89e 92 fDigitsArray = 0;
afc42102 93 fDefaults = 0;
792bb11c 94 fTrackHits = 0;
be5ffbfe 95 // fTrackHitsOld = 0;
8c2b3fd7 96#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
97 fHitType = 4; // ROOT containers
98#else
99 fHitType = 2; //default CONTAINERS - based on ROOT structure
100#endif
407ff276 101 fTPCParam = 0;
102 fNoiseTable = 0;
792bb11c 103 fActiveSectors =0;
407ff276 104
fe4da5cc 105}
106
107//_____________________________________________________________________________
108AliTPC::AliTPC(const char *name, const char *title)
109 : AliDetector(name,title)
110{
111 //
112 // Standard constructor
113 //
114
115 //
116 // Initialise arrays of hits and digits
117 fHits = new TClonesArray("AliTPChit", 176);
5d12ce38 118 gAlice->GetMCApp()->AddHitList(fHits);
cc80f89e 119 fDigitsArray = 0;
afc42102 120 fDefaults = 0;
fe4da5cc 121 //
792bb11c 122 fTrackHits = new AliTPCTrackHitsV2;
39c8eb58 123 fTrackHits->SetHitPrecision(0.002);
124 fTrackHits->SetStepPrecision(0.003);
792bb11c 125 fTrackHits->SetMaxDistance(100);
126
be5ffbfe 127 //fTrackHitsOld = new AliTPCTrackHits; //MI - 13.09.2000
128 //fTrackHitsOld->SetHitPrecision(0.002);
129 //fTrackHitsOld->SetStepPrecision(0.003);
130 //fTrackHitsOld->SetMaxDistance(100);
792bb11c 131
407ff276 132 fNoiseTable =0;
39c8eb58 133
8c2b3fd7 134#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
135 fHitType = 4; // ROOT containers
136#else
da33556f 137 fHitType = 2;
8c2b3fd7 138#endif
792bb11c 139 fActiveSectors = 0;
39c8eb58 140 //
fe4da5cc 141 // Initialise counters
cc80f89e 142 fNsectors = 0;
cc80f89e 143
fe4da5cc 144 //
145 fIshunt = 0;
146 //
147 // Initialise color attributes
148 SetMarkerColor(kYellow);
73042f01 149
150 //
151 // Set TPC parameters
152 //
153
afc42102 154
155 if (!strcmp(title,"Default")) {
156 fTPCParam = new AliTPCParamSR;
73042f01 157 } else {
8c2b3fd7 158 AliWarning("In Config.C you must set non-default parameters.");
73042f01 159 fTPCParam=0;
160 }
161
fe4da5cc 162}
163
164//_____________________________________________________________________________
9bdd974b 165AliTPC::AliTPC(const AliTPC& t):AliDetector(t){
166 //
167 // dummy copy constructor
168 //
169}
fe4da5cc 170AliTPC::~AliTPC()
171{
172 //
173 // TPC destructor
174 //
73042f01 175
fe4da5cc 176 fIshunt = 0;
177 delete fHits;
178 delete fDigits;
73042f01 179 delete fTPCParam;
39c8eb58 180 delete fTrackHits; //MI 15.09.2000
be5ffbfe 181 // delete fTrackHitsOld; //MI 10.12.2001
407ff276 182 if (fNoiseTable) delete [] fNoiseTable;
183
fe4da5cc 184}
185
fe4da5cc 186//_____________________________________________________________________________
187void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
188{
189 //
190 // Add a hit to the list
191 //
39c8eb58 192 if (fHitType&1){
193 TClonesArray &lhits = *fHits;
194 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
195 }
792bb11c 196 if (fHitType>1)
8c2b3fd7 197 AddHit2(track,vol,hits);
fe4da5cc 198}
88cb7938 199
fe4da5cc 200//_____________________________________________________________________________
201void AliTPC::BuildGeometry()
202{
cc80f89e 203
fe4da5cc 204 //
205 // Build TPC ROOT TNode geometry for the event display
206 //
73042f01 207 TNode *nNode, *nTop;
fe4da5cc 208 TTUBS *tubs;
209 Int_t i;
210 const int kColorTPC=19;
1283eee5 211 char name[5], title[25];
fe4da5cc 212 const Double_t kDegrad=TMath::Pi()/180;
1283eee5 213 const Double_t kRaddeg=180./TMath::Pi();
214
1283eee5 215
73042f01 216 Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
217 Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
1283eee5 218
73042f01 219 Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
220 Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
1283eee5 221
222 Int_t nLo = fTPCParam->GetNInnerSector()/2;
223 Int_t nHi = fTPCParam->GetNOuterSector()/2;
224
73042f01 225 const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
226 const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
227 const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
228 const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
1283eee5 229
230
73042f01 231 const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
232 const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
1283eee5 233
234 Double_t rl,ru;
235
236
fe4da5cc 237 //
238 // Get ALICE top node
fe4da5cc 239 //
1283eee5 240
73042f01 241 nTop=gAlice->GetGeometry()->GetNode("alice");
1283eee5 242
243 // inner sectors
244
cc80f89e 245 rl = fTPCParam->GetInnerRadiusLow();
246 ru = fTPCParam->GetInnerRadiusUp();
1283eee5 247
248
fe4da5cc 249 for(i=0;i<nLo;i++) {
250 sprintf(name,"LS%2.2d",i);
1283eee5 251 name[4]='\0';
252 sprintf(title,"TPC low sector %3d",i);
253 title[24]='\0';
254
73042f01 255 tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
256 kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
fe4da5cc 257 tubs->SetNumberOfDivisions(1);
73042f01 258 nTop->cd();
259 nNode = new TNode(name,title,name,0,0,0,"");
260 nNode->SetLineColor(kColorTPC);
261 fNodes->Add(nNode);
fe4da5cc 262 }
1283eee5 263
fe4da5cc 264 // Outer sectors
1283eee5 265
cc80f89e 266 rl = fTPCParam->GetOuterRadiusLow();
267 ru = fTPCParam->GetOuterRadiusUp();
1283eee5 268
fe4da5cc 269 for(i=0;i<nHi;i++) {
270 sprintf(name,"US%2.2d",i);
1283eee5 271 name[4]='\0';
fe4da5cc 272 sprintf(title,"TPC upper sector %d",i);
1283eee5 273 title[24]='\0';
73042f01 274 tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
275 khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
fe4da5cc 276 tubs->SetNumberOfDivisions(1);
73042f01 277 nTop->cd();
278 nNode = new TNode(name,title,name,0,0,0,"");
279 nNode->SetLineColor(kColorTPC);
280 fNodes->Add(nNode);
fe4da5cc 281 }
cc80f89e 282
73042f01 283}
1283eee5 284
fe4da5cc 285//_____________________________________________________________________________
286void AliTPC::CreateMaterials()
287{
8c555625 288 //-----------------------------------------------
37831078 289 // Create Materials for for TPC simulations
8c555625 290 //-----------------------------------------------
291
292 //-----------------------------------------------------------------
293 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
294 //-----------------------------------------------------------------
1283eee5 295
c1637e41 296 Int_t iSXFLD=gAlice->Field()->Integ();
73042f01 297 Float_t sXMGMX=gAlice->Field()->Max();
1283eee5 298
299 Float_t amat[5]; // atomic numbers
300 Float_t zmat[5]; // z
301 Float_t wmat[5]; // proportions
302
303 Float_t density;
37831078 304
1283eee5 305
1283eee5 306
c1637e41 307 //***************** Gases *************************
1283eee5 308
37831078 309
1283eee5 310 //--------------------------------------------------------------
c1637e41 311 // gases - air and CO2
1283eee5 312 //--------------------------------------------------------------
313
37831078 314 // CO2
1283eee5 315
316 amat[0]=12.011;
317 amat[1]=15.9994;
318
319 zmat[0]=6.;
320 zmat[1]=8.;
321
c1637e41 322 wmat[0]=0.2729;
323 wmat[1]=0.7271;
1283eee5 324
325 density=0.001977;
326
1283eee5 327
c1637e41 328 AliMixture(10,"CO2",amat,zmat,density,2,wmat);
329 //
330 // Air
331 //
332 amat[0]=15.9994;
333 amat[1]=14.007;
334 //
335 zmat[0]=8.;
336 zmat[1]=7.;
337 //
338 wmat[0]=0.233;
339 wmat[1]=0.767;
340 //
341 density=0.001205;
1283eee5 342
c1637e41 343 AliMixture(11,"Air",amat,zmat,density,2,wmat);
344
1283eee5 345 //----------------------------------------------------------------
c1637e41 346 // drift gases
1283eee5 347 //----------------------------------------------------------------
37831078 348
37831078 349
350 // Drift gases 1 - nonsensitive, 2 - sensitive
c1637e41 351 // Ne-CO2-N (85-10-5)
e4e763b9 352
353 amat[0]= 20.18;
354 amat[1]=12.011;
355 amat[2]=15.9994;
c1637e41 356 amat[3]=14.007;
e4e763b9 357
358 zmat[0]= 10.;
359 zmat[1]=6.;
360 zmat[2]=8.;
c1637e41 361 zmat[3]=7.;
e4e763b9 362
c1637e41 363 wmat[0]=0.7707;
364 wmat[1]=0.0539;
365 wmat[2]=0.1438;
366 wmat[3]=0.0316;
367
368 density=0.0010252;
e4e763b9 369
c1637e41 370 AliMixture(12,"Ne-CO2-N-1",amat,zmat,density,4,wmat);
371 AliMixture(13,"Ne-CO2-N-2",amat,zmat,density,4,wmat);
1283eee5 372
373 //----------------------------------------------------------------------
374 // solid materials
375 //----------------------------------------------------------------------
376
1283eee5 377
37831078 378 // Kevlar C14H22O2N2
1283eee5 379
37831078 380 amat[0] = 12.011;
381 amat[1] = 1.;
382 amat[2] = 15.999;
383 amat[3] = 14.006;
1283eee5 384
37831078 385 zmat[0] = 6.;
386 zmat[1] = 1.;
387 zmat[2] = 8.;
388 zmat[3] = 7.;
389
390 wmat[0] = 14.;
391 wmat[1] = 22.;
392 wmat[2] = 2.;
393 wmat[3] = 2.;
394
395 density = 1.45;
396
c1637e41 397 AliMixture(14,"Kevlar",amat,zmat,density,-4,wmat);
37831078 398
399 // NOMEX
1283eee5 400
37831078 401 amat[0] = 12.011;
402 amat[1] = 1.;
403 amat[2] = 15.999;
404 amat[3] = 14.006;
405
406 zmat[0] = 6.;
407 zmat[1] = 1.;
408 zmat[2] = 8.;
409 zmat[3] = 7.;
410
411 wmat[0] = 14.;
412 wmat[1] = 22.;
413 wmat[2] = 2.;
414 wmat[3] = 2.;
415
c1637e41 416 density = 0.029;
417
418 AliMixture(15,"NOMEX",amat,zmat,density,-4,wmat);
37831078 419
420 // Makrolon C16H18O3
421
422 amat[0] = 12.011;
423 amat[1] = 1.;
424 amat[2] = 15.999;
1283eee5 425
37831078 426 zmat[0] = 6.;
427 zmat[1] = 1.;
428 zmat[2] = 8.;
429
430 wmat[0] = 16.;
431 wmat[1] = 18.;
432 wmat[2] = 3.;
433
434 density = 1.2;
435
c1637e41 436 AliMixture(16,"Makrolon",amat,zmat,density,-3,wmat);
437
438 // Tedlar C2H3F
439
440 amat[0] = 12.011;
441 amat[1] = 1.;
442 amat[2] = 18.998;
443
444 zmat[0] = 6.;
445 zmat[1] = 1.;
446 zmat[2] = 9.;
447
448 wmat[0] = 2.;
449 wmat[1] = 3.;
450 wmat[2] = 1.;
451
452 density = 1.71;
453
454 AliMixture(17, "Tedlar",amat,zmat,density,-3,wmat);
37831078 455
1283eee5 456 // Mylar C5H4O2
457
458 amat[0]=12.011;
459 amat[1]=1.;
460 amat[2]=15.9994;
461
462 zmat[0]=6.;
463 zmat[1]=1.;
464 zmat[2]=8.;
465
466 wmat[0]=5.;
467 wmat[1]=4.;
468 wmat[2]=2.;
469
470 density = 1.39;
fe4da5cc 471
c1637e41 472 AliMixture(18, "Mylar",amat,zmat,density,-3,wmat);
473 // material for "prepregs"
474 // Epoxy - C14 H20 O3
475 // Quartz SiO2
476 // Carbon C
477 // prepreg1 60% C-fiber, 40% epoxy (vol)
478 amat[0]=12.011;
479 amat[1]=1.;
480 amat[2]=15.994;
1283eee5 481
c1637e41 482 zmat[0]=6.;
483 zmat[1]=1.;
484 zmat[2]=8.;
1283eee5 485
c1637e41 486 wmat[0]=0.923;
487 wmat[1]=0.023;
488 wmat[2]=0.054;
1283eee5 489
c1637e41 490 density=1.859;
1283eee5 491
c1637e41 492 AliMixture(19, "Prepreg1",amat,zmat,density,3,wmat);
1283eee5 493
c1637e41 494 //prepreg2 60% glass-fiber, 40% epoxy
1283eee5 495
c1637e41 496 amat[0]=12.01;
497 amat[1]=1.;
498 amat[2]=15.994;
499 amat[3]=28.086;
500
501 zmat[0]=6.;
502 zmat[1]=1.;
503 zmat[2]=8.;
504 zmat[3]=14.;
505
506 wmat[0]=0.194;
507 wmat[1]=0.023;
508 wmat[2]=0.443;
509 wmat[3]=0.340;
510
511 density=1.82;
512
513 AliMixture(20, "Prepreg2",amat,zmat,density,4,wmat);
514
515 //prepreg3 50% glass-fiber, 50% epoxy
516
517 amat[0]=12.01;
518 amat[1]=1.;
519 amat[2]=15.994;
520 amat[3]=28.086;
521
522 zmat[0]=6.;
523 zmat[1]=1.;
524 zmat[2]=8.;
525 zmat[3]=14.;
1283eee5 526
c1637e41 527 wmat[0]=0.225;
528 wmat[1]=0.03;
529 wmat[2]=0.443;
530 wmat[3]=0.3;
531
532 density=1.163;
533
534 AliMixture(21, "Prepreg3",amat,zmat,density,4,wmat);
535
536 // G10 60% SiO2 40% epoxy
537
538 amat[0]=12.01;
539 amat[1]=1.;
540 amat[2]=15.994;
541 amat[3]=28.086;
542
543 zmat[0]=6.;
544 zmat[1]=1.;
545 zmat[2]=8.;
546 zmat[3]=14.;
547
548 wmat[0]=0.194;
549 wmat[1]=0.023;
550 wmat[2]=0.443;
551 wmat[3]=0.340;
552
553 density=1.7;
554
555 AliMixture(22, "G10",amat,zmat,density,4,wmat);
556
37831078 557 // Al
1283eee5 558
37831078 559 amat[0] = 26.98;
560 zmat[0] = 13.;
1283eee5 561
37831078 562 density = 2.7;
1283eee5 563
c1637e41 564 AliMaterial(23,"Al",amat[0],zmat[0],density,999.,999.);
1283eee5 565
c1637e41 566 // Si (for electronics
1283eee5 567
37831078 568 amat[0] = 28.086;
9a3667bd 569 zmat[0] = 14.;
1283eee5 570
37831078 571 density = 2.33;
1283eee5 572
c1637e41 573 AliMaterial(24,"Si",amat[0],zmat[0],density,999.,999.);
1283eee5 574
37831078 575 // Cu
1283eee5 576
37831078 577 amat[0] = 63.546;
578 zmat[0] = 29.;
1283eee5 579
37831078 580 density = 8.96;
1283eee5 581
c1637e41 582 AliMaterial(25,"Cu",amat[0],zmat[0],density,999.,999.);
1283eee5 583
c1637e41 584 // Epoxy - C14 H20 O3
585
37831078 586 amat[0]=12.011;
587 amat[1]=1.;
588 amat[2]=15.9994;
1283eee5 589
37831078 590 zmat[0]=6.;
591 zmat[1]=1.;
592 zmat[2]=8.;
1283eee5 593
c1637e41 594 wmat[0]=14.;
595 wmat[1]=20.;
596 wmat[2]=3.;
1283eee5 597
c1637e41 598 density=1.25;
1283eee5 599
c1637e41 600 AliMixture(26,"Epoxy",amat,zmat,density,-3,wmat);
1283eee5 601
c1637e41 602 // Plexiglas C5H8O2
1283eee5 603
5a28a08f 604 amat[0]=12.011;
605 amat[1]=1.;
606 amat[2]=15.9994;
607
608 zmat[0]=6.;
609 zmat[1]=1.;
610 zmat[2]=8.;
611
c1637e41 612 wmat[0]=5.;
613 wmat[1]=8.;
614 wmat[2]=2.;
5a28a08f 615
c1637e41 616 density=1.18;
5a28a08f 617
c1637e41 618 AliMixture(27,"Plexiglas",amat,zmat,density,-3,wmat);
37831078 619
ba1d05f9 620 // Carbon
621
622 amat[0]=12.011;
623 zmat[0]=6.;
624 density= 2.265;
625
c1637e41 626 AliMaterial(28,"C",amat[0],zmat[0],density,999.,999.);
ba1d05f9 627
c1637e41 628 // Fe (steel for the inner heat screen)
e4e763b9 629
c1637e41 630 amat[0]=55.845;
e4e763b9 631
c1637e41 632 zmat[0]=26.;
e4e763b9 633
c1637e41 634 density=7.87;
ba1d05f9 635
c1637e41 636 AliMaterial(29,"Fe",amat[0],zmat[0],density,999.,999.);
ba1d05f9 637
37831078 638 //----------------------------------------------------------
639 // tracking media for gases
640 //----------------------------------------------------------
641
c1637e41 642 AliMedium(0, "Air", 11, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
643 AliMedium(1, "Ne-CO2-N-1", 12, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
644 AliMedium(2, "Ne-CO2-N-2", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
37831078 645 AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
646
647 //-----------------------------------------------------------
648 // tracking media for solids
649 //-----------------------------------------------------------
650
c1637e41 651 AliMedium(4,"Al",23,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
652 AliMedium(5,"Kevlar",14,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
653 AliMedium(6,"Nomex",15,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
654 AliMedium(7,"Makrolon",16,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
655 AliMedium(8,"Mylar",18,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
656 AliMedium(9,"Tedlar",17,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
657 //
658 AliMedium(10,"Prepreg1",19,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
659 AliMedium(11,"Prepreg2",20,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
660 AliMedium(12,"Prepreg3",21,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
661 AliMedium(13,"Epoxy",26,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
662
663 AliMedium(14,"Cu",25,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
664 AliMedium(15,"Si",24,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
665 AliMedium(16,"G10",22,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
666 AliMedium(17,"Plexiglas",27,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
667 AliMedium(18,"Steel",29,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
1283eee5 668
fe4da5cc 669}
670
407ff276 671void AliTPC::GenerNoise(Int_t tablesize)
672{
673 //
674 //Generate table with noise
675 //
676 if (fTPCParam==0) {
677 // error message
678 fNoiseDepth=0;
679 return;
680 }
681 if (fNoiseTable) delete[] fNoiseTable;
682 fNoiseTable = new Float_t[tablesize];
683 fNoiseDepth = tablesize;
684 fCurrentNoise =0; //!index of the noise in the noise table
685
686 Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
687 for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);
688}
689
690Float_t AliTPC::GetNoise()
691{
692 // get noise from table
693 // if ((fCurrentNoise%10)==0)
694 // fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
695 if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
696 return fNoiseTable[fCurrentNoise++];
697 //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
698}
699
700
9bdd974b 701Bool_t AliTPC::IsSectorActive(Int_t sec) const
792bb11c 702{
703 //
704 // check if the sector is active
705 if (!fActiveSectors) return kTRUE;
706 else return fActiveSectors[sec];
707}
708
709void AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
710{
711 // activate interesting sectors
88cb7938 712 SetTreeAddress();//just for security
792bb11c 713 if (fActiveSectors) delete [] fActiveSectors;
714 fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
715 for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
716 for (Int_t i=0;i<n;i++)
717 if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector()) fActiveSectors[sectors[i]]=kTRUE;
718
719}
720
12d8d654 721void AliTPC::SetActiveSectors(Int_t flag)
792bb11c 722{
723 //
724 // activate sectors which were hitted by tracks
725 //loop over tracks
88cb7938 726 SetTreeAddress();//just for security
792bb11c 727 if (fHitType==0) return; // if Clones hit - not short volume ID information
728 if (fActiveSectors) delete [] fActiveSectors;
729 fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
12d8d654 730 if (flag) {
731 for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
732 return;
733 }
792bb11c 734 for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
735 TBranch * branch=0;
88cb7938 736 if (TreeH() == 0x0)
737 {
8c2b3fd7 738 AliFatal("Can not find TreeH in folder");
88cb7938 739 return;
740 }
741 if (fHitType>1) branch = TreeH()->GetBranch("TPC2");
742 else branch = TreeH()->GetBranch("TPC");
743 Stat_t ntracks = TreeH()->GetEntries();
792bb11c 744 // loop over all hits
8c2b3fd7 745 AliDebug(1,Form("Got %d tracks",ntracks));
88cb7938 746
8c2b3fd7 747 for(Int_t track=0;track<ntracks;track++) {
792bb11c 748 ResetHits();
749 //
750 if (fTrackHits && fHitType&4) {
88cb7938 751 TBranch * br1 = TreeH()->GetBranch("fVolumes");
752 TBranch * br2 = TreeH()->GetBranch("fNVolumes");
792bb11c 753 br1->GetEvent(track);
754 br2->GetEvent(track);
755 Int_t *volumes = fTrackHits->GetVolumes();
756 for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++)
757 fActiveSectors[volumes[j]]=kTRUE;
758 }
759
760 //
be5ffbfe 761// if (fTrackHitsOld && fHitType&2) {
762// TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
763// br->GetEvent(track);
764// AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
765// for (UInt_t j=0;j<ar->GetSize();j++){
766// fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
767// }
768// }
792bb11c 769 }
792bb11c 770}
771
772
773
fe4da5cc 774
0421c3d1 775//_____________________________________________________________________________
776void AliTPC::Digits2Raw()
777{
778// convert digits of the current event to raw data
779
780 static const Int_t kThreshold = 0;
781 static const Bool_t kCompress = kTRUE;
782
783 fLoader->LoadDigits();
784 TTree* digits = fLoader->TreeD();
785 if (!digits) {
8c2b3fd7 786 AliError("No digits tree");
0421c3d1 787 return;
788 }
789
790 AliSimDigits digarr;
791 AliSimDigits* digrow = &digarr;
792 digits->GetBranch("Segment")->SetAddress(&digrow);
793
794 const char* fileName = "AliTPCDDL.dat";
795 AliTPCBuffer* buffer = new AliTPCBuffer(fileName);
796 //Verbose level
797 // 0: Silent
798 // 1: cout messages
799 // 2: txt files with digits
800 //BE CAREFUL, verbose level 2 MUST be used only for debugging and
801 //it is highly suggested to use this mode only for debugging digits files
802 //reasonably small, because otherwise the size of the txt files can reach
803 //quickly several MB wasting time and disk space.
804 buffer->SetVerbose(0);
805
806 Int_t nEntries = Int_t(digits->GetEntries());
807 Int_t previousSector = -1;
808 Int_t subSector = 0;
809 for (Int_t i = 0; i < nEntries; i++) {
810 digits->GetEntry(i);
811 Int_t sector, row;
812 fTPCParam->AdjustSectorRow(digarr.GetID(), sector, row);
813 if(previousSector != sector) {
814 subSector = 0;
815 previousSector = sector;
816 }
817
818 if (sector < 36) { //inner sector [0;35]
819 if (row != 30) {
820 //the whole row is written into the output file
821 buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
822 sector, subSector, row);
823 } else {
824 //only the pads in the range [37;48] are written into the output file
825 buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 1,
826 sector, subSector, row);
827 subSector = 1;
828 //only the pads outside the range [37;48] are written into the output file
829 buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 2,
830 sector, subSector, row);
831 }//end else
832
833 } else { //outer sector [36;71]
834 if (row == 54) subSector = 2;
835 if ((row != 27) && (row != 76)) {
836 buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
837 sector, subSector, row);
838 } else if (row == 27) {
8c2b3fd7 839 //only the pads outside the range [43;46] are written into the output file
840 buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 2,
0421c3d1 841 sector, subSector, row);
8c2b3fd7 842 subSector = 1;
843 //only the pads in the range [43;46] are written into the output file
844 buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 1,
0421c3d1 845 sector, subSector, row);
846 } else if (row == 76) {
8c2b3fd7 847 //only the pads outside the range [33;88] are written into the output file
848 buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 2,
849 sector, subSector, row);
850 subSector = 3;
851 //only the pads in the range [33;88] are written into the output file
852 buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 1,
0421c3d1 853 sector, subSector, row);
854 }
855 }//end else
856 }//end for
857
858 delete buffer;
859 fLoader->UnloadDigits();
860
861 AliTPCDDLRawData rawWriter;
862 rawWriter.SetVerbose(0);
863
864 rawWriter.RawData(fileName);
865 gSystem->Unlink(fileName);
866
867 if (kCompress) {
8c2b3fd7 868 AliInfo("Compressing raw data");
0421c3d1 869 rawWriter.RawDataCompDecompress(kTRUE);
870 gSystem->Unlink("Statistics");
871 }
872}
873
874
0a61bf9d 875
85a5290f 876//______________________________________________________________________
9bdd974b 877AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
85a5290f 878{
879 return new AliTPCDigitizer(manager);
880}
0a61bf9d 881//__
176aff27 882void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)
0a61bf9d 883{
884 //create digits from summable digits
407ff276 885 GenerNoise(500000); //create teble with noise
0a61bf9d 886
887 //conect tree with sSDigits
88cb7938 888 TTree *t = fLoader->TreeS();
889
8c2b3fd7 890 if (t == 0x0) {
891 fLoader->LoadSDigits("READ");
892 t = fLoader->TreeS();
893 if (t == 0x0) {
894 AliError("Can not get input TreeS");
895 return;
896 }
897 }
88cb7938 898
899 if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D");
900
0a61bf9d 901 AliSimDigits digarr, *dummy=&digarr;
88cb7938 902 TBranch* sdb = t->GetBranch("Segment");
8c2b3fd7 903 if (sdb == 0x0) {
904 AliError("Can not find branch with segments in TreeS.");
905 return;
906 }
88cb7938 907
908 sdb->SetAddress(&dummy);
909
0a61bf9d 910 Stat_t nentries = t->GetEntries();
911
5f16d237 912 // set zero suppression
913
914 fTPCParam->SetZeroSup(2);
915
916 // get zero suppression
917
918 Int_t zerosup = fTPCParam->GetZeroSup();
919
920 //make tree with digits
921
0a61bf9d 922 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
923 arr->SetClass("AliSimDigits");
924 arr->Setup(fTPCParam);
88cb7938 925 arr->MakeTree(fLoader->TreeD());
0a61bf9d 926
88cb7938 927 AliTPCParam * par = fTPCParam;
5f16d237 928
0a61bf9d 929 //Loop over segments of the TPC
5f16d237 930
0a61bf9d 931 for (Int_t n=0; n<nentries; n++) {
932 t->GetEvent(n);
933 Int_t sec, row;
934 if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
8c2b3fd7 935 AliWarning(Form("Invalid segment ID ! %d",digarr.GetID()));
0a61bf9d 936 continue;
937 }
8c2b3fd7 938 if (!IsSectorActive(sec)) continue;
939
0a61bf9d 940 AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
941 Int_t nrows = digrow->GetNRows();
942 Int_t ncols = digrow->GetNCols();
943
944 digrow->ExpandBuffer();
945 digarr.ExpandBuffer();
946 digrow->ExpandTrackBuffer();
947 digarr.ExpandTrackBuffer();
948
5f16d237 949
407ff276 950 Short_t * pamp0 = digarr.GetDigits();
951 Int_t * ptracks0 = digarr.GetTracks();
952 Short_t * pamp1 = digrow->GetDigits();
953 Int_t * ptracks1 = digrow->GetTracks();
954 Int_t nelems =nrows*ncols;
955 Int_t saturation = fTPCParam->GetADCSat();
956 //use internal structure of the AliDigits - for speed reason
957 //if you cahnge implementation
958 //of the Alidigits - it must be rewriten -
959 for (Int_t i= 0; i<nelems; i++){
407ff276 960 Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
961 if (q>zerosup){
962 if (q>saturation) q=saturation;
963 *pamp1=(Short_t)q;
8c2b3fd7 964
407ff276 965 ptracks1[0]=ptracks0[0];
966 ptracks1[nelems]=ptracks0[nelems];
967 ptracks1[2*nelems]=ptracks0[2*nelems];
968 }
969 pamp0++;
970 pamp1++;
971 ptracks0++;
972 ptracks1++;
973 }
5f16d237 974
5f16d237 975 arr->StoreRow(sec,row);
976 arr->ClearRow(sec,row);
5f16d237 977 }
0a61bf9d 978
0a61bf9d 979
5f16d237 980 //write results
88cb7938 981 fLoader->WriteDigits("OVERWRITE");
5f16d237 982
88cb7938 983 delete arr;
afc42102 984}
985//__________________________________________________________________
986void AliTPC::SetDefaults(){
9bdd974b 987 //
988 // setting the defaults
989 //
afc42102 990
afc42102 991 // Set response functions
992
88cb7938 993 //
024a7e64 994 AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
88cb7938 995 rl->CdGAFile();
2ab0c725 996 AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
7a09f434 997 if(param){
8c2b3fd7 998 AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
7a09f434 999 delete param;
1000 param = new AliTPCParamSR();
1001 }
1002 else {
1003 param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1004 }
1005 if(!param){
8c2b3fd7 1006 AliFatal("No TPC parameters found");
7a09f434 1007 }
1008
1009
2ab0c725 1010 AliTPCPRF2D * prfinner = new AliTPCPRF2D;
f03e3423 1011 AliTPCPRF2D * prfouter1 = new AliTPCPRF2D;
1012 AliTPCPRF2D * prfouter2 = new AliTPCPRF2D;
2ab0c725 1013 AliTPCRF1D * rf = new AliTPCRF1D(kTRUE);
2ab0c725 1014 rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1015 rf->SetOffset(3*param->GetZSigma());
1016 rf->Update();
afc42102 1017
1018 TDirectory *savedir=gDirectory;
2ab0c725 1019 TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
8c2b3fd7 1020 if (!f->IsOpen())
1021 AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !");
2a336e15 1022
1023 TString s;
2ab0c725 1024 prfinner->Read("prf_07504_Gati_056068_d02");
2a336e15 1025 //PH Set different names
1026 s=prfinner->GetGRF()->GetName();
1027 s+="in";
1028 prfinner->GetGRF()->SetName(s.Data());
1029
f03e3423 1030 prfouter1->Read("prf_10006_Gati_047051_d03");
2a336e15 1031 s=prfouter1->GetGRF()->GetName();
1032 s+="out1";
1033 prfouter1->GetGRF()->SetName(s.Data());
1034
f03e3423 1035 prfouter2->Read("prf_15006_Gati_047051_d03");
2a336e15 1036 s=prfouter2->GetGRF()->GetName();
1037 s+="out2";
1038 prfouter2->GetGRF()->SetName(s.Data());
1039
2ab0c725 1040 f->Close();
afc42102 1041 savedir->cd();
1042
2ab0c725 1043 param->SetInnerPRF(prfinner);
f03e3423 1044 param->SetOuter1PRF(prfouter1);
1045 param->SetOuter2PRF(prfouter2);
2ab0c725 1046 param->SetTimeRF(rf);
1047
afc42102 1048 // set fTPCParam
1049
2ab0c725 1050 SetParam(param);
1051
afc42102 1052
1053 fDefaults = 1;
1054
1055}
1056//__________________________________________________________________
85a5290f 1057void AliTPC::Hits2Digits()
1058{
9bdd974b 1059 //
1060 // creates digits from hits
1061 //
1062
85a5290f 1063 fLoader->LoadHits("read");
1064 fLoader->LoadDigits("recreate");
1065 AliRunLoader* runLoader = fLoader->GetRunLoader();
1066
1067 for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1068 runLoader->GetEvent(iEvent);
1069 SetActiveSectors();
1070 Hits2Digits(iEvent);
1071 }
1072
1073 fLoader->UnloadHits();
1074 fLoader->UnloadDigits();
1075}
1076//__________________________________________________________________
afc42102 1077void AliTPC::Hits2Digits(Int_t eventnumber)
1078{
1079 //----------------------------------------------------
1080 // Loop over all sectors for a single event
1081 //----------------------------------------------------
024a7e64 1082 AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
88cb7938 1083 rl->GetEvent(eventnumber);
8c2b3fd7 1084 if (fLoader->TreeH() == 0x0) {
1085 if(fLoader->LoadHits()) {
1086 AliError("Can not load hits.");
1087 }
1088 }
88cb7938 1089 SetTreeAddress();
1090
8c2b3fd7 1091 if (fLoader->TreeD() == 0x0 ) {
1092 fLoader->MakeTree("D");
1093 if (fLoader->TreeD() == 0x0 ) {
1094 AliError("Can not get TreeD");
1095 return;
1096 }
1097 }
afc42102 1098
1099 if(fDefaults == 0) SetDefaults(); // check if the parameters are set
407ff276 1100 GenerNoise(500000); //create teble with noise
2ab0c725 1101
1102 //setup TPCDigitsArray
afc42102 1103
1104 if(GetDigitsArray()) delete GetDigitsArray();
8c2b3fd7 1105
2ab0c725 1106 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1107 arr->SetClass("AliSimDigits");
afc42102 1108 arr->Setup(fTPCParam);
88cb7938 1109
1110 arr->MakeTree(fLoader->TreeD());
2ab0c725 1111 SetDigitsArray(arr);
1112
afc42102 1113 fDigitsSwitch=0; // standard digits
2ab0c725 1114
8c2b3fd7 1115 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1116 if (IsSectorActive(isec)) {
1117 AliDebug(1,Form("Hits2Digits","Sector %d is active.",isec));
1118 Hits2DigitsSector(isec);
1119 }
1120 else {
1121 AliDebug(1,Form("Hits2Digits","Sector %d is NOT active.",isec));
1122 }
1123
88cb7938 1124 fLoader->WriteDigits("OVERWRITE");
f2a509af 1125
1126//this line prevents the crash in the similar one
1127//on the beginning of this method
1128//destructor attempts to reset the tree, which is deleted by the loader
1129//need to be redesign
8c2b3fd7 1130 if(GetDigitsArray()) delete GetDigitsArray();
1131 SetDigitsArray(0x0);
f2a509af 1132
8c555625 1133}
1134
f8cf550c 1135//__________________________________________________________________
0a61bf9d 1136void AliTPC::Hits2SDigits2(Int_t eventnumber)
f8cf550c 1137{
1138
1139 //-----------------------------------------------------------
1140 // summable digits - 16 bit "ADC", no noise, no saturation
1141 //-----------------------------------------------------------
1142
8c2b3fd7 1143 //----------------------------------------------------
1144 // Loop over all sectors for a single event
1145 //----------------------------------------------------
88cb7938 1146
1147 AliRunLoader* rl = fLoader->GetRunLoader();
1148
1149 rl->GetEvent(eventnumber);
8c2b3fd7 1150 if (fLoader->TreeH() == 0x0) {
1151 if(fLoader->LoadHits()) {
1152 AliError("Can not load hits.");
1153 return;
1154 }
1155 }
88cb7938 1156 SetTreeAddress();
f8cf550c 1157
afc42102 1158
8c2b3fd7 1159 if (fLoader->TreeS() == 0x0 ) {
1160 fLoader->MakeTree("S");
1161 }
88cb7938 1162
afc42102 1163 if(fDefaults == 0) SetDefaults();
88cb7938 1164
407ff276 1165 GenerNoise(500000); //create table with noise
afc42102 1166 //setup TPCDigitsArray
1167
1168 if(GetDigitsArray()) delete GetDigitsArray();
1169
88cb7938 1170
afc42102 1171 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1172 arr->SetClass("AliSimDigits");
1173 arr->Setup(fTPCParam);
88cb7938 1174 arr->MakeTree(fLoader->TreeS());
1175
afc42102 1176 SetDigitsArray(arr);
1177
afc42102 1178 fDigitsSwitch=1; // summable digits
5f16d237 1179
1180 // set zero suppression to "0"
1181
1182 fTPCParam->SetZeroSup(0);
f8cf550c 1183
8c2b3fd7 1184 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1185 if (IsSectorActive(isec)) {
1186 Hits2DigitsSector(isec);
1187 }
88cb7938 1188
8c2b3fd7 1189 fLoader->WriteSDigits("OVERWRITE");
88cb7938 1190
1191//this line prevents the crash in the similar one
1192//on the beginning of this method
1193//destructor attempts to reset the tree, which is deleted by the loader
1194//need to be redesign
8c2b3fd7 1195 if(GetDigitsArray()) delete GetDigitsArray();
1196 SetDigitsArray(0x0);
f8cf550c 1197}
0a61bf9d 1198//__________________________________________________________________
88cb7938 1199
0a61bf9d 1200void AliTPC::Hits2SDigits()
1201{
1202
1203 //-----------------------------------------------------------
1204 // summable digits - 16 bit "ADC", no noise, no saturation
1205 //-----------------------------------------------------------
1206
85a5290f 1207 fLoader->LoadHits("read");
1208 fLoader->LoadSDigits("recreate");
1209 AliRunLoader* runLoader = fLoader->GetRunLoader();
0a61bf9d 1210
85a5290f 1211 for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1212 runLoader->GetEvent(iEvent);
1213 SetTreeAddress();
1214 SetActiveSectors();
1215 Hits2SDigits2(iEvent);
1216 }
0a61bf9d 1217
85a5290f 1218 fLoader->UnloadHits();
1219 fLoader->UnloadSDigits();
0a61bf9d 1220}
fe4da5cc 1221//_____________________________________________________________________________
88cb7938 1222
8c555625 1223void AliTPC::Hits2DigitsSector(Int_t isec)
fe4da5cc 1224{
8c555625 1225 //-------------------------------------------------------------------
fe4da5cc 1226 // TPC conversion from hits to digits.
8c555625 1227 //-------------------------------------------------------------------
1228
1229 //-----------------------------------------------------------------
1230 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1231 //-----------------------------------------------------------------
1232
fe4da5cc 1233 //-------------------------------------------------------
8c555625 1234 // Get the access to the track hits
fe4da5cc 1235 //-------------------------------------------------------
8c555625 1236
afc42102 1237 // check if the parameters are set - important if one calls this method
1238 // directly, not from the Hits2Digits
1239
1240 if(fDefaults == 0) SetDefaults();
cc80f89e 1241
88cb7938 1242 TTree *tH = TreeH(); // pointer to the hits tree
8c2b3fd7 1243 if (tH == 0x0) {
1244 AliFatal("Can not find TreeH in folder");
1245 return;
1246 }
88cb7938 1247
73042f01 1248 Stat_t ntracks = tH->GetEntries();
8c555625 1249
1250 if( ntracks > 0){
1251
1252 //-------------------------------------------
1253 // Only if there are any tracks...
1254 //-------------------------------------------
1255
8c555625 1256 TObjArray **row;
fe4da5cc 1257
8c2b3fd7 1258 Int_t nrows =fTPCParam->GetNRow(isec);
8c555625 1259
8c2b3fd7 1260 row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
fe4da5cc 1261
8c2b3fd7 1262 MakeSector(isec,nrows,tH,ntracks,row);
8c555625 1263
8c2b3fd7 1264 //--------------------------------------------------------
1265 // Digitize this sector, row by row
1266 // row[i] is the pointer to the TObjArray of TVectors,
1267 // each one containing electrons accepted on this
1268 // row, assigned into tracks
1269 //--------------------------------------------------------
8c555625 1270
8c2b3fd7 1271 Int_t i;
8c555625 1272
8c2b3fd7 1273 if (fDigitsArray->GetTree()==0) {
1274 AliFatal("Tree not set in fDigitsArray");
1275 }
8c555625 1276
8c2b3fd7 1277 for (i=0;i<nrows;i++){
1278
1279 AliDigits * dig = fDigitsArray->CreateRow(isec,i);
8c555625 1280
8c2b3fd7 1281 DigitizeRow(i,isec,row);
8c555625 1282
8c2b3fd7 1283 fDigitsArray->StoreRow(isec,i);
8c555625 1284
8c2b3fd7 1285 Int_t ndig = dig->GetDigitSize();
88cb7938 1286
8c2b3fd7 1287 AliDebug(10,
1288 Form("*** Sector, row, compressed digits %d %d %d ***\n",
1289 isec,i,ndig));
792bb11c 1290
8c2b3fd7 1291 fDigitsArray->ClearRow(isec,i);
8c555625 1292
cc80f89e 1293
8c2b3fd7 1294 } // end of the sector digitization
8c555625 1295
8c2b3fd7 1296 for(i=0;i<nrows+2;i++){
1297 row[i]->Delete();
1298 delete row[i];
1299 }
cc80f89e 1300
8c2b3fd7 1301 delete [] row; // delete the array of pointers to TObjArray-s
8c555625 1302
1303 } // ntracks >0
8c555625 1304
cc80f89e 1305} // end of Hits2DigitsSector
8c555625 1306
8c555625 1307
8c555625 1308//_____________________________________________________________________________
cc80f89e 1309void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
8c555625 1310{
1311 //-----------------------------------------------------------
1312 // Single row digitization, coupling from the neighbouring
1313 // rows taken into account
1314 //-----------------------------------------------------------
1315
1316 //-----------------------------------------------------------------
1317 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
cc80f89e 1318 // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
8c555625 1319 //-----------------------------------------------------------------
1320
8c555625 1321 Float_t zerosup = fTPCParam->GetZeroSup();
8c2b3fd7 1322
cc80f89e 1323 fCurrentIndex[1]= isec;
8c555625 1324
8c555625 1325
73042f01 1326 Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1327 Int_t nofTbins = fTPCParam->GetMaxTBin();
1328 Int_t indexRange[4];
8c555625 1329 //
1330 // Integrated signal for this row
1331 // and a single track signal
cc80f89e 1332 //
de61d5d5 1333
e8d02863 1334 TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated
1335 TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single
8c555625 1336 //
e8d02863 1337 TMatrixF &total = *m1;
8c555625 1338
1339 // Array of pointers to the label-signal list
1340
73042f01 1341 Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1342 Float_t **pList = new Float_t* [nofDigits];
8c555625 1343
1344 Int_t lp;
cc80f89e 1345 Int_t i1;
73042f01 1346 for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
8c555625 1347 //
cc80f89e 1348 //calculate signal
1349 //
b584c7dd 1350 Int_t row1=irow;
1351 Int_t row2=irow+2;
cc80f89e 1352 for (Int_t row= row1;row<=row2;row++){
1353 Int_t nTracks= rows[row]->GetEntries();
1354 for (i1=0;i1<nTracks;i1++){
1355 fCurrentIndex[2]= row;
b584c7dd 1356 fCurrentIndex[3]=irow+1;
1357 if (row==irow+1){
cc80f89e 1358 m2->Zero(); // clear single track signal matrix
73042f01 1359 Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1360 GetList(trackLabel,nofPads,m2,indexRange,pList);
cc80f89e 1361 }
73042f01 1362 else GetSignal(rows[row],i1,0,m1,indexRange);
cc80f89e 1363 }
8c555625 1364 }
cc80f89e 1365
8c555625 1366 Int_t tracks[3];
8c555625 1367
cc80f89e 1368 AliDigits *dig = fDigitsArray->GetRow(isec,irow);
de61d5d5 1369 Int_t gi=-1;
1370 Float_t fzerosup = zerosup+0.5;
1371 for(Int_t it=0;it<nofTbins;it++){
de61d5d5 1372 for(Int_t ip=0;ip<nofPads;ip++){
1373 gi++;
8c2b3fd7 1374 Float_t q=total(ip,it);
f8cf550c 1375 if(fDigitsSwitch == 0){
407ff276 1376 q+=GetNoise();
de61d5d5 1377 if(q <=fzerosup) continue; // do not fill zeros
68771f83 1378 q = TMath::Nint(q);
f8cf550c 1379 if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation
cc80f89e 1380
f8cf550c 1381 }
1382
1383 else {
8c2b3fd7 1384 if(q <= 0.) continue; // do not fill zeros
1385 if(q>2000.) q=2000.;
1386 q *= 16.;
1387 q = TMath::Nint(q);
f8cf550c 1388 }
8c555625 1389
1390 //
1391 // "real" signal or electronic noise (list = -1)?
1392 //
1393
1394 for(Int_t j1=0;j1<3;j1++){
407ff276 1395 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
8c555625 1396 }
1397
cc80f89e 1398//Begin_Html
1399/*
1400 <A NAME="AliDigits"></A>
1401 using of AliDigits object
1402*/
1403//End_Html
1404 dig->SetDigitFast((Short_t)q,it,ip);
8c2b3fd7 1405 if (fDigitsArray->IsSimulated()) {
1406 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1407 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1408 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1409 }
8c555625 1410
1411 } // end of loop over time buckets
1412 } // end of lop over pads
1413
1414 //
1415 // This row has been digitized, delete nonused stuff
1416 //
1417
73042f01 1418 for(lp=0;lp<nofDigits;lp++){
8c555625 1419 if(pList[lp]) delete [] pList[lp];
1420 }
1421
1422 delete [] pList;
1423
1424 delete m1;
1425 delete m2;
8c555625 1426
1427} // end of DigitizeRow
cc80f89e 1428
8c555625 1429//_____________________________________________________________________________
cc80f89e 1430
de61d5d5 1431Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
e8d02863 1432 TMatrixF *m1, TMatrixF *m2,Int_t *indexRange)
8c555625 1433{
1434
1435 //---------------------------------------------------------------
1436 // Calculates 2-D signal (pad,time) for a single track,
1437 // returns a pointer to the signal matrix and the track label
1438 // No digitization is performed at this level!!!
1439 //---------------------------------------------------------------
1440
1441 //-----------------------------------------------------------------
1442 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
cc80f89e 1443 // Modified: Marian Ivanov
8c555625 1444 //-----------------------------------------------------------------
1445
8c2b3fd7 1446 TVector *tv;
de61d5d5 1447
8c2b3fd7 1448 tv = (TVector*)p1->At(ntr); // pointer to a track
1449 TVector &v = *tv;
8c555625 1450
1451 Float_t label = v(0);
b584c7dd 1452 Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1)-1)/2;
8c555625 1453
e61fd20d 1454 Int_t nElectrons = (tv->GetNrows()-1)/5;
73042f01 1455 indexRange[0]=9999; // min pad
1456 indexRange[1]=-1; // max pad
1457 indexRange[2]=9999; //min time
1458 indexRange[3]=-1; // max time
8c555625 1459
e8d02863 1460 TMatrixF &signal = *m1;
1461 TMatrixF &total = *m2;
8c555625 1462 //
1463 // Loop over all electrons
1464 //
8c555625 1465 for(Int_t nel=0; nel<nElectrons; nel++){
e61fd20d 1466 Int_t idx=nel*5;
cc80f89e 1467 Float_t aval = v(idx+4);
1468 Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
e61fd20d 1469 Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)};
de61d5d5 1470 Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
1471
1472 Int_t *index = fTPCParam->GetResBin(0);
1473 Float_t *weight = & (fTPCParam->GetResWeight(0));
1474
1475 if (n>0) for (Int_t i =0; i<n; i++){
8c2b3fd7 1476 Int_t pad=index[1]+centralPad; //in digit coordinates central pad has coordinate 0
1477
1478 if (pad>=0){
1479 Int_t time=index[2];
1480 Float_t qweight = *(weight)*eltoadcfac;
1481
1482 if (m1!=0) signal(pad,time)+=qweight;
1483 total(pad,time)+=qweight;
1484 if (indexRange[0]>pad) indexRange[0]=pad;
1485 if (indexRange[1]<pad) indexRange[1]=pad;
1486 if (indexRange[2]>time) indexRange[2]=time;
1487 if (indexRange[3]<time) indexRange[3]=time;
1488
1489 index+=3;
1490 weight++;
1491
1492 }
cc80f89e 1493 }
8c555625 1494 } // end of loop over electrons
cc80f89e 1495
8c555625 1496 return label; // returns track label when finished
1497}
1498
1499//_____________________________________________________________________________
e8d02863 1500void AliTPC::GetList(Float_t label,Int_t np,TMatrixF *m,
de61d5d5 1501 Int_t *indexRange, Float_t **pList)
8c555625 1502{
1503 //----------------------------------------------------------------------
1504 // Updates the list of tracks contributing to digits for a given row
1505 //----------------------------------------------------------------------
1506
1507 //-----------------------------------------------------------------
1508 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1509 //-----------------------------------------------------------------
1510
e8d02863 1511 TMatrixF &signal = *m;
8c555625 1512
1513 // lop over nonzero digits
1514
73042f01 1515 for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1516 for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
8c555625 1517
1518
8c2b3fd7 1519 // accept only the contribution larger than 500 electrons (1/2 s_noise)
921bf71a 1520
8c2b3fd7 1521 if(signal(ip,it)<0.5) continue;
921bf71a 1522
8c2b3fd7 1523 Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
8c555625 1524
8c2b3fd7 1525 if(!pList[globalIndex]){
8c555625 1526
8c2b3fd7 1527 //
1528 // Create new list (6 elements - 3 signals and 3 labels),
1529 //
8c555625 1530
8c2b3fd7 1531 pList[globalIndex] = new Float_t [6];
8c555625 1532
8c2b3fd7 1533 // set list to -1
1534
1535 *pList[globalIndex] = -1.;
1536 *(pList[globalIndex]+1) = -1.;
1537 *(pList[globalIndex]+2) = -1.;
1538 *(pList[globalIndex]+3) = -1.;
1539 *(pList[globalIndex]+4) = -1.;
1540 *(pList[globalIndex]+5) = -1.;
1541
1542 *pList[globalIndex] = label;
1543 *(pList[globalIndex]+3) = signal(ip,it);
1544 }
1545 else {
8c555625 1546
8c2b3fd7 1547 // check the signal magnitude
8c555625 1548
8c2b3fd7 1549 Float_t highest = *(pList[globalIndex]+3);
1550 Float_t middle = *(pList[globalIndex]+4);
1551 Float_t lowest = *(pList[globalIndex]+5);
1552
1553 //
1554 // compare the new signal with already existing list
1555 //
1556
1557 if(signal(ip,it)<lowest) continue; // neglect this track
8c555625 1558
8c2b3fd7 1559 //
8c555625 1560
8c2b3fd7 1561 if (signal(ip,it)>highest){
1562 *(pList[globalIndex]+5) = middle;
1563 *(pList[globalIndex]+4) = highest;
1564 *(pList[globalIndex]+3) = signal(ip,it);
1565
1566 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1567 *(pList[globalIndex]+1) = *pList[globalIndex];
1568 *pList[globalIndex] = label;
1569 }
1570 else if (signal(ip,it)>middle){
1571 *(pList[globalIndex]+5) = middle;
1572 *(pList[globalIndex]+4) = signal(ip,it);
1573
1574 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1575 *(pList[globalIndex]+1) = label;
1576 }
1577 else{
1578 *(pList[globalIndex]+5) = signal(ip,it);
1579 *(pList[globalIndex]+2) = label;
1580 }
1581 }
1582
8c555625 1583 } // end of loop over pads
1584 } // end of loop over time bins
1585
8c555625 1586}//end of GetList
1587//___________________________________________________________________
1588void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1589 Stat_t ntracks,TObjArray **row)
1590{
1591
1592 //-----------------------------------------------------------------
1593 // Prepares the sector digitization, creates the vectors of
1594 // tracks for each row of this sector. The track vector
1595 // contains the track label and the position of electrons.
1596 //-----------------------------------------------------------------
1597
1598 //-----------------------------------------------------------------
1599 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1600 //-----------------------------------------------------------------
1601
cc80f89e 1602 Float_t gasgain = fTPCParam->GetGasGain();
8c555625 1603 Int_t i;
e61fd20d 1604 Float_t xyz[5];
8c555625 1605
1606 AliTPChit *tpcHit; // pointer to a sigle TPC hit
39c8eb58 1607 //MI change
1608 TBranch * branch=0;
792bb11c 1609 if (fHitType>1) branch = TH->GetBranch("TPC2");
39c8eb58 1610 else branch = TH->GetBranch("TPC");
1611
8c555625 1612
1613 //----------------------------------------------
1614 // Create TObjArray-s, one for each row,
8c2b3fd7 1615 // each TObjArray will store the TVectors
1616 // of electrons, one TVectors per each track.
8c555625 1617 //----------------------------------------------
1618
b584c7dd 1619 Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
8c2b3fd7 1620 TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
de61d5d5 1621
b584c7dd 1622 for(i=0; i<nrows+2; i++){
8c555625 1623 row[i] = new TObjArray;
f74bb6f5 1624 nofElectrons[i]=0;
1625 tracks[i]=0;
8c555625 1626 }
8c555625 1627
37831078 1628
1629
8c555625 1630 //--------------------------------------------------------------------
1631 // Loop over tracks, the "track" contains the full history
1632 //--------------------------------------------------------------------
8c2b3fd7 1633
8c555625 1634 Int_t previousTrack,currentTrack;
1635 previousTrack = -1; // nothing to store so far!
1636
1637 for(Int_t track=0;track<ntracks;track++){
39c8eb58 1638 Bool_t isInSector=kTRUE;
8c555625 1639 ResetHits();
792bb11c 1640 isInSector = TrackInVolume(isec,track);
39c8eb58 1641 if (!isInSector) continue;
1642 //MI change
1643 branch->GetEntry(track); // get next track
8c2b3fd7 1644
39c8eb58 1645 //M.I. changes
8c555625 1646
39c8eb58 1647 tpcHit = (AliTPChit*)FirstHit(-1);
8c555625 1648
1649 //--------------------------------------------------------------
1650 // Loop over hits
1651 //--------------------------------------------------------------
1652
8c555625 1653
39c8eb58 1654 while(tpcHit){
8c555625 1655
1656 Int_t sector=tpcHit->fSector; // sector number
39c8eb58 1657 if(sector != isec){
1658 tpcHit = (AliTPChit*) NextHit();
1659 continue;
1660 }
8c555625 1661
8c2b3fd7 1662 currentTrack = tpcHit->Track(); // track number
39c8eb58 1663
8c2b3fd7 1664 if(currentTrack != previousTrack){
8c555625 1665
8c2b3fd7 1666 // store already filled fTrack
8c555625 1667
8c2b3fd7 1668 for(i=0;i<nrows+2;i++){
1669 if(previousTrack != -1){
1670 if(nofElectrons[i]>0){
1671 TVector &v = *tracks[i];
1672 v(0) = previousTrack;
e61fd20d 1673 tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
8c2b3fd7 1674 row[i]->Add(tracks[i]);
1675 }
1676 else {
1677 delete tracks[i]; // delete empty TVector
1678 tracks[i]=0;
1679 }
1680 }
1681
1682 nofElectrons[i]=0;
e61fd20d 1683 tracks[i] = new TVector(601); // TVectors for the next fTrack
8c2b3fd7 1684
1685 } // end of loop over rows
8c555625 1686
8c2b3fd7 1687 previousTrack=currentTrack; // update track label
1688 }
8c555625 1689
8c2b3fd7 1690 Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
8c555625 1691
8c2b3fd7 1692 //---------------------------------------------------
1693 // Calculate the electron attachment probability
1694 //---------------------------------------------------
8c555625 1695
cc80f89e 1696
8c2b3fd7 1697 Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
1698 /fTPCParam->GetDriftV();
1699 // in microseconds!
1700 Float_t attProb = fTPCParam->GetAttCoef()*
1701 fTPCParam->GetOxyCont()*time; // fraction!
8c555625 1702
8c2b3fd7 1703 //-----------------------------------------------
1704 // Loop over electrons
1705 //-----------------------------------------------
1706 Int_t index[3];
1707 index[1]=isec;
1708 for(Int_t nel=0;nel<qI;nel++){
1709 // skip if electron lost due to the attachment
1710 if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1711 xyz[0]=tpcHit->X();
1712 xyz[1]=tpcHit->Y();
1713 xyz[2]=tpcHit->Z();
1714 //
1715 // protection for the nonphysical avalanche size (10**6 maximum)
1716 //
1717 Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
1718 xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
1719 index[0]=1;
1720
1721 TransportElectron(xyz,index);
1722 Int_t rowNumber;
1723 fTPCParam->GetPadRow(xyz,index);
e61fd20d 1724 // Electron track time (for pileup simulation)
1725 xyz[4] = tpcHit->Time()/fTPCParam->GetTSample();
8c2b3fd7 1726 // row 0 - cross talk from the innermost row
1727 // row fNRow+1 cross talk from the outermost row
1728 rowNumber = index[2]+1;
1729 //transform position to local digit coordinates
1730 //relative to nearest pad row
1731 if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
1732 Float_t x1,y1;
1733 if (isec <fTPCParam->GetNInnerSector()) {
1734 x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
1735 y1 = fTPCParam->GetYInner(rowNumber);
1736 }
1737 else{
1738 x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
1739 y1 = fTPCParam->GetYOuter(rowNumber);
1740 }
1741 // gain inefficiency at the wires edges - linear
1742 x1=TMath::Abs(x1);
1743 y1-=1.;
1744 if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));
1745
1746 nofElectrons[rowNumber]++;
1747 //----------------------------------
1748 // Expand vector if necessary
1749 //----------------------------------
1750 if(nofElectrons[rowNumber]>120){
1751 Int_t range = tracks[rowNumber]->GetNrows();
e61fd20d 1752 if((nofElectrons[rowNumber])>(range-1)/5){
8c2b3fd7 1753
e61fd20d 1754 tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons
fe4da5cc 1755 }
8c2b3fd7 1756 }
1757
1758 TVector &v = *tracks[rowNumber];
e61fd20d 1759 Int_t idx = 5*nofElectrons[rowNumber]-4;
8c2b3fd7 1760 Real_t * position = &(((TVector&)v)(idx)); //make code faster
e61fd20d 1761 memcpy(position,xyz,5*sizeof(Float_t));
8c2b3fd7 1762
1763 } // end of loop over electrons
39c8eb58 1764
8c2b3fd7 1765 tpcHit = (AliTPChit*)NextHit();
1766
1767 } // end of loop over hits
1768 } // end of loop over tracks
8c555625 1769
1770 //
1771 // store remaining track (the last one) if not empty
1772 //
8c2b3fd7 1773
1774 for(i=0;i<nrows+2;i++){
1775 if(nofElectrons[i]>0){
1776 TVector &v = *tracks[i];
1777 v(0) = previousTrack;
e61fd20d 1778 tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
8c2b3fd7 1779 row[i]->Add(tracks[i]);
1780 }
1781 else{
1782 delete tracks[i];
1783 tracks[i]=0;
1784 }
1785 }
1786
1787 delete [] tracks;
1788 delete [] nofElectrons;
8c555625 1789
cc80f89e 1790} // end of MakeSector
8c555625 1791
fe4da5cc 1792
1793//_____________________________________________________________________________
1794void AliTPC::Init()
1795{
1796 //
1797 // Initialise TPC detector after definition of geometry
1798 //
8c2b3fd7 1799 AliDebug(1,"*********************************************");
fe4da5cc 1800}
1801
1802//_____________________________________________________________________________
88cb7938 1803void AliTPC::MakeBranch(Option_t* option)
fe4da5cc 1804{
1805 //
1806 // Create Tree branches for the TPC.
1807 //
8c2b3fd7 1808 AliDebug(1,"");
fe4da5cc 1809 Int_t buffersize = 4000;
1810 char branchname[10];
1811 sprintf(branchname,"%s",GetName());
88cb7938 1812
1813 const char *h = strstr(option,"H");
8c2b3fd7 1814
88cb7938 1815 if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
1816
1817 AliDetector::MakeBranch(option);
8c2b3fd7 1818
5cf7bbad 1819 const char *d = strstr(option,"D");
8c2b3fd7 1820
1821 if (fDigits && fLoader->TreeD() && d) {
1822 MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0);
1823 }
fe4da5cc 1824
88cb7938 1825 if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000
fe4da5cc 1826}
1827
1828//_____________________________________________________________________________
1829void AliTPC::ResetDigits()
1830{
1831 //
1832 // Reset number of digits and the digits array for this detector
fe4da5cc 1833 //
1834 fNdigits = 0;
cc80f89e 1835 if (fDigits) fDigits->Clear();
fe4da5cc 1836}
1837
fe4da5cc 1838
fe4da5cc 1839
1840//_____________________________________________________________________________
1841void AliTPC::SetSens(Int_t sens)
1842{
8c555625 1843
1844 //-------------------------------------------------------------
1845 // Activates/deactivates the sensitive strips at the center of
1846 // the pad row -- this is for the space-point resolution calculations
1847 //-------------------------------------------------------------
1848
1849 //-----------------------------------------------------------------
1850 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1851 //-----------------------------------------------------------------
1852
fe4da5cc 1853 fSens = sens;
1854}
2b06d5c3 1855
4b0fdcad 1856
73042f01 1857void AliTPC::SetSide(Float_t side=0.)
4b0fdcad 1858{
73042f01 1859 // choice of the TPC side
1860
4b0fdcad 1861 fSide = side;
1862
1863}
cc80f89e 1864//_____________________________________________________________________________
1865
1866void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
1867{
1868 //
1869 // electron transport taking into account:
1870 // 1. diffusion,
1871 // 2.ExB at the wires
1872 // 3. nonisochronity
1873 //
1874 // xyz and index must be already transformed to system 1
1875 //
1876
1877 fTPCParam->Transform1to2(xyz,index);
1878
1879 //add diffusion
1880 Float_t driftl=xyz[2];
1881 if(driftl<0.01) driftl=0.01;
1882 driftl=TMath::Sqrt(driftl);
73042f01 1883 Float_t sigT = driftl*(fTPCParam->GetDiffT());
1884 Float_t sigL = driftl*(fTPCParam->GetDiffL());
1885 xyz[0]=gRandom->Gaus(xyz[0],sigT);
1886 xyz[1]=gRandom->Gaus(xyz[1],sigT);
1887 xyz[2]=gRandom->Gaus(xyz[2],sigL);
cc80f89e 1888
1889 // ExB
1890
1891 if (fTPCParam->GetMWPCReadout()==kTRUE){
b584c7dd 1892 Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
cc80f89e 1893 xyz[1]+=dx*(fTPCParam->GetOmegaTau());
1894 }
b584c7dd 1895 //add nonisochronity (not implemented yet)
1283eee5 1896}
fe4da5cc 1897
fe4da5cc 1898ClassImp(AliTPChit)
1899
1900//_____________________________________________________________________________
1901AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1902AliHit(shunt,track)
1903{
1904 //
1905 // Creates a TPC hit object
1906 //
1907 fSector = vol[0];
1908 fPadRow = vol[1];
1909 fX = hits[0];
1910 fY = hits[1];
1911 fZ = hits[2];
1912 fQ = hits[3];
e61fd20d 1913 fTime = hits[4];
fe4da5cc 1914}
1915
39c8eb58 1916//________________________________________________________________________
792bb11c 1917// Additional code because of the AliTPCTrackHitsV2
39c8eb58 1918
176aff27 1919void AliTPC::MakeBranch2(Option_t *option,const char */*file*/)
39c8eb58 1920{
1921 //
1922 // Create a new branch in the current Root Tree
1923 // The branch of fHits is automatically split
1924 // MI change 14.09.2000
8c2b3fd7 1925 AliDebug(1,"");
792bb11c 1926 if (fHitType<2) return;
39c8eb58 1927 char branchname[10];
1928 sprintf(branchname,"%s2",GetName());
1929 //
1930 // Get the pointer to the header
5cf7bbad 1931 const char *cH = strstr(option,"H");
39c8eb58 1932 //
8c2b3fd7 1933 if (fTrackHits && TreeH() && cH && fHitType&4) {
1934 AliDebug(1,"Making branch for Type 4 Hits");
88cb7938 1935 TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
8c2b3fd7 1936 }
792bb11c 1937
be5ffbfe 1938// if (fTrackHitsOld && TreeH() && cH && fHitType&2) {
1939// AliDebug(1,"Making branch for Type 2 Hits");
1940// AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
1941// TreeH(),fBufferSize,99);
1942// TreeH()->GetListOfBranches()->Add(branch);
1943// }
39c8eb58 1944}
1945
1946void AliTPC::SetTreeAddress()
1947{
8c2b3fd7 1948 //Sets tree address for hits
1949 if (fHitType<=1) {
1950 if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
1951 AliDetector::SetTreeAddress();
1952 }
792bb11c 1953 if (fHitType>1) SetTreeAddress2();
39c8eb58 1954}
1955
1956void AliTPC::SetTreeAddress2()
1957{
1958 //
1959 // Set branch address for the TrackHits Tree
1960 //
8c2b3fd7 1961 AliDebug(1,"");
88cb7938 1962
39c8eb58 1963 TBranch *branch;
1964 char branchname[20];
1965 sprintf(branchname,"%s2",GetName());
1966 //
1967 // Branch address for hit tree
88cb7938 1968 TTree *treeH = TreeH();
792bb11c 1969 if ((treeH)&&(fHitType&4)) {
39c8eb58 1970 branch = treeH->GetBranch(branchname);
8c2b3fd7 1971 if (branch) {
1972 branch->SetAddress(&fTrackHits);
1973 AliDebug(1,"fHitType&4 Setting");
1974 }
f2a509af 1975 else
8c2b3fd7 1976 AliDebug(1,"fHitType&4 Failed (can not find branch)");
f2a509af 1977
39c8eb58 1978 }
be5ffbfe 1979 // if ((treeH)&&(fHitType&2)) {
1980// branch = treeH->GetBranch(branchname);
1981// if (branch) {
1982// branch->SetAddress(&fTrackHitsOld);
1983// AliDebug(1,"fHitType&2 Setting");
1984// }
1985// else
1986// AliDebug(1,"fHitType&2 Failed (can not find branch)");
1987// }
b6e0d3fe 1988 //set address to TREETR
8c2b3fd7 1989
88cb7938 1990 TTree *treeTR = TreeTR();
b6e0d3fe 1991 if (treeTR && fTrackReferences) {
1992 branch = treeTR->GetBranch(GetName());
1993 if (branch) branch->SetAddress(&fTrackReferences);
1994 }
1995
39c8eb58 1996}
1997
1998void AliTPC::FinishPrimary()
1999{
792bb11c 2000 if (fTrackHits &&fHitType&4) fTrackHits->FlushHitStack();
be5ffbfe 2001 // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->FlushHitStack();
39c8eb58 2002}
2003
2004
2005void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2006{
2007 //
2008 // add hit to the list
39c8eb58 2009 Int_t rtrack;
2010 if (fIshunt) {
5d12ce38 2011 int primary = gAlice->GetMCApp()->GetPrimary(track);
2012 gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
39c8eb58 2013 rtrack=primary;
2014 } else {
2015 rtrack=track;
5d12ce38 2016 gAlice->GetMCApp()->FlagTrack(track);
39c8eb58 2017 }
792bb11c 2018 if (fTrackHits && fHitType&4)
39c8eb58 2019 fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
e61fd20d 2020 hits[1],hits[2],(Int_t)hits[3],hits[4]);
be5ffbfe 2021 // if (fTrackHitsOld &&fHitType&2 )
2022// fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2023// hits[1],hits[2],(Int_t)hits[3]);
792bb11c 2024
39c8eb58 2025}
2026
2027void AliTPC::ResetHits()
88cb7938 2028{
39c8eb58 2029 if (fHitType&1) AliDetector::ResetHits();
792bb11c 2030 if (fHitType>1) ResetHits2();
39c8eb58 2031}
2032
2033void AliTPC::ResetHits2()
2034{
2035 //
2036 //reset hits
792bb11c 2037 if (fTrackHits && fHitType&4) fTrackHits->Clear();
be5ffbfe 2038 // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
792bb11c 2039
39c8eb58 2040}
2041
2042AliHit* AliTPC::FirstHit(Int_t track)
2043{
792bb11c 2044 if (fHitType>1) return FirstHit2(track);
39c8eb58 2045 return AliDetector::FirstHit(track);
2046}
2047AliHit* AliTPC::NextHit()
2048{
9bdd974b 2049 //
2050 // gets next hit
2051 //
792bb11c 2052 if (fHitType>1) return NextHit2();
2053
39c8eb58 2054 return AliDetector::NextHit();
2055}
2056
2057AliHit* AliTPC::FirstHit2(Int_t track)
2058{
2059 //
2060 // Initialise the hit iterator
2061 // Return the address of the first hit for track
2062 // If track>=0 the track is read from disk
2063 // while if track<0 the first hit of the current
2064 // track is returned
2065 //
2066 if(track>=0) {
2067 gAlice->ResetHits();
88cb7938 2068 TreeH()->GetEvent(track);
39c8eb58 2069 }
2070 //
792bb11c 2071 if (fTrackHits && fHitType&4) {
39c8eb58 2072 fTrackHits->First();
2073 return fTrackHits->GetHit();
2074 }
be5ffbfe 2075 // if (fTrackHitsOld && fHitType&2) {
2076// fTrackHitsOld->First();
2077// return fTrackHitsOld->GetHit();
2078// }
792bb11c 2079
39c8eb58 2080 else return 0;
2081}
2082
2083AliHit* AliTPC::NextHit2()
2084{
2085 //
2086 //Return the next hit for the current track
2087
792bb11c 2088
be5ffbfe 2089// if (fTrackHitsOld && fHitType&2) {
2090// fTrackHitsOld->Next();
2091// return fTrackHitsOld->GetHit();
2092// }
39c8eb58 2093 if (fTrackHits) {
2094 fTrackHits->Next();
2095 return fTrackHits->GetHit();
2096 }
2097 else
2098 return 0;
2099}
2100
2101void AliTPC::LoadPoints(Int_t)
2102{
2103 //
2104 Int_t a = 0;
8c2b3fd7 2105
f2e8b846 2106 if(fHitType==1) AliDetector::LoadPoints(a);
2107 else LoadPoints2(a);
39c8eb58 2108}
2109
2110
2111void AliTPC::RemapTrackHitIDs(Int_t *map)
2112{
9bdd974b 2113 //
2114 // remapping
2115 //
39c8eb58 2116 if (!fTrackHits) return;
39c8eb58 2117
be5ffbfe 2118// if (fTrackHitsOld && fHitType&2){
2119// AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2120// for (UInt_t i=0;i<arr->GetSize();i++){
2121// AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2122// info->fTrackID = map[info->fTrackID];
2123// }
2124// }
2125// if (fTrackHitsOld && fHitType&4){
2126 if (fTrackHits && fHitType&4){
792bb11c 2127 TClonesArray * arr = fTrackHits->GetArray();;
2128 for (Int_t i=0;i<arr->GetEntriesFast();i++){
2129 AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2130 info->fTrackID = map[info->fTrackID];
2131 }
2132 }
39c8eb58 2133}
2134
792bb11c 2135Bool_t AliTPC::TrackInVolume(Int_t id,Int_t track)
2136{
2137 //return bool information - is track in given volume
2138 //load only part of the track information
2139 //return true if current track is in volume
2140 //
2141 // return kTRUE;
be5ffbfe 2142 // if (fTrackHitsOld && fHitType&2) {
2143// TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
2144// br->GetEvent(track);
2145// AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2146// for (UInt_t j=0;j<ar->GetSize();j++){
2147// if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2148// }
2149// }
792bb11c 2150
2151 if (fTrackHits && fHitType&4) {
88cb7938 2152 TBranch * br1 = TreeH()->GetBranch("fVolumes");
2153 TBranch * br2 = TreeH()->GetBranch("fNVolumes");
792bb11c 2154 br2->GetEvent(track);
2155 br1->GetEvent(track);
2156 Int_t *volumes = fTrackHits->GetVolumes();
2157 Int_t nvolumes = fTrackHits->GetNVolumes();
2158 if (!volumes && nvolumes>0) {
8c2b3fd7 2159 AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes));
792bb11c 2160 return kFALSE;
2161 }
2162 for (Int_t j=0;j<nvolumes; j++)
2163 if (volumes[j]==id) return kTRUE;
2164 }
2165
2166 if (fHitType&1) {
88cb7938 2167 TBranch * br = TreeH()->GetBranch("fSector");
792bb11c 2168 br->GetEvent(track);
2169 for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2170 if ( ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2171 }
2172 }
2173 return kFALSE;
2174
2175}
39c8eb58 2176
2177//_____________________________________________________________________________
2178void AliTPC::LoadPoints2(Int_t)
2179{
2180 //
2181 // Store x, y, z of all hits in memory
2182 //
be5ffbfe 2183 // if (fTrackHits == 0 && fTrackHitsOld==0) return;
2184 if (fTrackHits == 0 ) return;
39c8eb58 2185 //
792bb11c 2186 Int_t nhits =0;
2187 if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
be5ffbfe 2188 // if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
792bb11c 2189
39c8eb58 2190 if (nhits == 0) return;
5d12ce38 2191 Int_t tracks = gAlice->GetMCApp()->GetNtrack();
39c8eb58 2192 if (fPoints == 0) fPoints = new TObjArray(tracks);
2193 AliHit *ahit;
2194 //
2195 Int_t *ntrk=new Int_t[tracks];
2196 Int_t *limi=new Int_t[tracks];
2197 Float_t **coor=new Float_t*[tracks];
2198 for(Int_t i=0;i<tracks;i++) {
2199 ntrk[i]=0;
2200 coor[i]=0;
2201 limi[i]=0;
2202 }
2203 //
2204 AliPoints *points = 0;
2205 Float_t *fp=0;
2206 Int_t trk;
2207 Int_t chunk=nhits/4+1;
2208 //
2209 // Loop over all the hits and store their position
2210 //
2211 ahit = FirstHit2(-1);
39c8eb58 2212 while (ahit){
39c8eb58 2213 trk=ahit->GetTrack();
2214 if(ntrk[trk]==limi[trk]) {
2215 //
2216 // Initialise a new track
2217 fp=new Float_t[3*(limi[trk]+chunk)];
2218 if(coor[trk]) {
2219 memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2220 delete [] coor[trk];
2221 }
2222 limi[trk]+=chunk;
2223 coor[trk] = fp;
2224 } else {
2225 fp = coor[trk];
2226 }
2227 fp[3*ntrk[trk] ] = ahit->X();
2228 fp[3*ntrk[trk]+1] = ahit->Y();
2229 fp[3*ntrk[trk]+2] = ahit->Z();
2230 ntrk[trk]++;
2231 ahit = NextHit2();
2232 }
792bb11c 2233
2234
2235
39c8eb58 2236 //
2237 for(trk=0; trk<tracks; ++trk) {
2238 if(ntrk[trk]) {
2239 points = new AliPoints();
2240 points->SetMarkerColor(GetMarkerColor());
2241 points->SetMarkerSize(GetMarkerSize());
2242 points->SetDetector(this);
2243 points->SetParticle(trk);
2244 points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
2245 fPoints->AddAt(points,trk);
2246 delete [] coor[trk];
2247 coor[trk]=0;
2248 }
2249 }
2250 delete [] coor;
2251 delete [] ntrk;
2252 delete [] limi;
2253}
2254
2255
2256//_____________________________________________________________________________
2257void AliTPC::LoadPoints3(Int_t)
2258{
2259 //
2260 // Store x, y, z of all hits in memory
2261 // - only intersection point with pad row
2262 if (fTrackHits == 0) return;
2263 //
2264 Int_t nhits = fTrackHits->GetEntriesFast();
2265 if (nhits == 0) return;
5d12ce38 2266 Int_t tracks = gAlice->GetMCApp()->GetNtrack();
39c8eb58 2267 if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2268 fPoints->Expand(2*tracks);
2269 AliHit *ahit;
2270 //
2271 Int_t *ntrk=new Int_t[tracks];
2272 Int_t *limi=new Int_t[tracks];
2273 Float_t **coor=new Float_t*[tracks];
2274 for(Int_t i=0;i<tracks;i++) {
2275 ntrk[i]=0;
2276 coor[i]=0;
2277 limi[i]=0;
2278 }
2279 //
2280 AliPoints *points = 0;
2281 Float_t *fp=0;
2282 Int_t trk;
2283 Int_t chunk=nhits/4+1;
2284 //
2285 // Loop over all the hits and store their position
2286 //
2287 ahit = FirstHit2(-1);
39c8eb58 2288
2289 Int_t lastrow = -1;
2290 while (ahit){
39c8eb58 2291 trk=ahit->GetTrack();
2292 Float_t x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2293 Int_t index[3]={1,((AliTPChit*)ahit)->fSector,0};
2294 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
2295 if (currentrow!=lastrow){
2296 lastrow = currentrow;
2297 //later calculate intersection point
2298 if(ntrk[trk]==limi[trk]) {
2299 //
2300 // Initialise a new track
2301 fp=new Float_t[3*(limi[trk]+chunk)];
2302 if(coor[trk]) {
2303 memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2304 delete [] coor[trk];
2305 }
2306 limi[trk]+=chunk;
2307 coor[trk] = fp;
2308 } else {
2309 fp = coor[trk];
2310 }
2311 fp[3*ntrk[trk] ] = ahit->X();
2312 fp[3*ntrk[trk]+1] = ahit->Y();
2313 fp[3*ntrk[trk]+2] = ahit->Z();
2314 ntrk[trk]++;
2315 }
2316 ahit = NextHit2();
2317 }
2318
2319 //
2320 for(trk=0; trk<tracks; ++trk) {
2321 if(ntrk[trk]) {
2322 points = new AliPoints();
2323 points->SetMarkerColor(GetMarkerColor()+1);
2324 points->SetMarkerStyle(5);
2325 points->SetMarkerSize(0.2);
2326 points->SetDetector(this);
2327 points->SetParticle(trk);
39c8eb58 2328 points->SetPolyMarker(ntrk[trk],coor[trk],30);
2329 fPoints->AddAt(points,tracks+trk);
2330 delete [] coor[trk];
2331 coor[trk]=0;
2332 }
2333 }
2334 delete [] coor;
2335 delete [] ntrk;
2336 delete [] limi;
2337}
2338
2339
2340
88cb7938 2341AliLoader* AliTPC::MakeLoader(const char* topfoldername)
2342{
8c2b3fd7 2343 //Makes TPC loader
2344 fLoader = new AliTPCLoader(GetName(),topfoldername);
2345 return fLoader;
88cb7938 2346}
2347
42157e55 2348////////////////////////////////////////////////////////////////////////
2349AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
2350//
2351// load TPC paarmeters from a given file or create new if the object
2352// is not found there
88cb7938 2353// 12/05/2003 This method should be moved to the AliTPCLoader
2354// and one has to decide where to store the TPC parameters
2355// M.Kowalski
42157e55 2356 char paramName[50];
2357 sprintf(paramName,"75x40_100x60_150x60");
2358 AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2359 if (paramTPC) {
8c2b3fd7 2360 AliDebugClass(1,Form("TPC parameters %s found.",paramName));
42157e55 2361 } else {
8c2b3fd7 2362 AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
42157e55 2363 paramTPC = new AliTPCParamSR;
2364 }
2365 return paramTPC;
2366
2367// the older version of parameters can be accessed with this code.
2368// In some cases, we have old parameters saved in the file but
2369// digits were created with new parameters, it can be distinguish
2370// by the name of TPC TreeD. The code here is just for the case
2371// we would need to compare with old data, uncomment it if needed.
2372//
2373// char paramName[50];
2374// sprintf(paramName,"75x40_100x60");
2375// AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2376// if (paramTPC) {
2377// cout<<"TPC parameters "<<paramName<<" found."<<endl;
2378// } else {
2379// sprintf(paramName,"75x40_100x60_150x60");
2380// paramTPC=(AliTPCParam*)in->Get(paramName);
2381// if (paramTPC) {
2382// cout<<"TPC parameters "<<paramName<<" found."<<endl;
2383// } else {
2384// cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2385// <<endl;
2386// paramTPC = new AliTPCParamSR;
2387// }
2388// }
2389// return paramTPC;
2390
2391}
2392
85a5290f 2393