]> 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
Fatsre algorithms
[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$
68771f83 18Revision 1.46 2001/11/07 06:47:12 kowal2
19Removed printouts
20
9fe90376 21Revision 1.45 2001/11/03 13:33:48 kowal2
22Updated algorithms in Hits2SDigits, SDigits2Digits,
23Hits2ExactClusters.
24Added method Merge
25
5f16d237 26Revision 1.44 2001/08/30 09:28:48 hristov
27TTree names are explicitly set via SetName(name) and then Write() is called
28
6bc3458d 29Revision 1.43 2001/07/28 12:02:54 hristov
30Branch split level set to 99
31
8f90302f 32Revision 1.42 2001/07/28 11:38:52 hristov
33Loop variable declared once
34
43272665 35Revision 1.41 2001/07/28 10:53:50 hristov
36Digitisation done according to the general scheme (M.Ivanov)
37
0a61bf9d 38Revision 1.40 2001/07/27 13:03:14 hristov
39Default Branch split level set to 99
40
d0f40f23 41Revision 1.39 2001/07/26 09:09:34 kowal2
42Hits2Reco method added
43
a62d28a2 44Revision 1.38 2001/07/20 14:32:43 kowal2
45Processing of many events possible now
46
afc42102 47Revision 1.37 2001/06/12 07:17:18 kowal2
48Hits2SDigits method implemented (summable digits)
49
f8cf550c 50Revision 1.36 2001/05/16 14:57:25 alibrary
51New files for folders and Stack
52
9e1a0ddb 53Revision 1.35 2001/05/08 16:02:22 kowal2
54Updated material specifications
55
ba1d05f9 56Revision 1.34 2001/05/08 15:00:15 hristov
57Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
58
b9de75e1 59Revision 1.33 2001/04/03 12:40:43 kowal2
60Removed printouts
61
5a28a08f 62Revision 1.32 2001/03/12 17:47:36 hristov
63Changes needed on Sun with CC 5.0
64
5cf7bbad 65Revision 1.31 2001/03/12 08:21:50 kowal2
66Corrected C++ bug in the material definitions
67
9a3667bd 68Revision 1.30 2001/03/01 17:34:47 kowal2
69Correction due to the accuracy problem
70
710dacf5 71Revision 1.29 2001/02/28 16:34:40 kowal2
72Protection against nonphysical values of the avalanche size,
7310**6 is the maximum
74
1705494b 75Revision 1.28 2001/01/26 19:57:19 hristov
76Major upgrade of AliRoot code
77
2ab0c725 78Revision 1.27 2001/01/13 17:29:33 kowal2
79Sun compiler correction
80
f2e8b846 81Revision 1.26 2001/01/10 07:59:43 kowal2
82Corrections to load points from the noncompressed hits.
83
53048f2e 84Revision 1.25 2000/11/02 07:25:31 kowal2
85Changes due to the new hit structure.
86Memory leak removed.
87
39c8eb58 88Revision 1.24 2000/10/05 16:06:09 kowal2
89Forward declarations. Changes due to a new class AliComplexCluster.
90
2b06d5c3 91Revision 1.23 2000/10/02 21:28:18 fca
92Removal of useless dependecies via forward declarations
93
94de3818 94Revision 1.22 2000/07/10 20:57:39 hristov
95Update of TPC code and macros by M.Kowalski
96
73042f01 97Revision 1.19.2.4 2000/06/26 07:39:42 kowal2
98Changes to obey the coding rules
99
100Revision 1.19.2.3 2000/06/25 08:38:41 kowal2
101Splitted from AliTPCtracking
102
103Revision 1.19.2.2 2000/06/16 12:59:28 kowal2
104Changed parameter settings
105
106Revision 1.19.2.1 2000/06/09 07:15:07 kowal2
107
108Defaults loaded automatically (hard-wired)
109Optional parameters can be set via macro called in the constructor
110
111Revision 1.19 2000/04/18 19:00:59 fca
112Small bug fixes to TPC files
113
4d68a14a 114Revision 1.18 2000/04/17 09:37:33 kowal2
115removed obsolete AliTPCDigitsDisplay.C
116
cc80f89e 117Revision 1.17.2.2 2000/04/10 08:15:12 kowal2
118
119New, experimental data structure from M. Ivanov
120New tracking algorithm
121Different pad geometry for different sectors
122Digitization rewritten
123
124Revision 1.17.2.1 2000/04/10 07:56:53 kowal2
125Not used anymore - removed
126
127Revision 1.17 2000/01/19 17:17:30 fca
128Introducing a list of lists of hits -- more hits allowed for detector now
129
1cedd08a 130Revision 1.16 1999/11/05 09:29:23 fca
131Accept only signals > 0
132
921bf71a 133Revision 1.15 1999/10/08 06:26:53 fca
134Removed ClustersIndex - not used anymore
135
e674b993 136Revision 1.14 1999/09/29 09:24:33 fca
137Introduction of the Copyright and cvs Log
138
4c039060 139*/
140
fe4da5cc 141///////////////////////////////////////////////////////////////////////////////
142// //
143// Time Projection Chamber //
144// This class contains the basic functions for the Time Projection Chamber //
145// detector. Functions specific to one particular geometry are //
146// contained in the derived classes //
147// //
148//Begin_Html
149/*
1439f98e 150<img src="picts/AliTPCClass.gif">
fe4da5cc 151*/
152//End_Html
153// //
8c555625 154// //
fe4da5cc 155///////////////////////////////////////////////////////////////////////////////
156
73042f01 157//
158
fe4da5cc 159#include <TMath.h>
160#include <TRandom.h>
161#include <TVector.h>
8c555625 162#include <TMatrix.h>
fe4da5cc 163#include <TGeometry.h>
164#include <TNode.h>
165#include <TTUBS.h>
166#include <TObjectTable.h>
1578254f 167#include "TParticle.h"
fe4da5cc 168#include "AliTPC.h"
afc42102 169#include <TFile.h>
170#include <TROOT.h>
171#include <TSystem.h>
fe4da5cc 172#include "AliRun.h"
173#include <iostream.h>
2ab0c725 174#include <stdlib.h>
fe4da5cc 175#include <fstream.h>
176#include "AliMC.h"
94de3818 177#include "AliMagF.h"
fe4da5cc 178
cc80f89e 179
73042f01 180#include "AliTPCParamSR.h"
8c555625 181#include "AliTPCPRF2D.h"
182#include "AliTPCRF1D.h"
cc80f89e 183#include "AliDigits.h"
184#include "AliSimDigits.h"
39c8eb58 185#include "AliTPCTrackHits.h"
186#include "AliPoints.h"
187#include "AliArrayBranch.h"
188
cc80f89e 189
190#include "AliTPCDigitsArray.h"
2b06d5c3 191#include "AliComplexCluster.h"
cc80f89e 192#include "AliClusters.h"
193#include "AliTPCClustersRow.h"
194#include "AliTPCClustersArray.h"
8c555625 195
73042f01 196#include "AliTPCcluster.h"
197#include "AliTPCclusterer.h"
198#include "AliTPCtracker.h"
8c555625 199
73042f01 200#include <TInterpreter.h>
2b06d5c3 201#include <TTree.h>
1283eee5 202
39c8eb58 203
204
fe4da5cc 205ClassImp(AliTPC)
206
207//_____________________________________________________________________________
208AliTPC::AliTPC()
209{
210 //
211 // Default constructor
212 //
213 fIshunt = 0;
fe4da5cc 214 fHits = 0;
215 fDigits = 0;
fe4da5cc 216 fNsectors = 0;
cc80f89e 217 //MI changes
218 fDigitsArray = 0;
219 fClustersArray = 0;
afc42102 220 fDefaults = 0;
39c8eb58 221 fTrackHits = 0;
222 fHitType = 2;
223 fTPCParam = 0;
fe4da5cc 224}
225
226//_____________________________________________________________________________
227AliTPC::AliTPC(const char *name, const char *title)
228 : AliDetector(name,title)
229{
230 //
231 // Standard constructor
232 //
233
234 //
235 // Initialise arrays of hits and digits
236 fHits = new TClonesArray("AliTPChit", 176);
4d68a14a 237 gAlice->AddHitList(fHits);
cc80f89e 238 //MI change
239 fDigitsArray = 0;
240 fClustersArray= 0;
afc42102 241 fDefaults = 0;
fe4da5cc 242 //
39c8eb58 243 fTrackHits = new AliTPCTrackHits; //MI - 13.09.2000
244 fTrackHits->SetHitPrecision(0.002);
245 fTrackHits->SetStepPrecision(0.003);
246 fTrackHits->SetMaxDistance(100);
247
248 fHitType = 2;
249 //
fe4da5cc 250 // Initialise counters
cc80f89e 251 fNsectors = 0;
cc80f89e 252
fe4da5cc 253 //
254 fIshunt = 0;
255 //
256 // Initialise color attributes
257 SetMarkerColor(kYellow);
73042f01 258
259 //
260 // Set TPC parameters
261 //
262
afc42102 263
264 if (!strcmp(title,"Default")) {
265 fTPCParam = new AliTPCParamSR;
73042f01 266 } else {
267 cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
268 fTPCParam=0;
269 }
270
fe4da5cc 271}
272
273//_____________________________________________________________________________
274AliTPC::~AliTPC()
275{
276 //
277 // TPC destructor
278 //
73042f01 279
fe4da5cc 280 fIshunt = 0;
281 delete fHits;
282 delete fDigits;
73042f01 283 delete fTPCParam;
39c8eb58 284 delete fTrackHits; //MI 15.09.2000
fe4da5cc 285}
286
fe4da5cc 287//_____________________________________________________________________________
288void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
289{
290 //
291 // Add a hit to the list
292 //
39c8eb58 293 // TClonesArray &lhits = *fHits;
294 // new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
295 if (fHitType&1){
296 TClonesArray &lhits = *fHits;
297 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
298 }
299 if (fHitType&2)
300 AddHit2(track,vol,hits);
fe4da5cc 301}
302
fe4da5cc 303//_____________________________________________________________________________
304void AliTPC::BuildGeometry()
305{
cc80f89e 306
fe4da5cc 307 //
308 // Build TPC ROOT TNode geometry for the event display
309 //
73042f01 310 TNode *nNode, *nTop;
fe4da5cc 311 TTUBS *tubs;
312 Int_t i;
313 const int kColorTPC=19;
1283eee5 314 char name[5], title[25];
fe4da5cc 315 const Double_t kDegrad=TMath::Pi()/180;
1283eee5 316 const Double_t kRaddeg=180./TMath::Pi();
317
1283eee5 318
73042f01 319 Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
320 Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
1283eee5 321
73042f01 322 Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
323 Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
1283eee5 324
325 Int_t nLo = fTPCParam->GetNInnerSector()/2;
326 Int_t nHi = fTPCParam->GetNOuterSector()/2;
327
73042f01 328 const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
329 const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
330 const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
331 const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
1283eee5 332
333
73042f01 334 const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
335 const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
1283eee5 336
337 Double_t rl,ru;
338
339
fe4da5cc 340 //
341 // Get ALICE top node
fe4da5cc 342 //
1283eee5 343
73042f01 344 nTop=gAlice->GetGeometry()->GetNode("alice");
1283eee5 345
346 // inner sectors
347
cc80f89e 348 rl = fTPCParam->GetInnerRadiusLow();
349 ru = fTPCParam->GetInnerRadiusUp();
1283eee5 350
351
fe4da5cc 352 for(i=0;i<nLo;i++) {
353 sprintf(name,"LS%2.2d",i);
1283eee5 354 name[4]='\0';
355 sprintf(title,"TPC low sector %3d",i);
356 title[24]='\0';
357
73042f01 358 tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
359 kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
fe4da5cc 360 tubs->SetNumberOfDivisions(1);
73042f01 361 nTop->cd();
362 nNode = new TNode(name,title,name,0,0,0,"");
363 nNode->SetLineColor(kColorTPC);
364 fNodes->Add(nNode);
fe4da5cc 365 }
1283eee5 366
fe4da5cc 367 // Outer sectors
1283eee5 368
cc80f89e 369 rl = fTPCParam->GetOuterRadiusLow();
370 ru = fTPCParam->GetOuterRadiusUp();
1283eee5 371
fe4da5cc 372 for(i=0;i<nHi;i++) {
373 sprintf(name,"US%2.2d",i);
1283eee5 374 name[4]='\0';
fe4da5cc 375 sprintf(title,"TPC upper sector %d",i);
1283eee5 376 title[24]='\0';
73042f01 377 tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
378 khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
fe4da5cc 379 tubs->SetNumberOfDivisions(1);
73042f01 380 nTop->cd();
381 nNode = new TNode(name,title,name,0,0,0,"");
382 nNode->SetLineColor(kColorTPC);
383 fNodes->Add(nNode);
fe4da5cc 384 }
cc80f89e 385
73042f01 386}
1283eee5 387
fe4da5cc 388//_____________________________________________________________________________
389Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
390{
391 //
392 // Calculate distance from TPC to mouse on the display
393 // Dummy procedure
394 //
395 return 9999;
396}
397
73042f01 398void AliTPC::Clusters2Tracks(TFile *of) {
3c0f9266 399 //-----------------------------------------------------------------
400 // This is a track finder.
3c0f9266 401 //-----------------------------------------------------------------
b9de75e1 402 AliTPCtracker tracker(fTPCParam);
403 tracker.Clusters2Tracks(gFile,of);
fe4da5cc 404}
405
406//_____________________________________________________________________________
407void AliTPC::CreateMaterials()
408{
8c555625 409 //-----------------------------------------------
37831078 410 // Create Materials for for TPC simulations
8c555625 411 //-----------------------------------------------
412
413 //-----------------------------------------------------------------
414 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
415 //-----------------------------------------------------------------
1283eee5 416
73042f01 417 Int_t iSXFLD=gAlice->Field()->Integ();
418 Float_t sXMGMX=gAlice->Field()->Max();
1283eee5 419
420 Float_t amat[5]; // atomic numbers
421 Float_t zmat[5]; // z
422 Float_t wmat[5]; // proportions
423
424 Float_t density;
37831078 425 Float_t apure[2];
1283eee5 426
1283eee5 427
37831078 428 //***************** Gases *************************
429
430 //-------------------------------------------------
1283eee5 431 // pure gases
37831078 432 //-------------------------------------------------
1283eee5 433
37831078 434 // Neon
1283eee5 435
436
37831078 437 amat[0]= 20.18;
438 zmat[0]= 10.;
1283eee5 439 density = 0.0009;
37831078 440
441 apure[0]=amat[0];
1283eee5 442
37831078 443 AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
1283eee5 444
37831078 445 // Argon
1283eee5 446
37831078 447 amat[0]= 39.948;
448 zmat[0]= 18.;
449 density = 0.001782;
1283eee5 450
37831078 451 apure[1]=amat[0];
1283eee5 452
37831078 453 AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
454
1283eee5 455
456 //--------------------------------------------------------------
457 // gases - compounds
458 //--------------------------------------------------------------
459
460 Float_t amol[3];
461
37831078 462 // CO2
1283eee5 463
464 amat[0]=12.011;
465 amat[1]=15.9994;
466
467 zmat[0]=6.;
468 zmat[1]=8.;
469
470 wmat[0]=1.;
471 wmat[1]=2.;
472
473 density=0.001977;
474
475 amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
476
477 AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
37831078 478
1283eee5 479 // CF4
480
481 amat[0]=12.011;
482 amat[1]=18.998;
483
484 zmat[0]=6.;
485 zmat[1]=9.;
486
487 wmat[0]=1.;
488 wmat[1]=4.;
489
490 density=0.003034;
491
492 amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
493
494 AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
495
37831078 496
1283eee5 497 // CH4
498
499 amat[0]=12.011;
500 amat[1]=1.;
501
502 zmat[0]=6.;
503 zmat[1]=1.;
504
505 wmat[0]=1.;
506 wmat[1]=4.;
507
508 density=0.000717;
509
510 amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
511
512 AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
513
514 //----------------------------------------------------------------
515 // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
516 //----------------------------------------------------------------
37831078 517
518 char namate[21];
1283eee5 519 density = 0.;
520 Float_t am=0;
521 Int_t nc;
37831078 522 Float_t rho,absl,X0,buf[1];
1283eee5 523 Int_t nbuf;
37831078 524 Float_t a,z;
1283eee5 525
526 for(nc = 0;nc<fNoComp;nc++)
527 {
528
529 // retrive material constants
530
37831078 531 gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
1283eee5 532
533 amat[nc] = a;
534 zmat[nc] = z;
535
536 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
537
37831078 538 am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
1283eee5 539 density += fMixtProp[nc]*rho; // density of the mixture
540
541 }
542
543 // mixture proportions by weight!
544
545 for(nc = 0;nc<fNoComp;nc++)
546 {
547
548 Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
549
37831078 550 wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
551 apure[nnc] : amol[nnc])/am;
552
553 }
554
555 // Drift gases 1 - nonsensitive, 2 - sensitive
1283eee5 556
1283eee5 557 AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
558 AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
1283eee5 559
1283eee5 560
37831078 561 // Air
562
563 amat[0] = 14.61;
564 zmat[0] = 7.3;
565 density = 0.001205;
1283eee5 566
37831078 567 AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
1283eee5 568
1283eee5 569
570 //----------------------------------------------------------------------
571 // solid materials
572 //----------------------------------------------------------------------
573
1283eee5 574
37831078 575 // Kevlar C14H22O2N2
1283eee5 576
37831078 577 amat[0] = 12.011;
578 amat[1] = 1.;
579 amat[2] = 15.999;
580 amat[3] = 14.006;
1283eee5 581
37831078 582 zmat[0] = 6.;
583 zmat[1] = 1.;
584 zmat[2] = 8.;
585 zmat[3] = 7.;
586
587 wmat[0] = 14.;
588 wmat[1] = 22.;
589 wmat[2] = 2.;
590 wmat[3] = 2.;
591
592 density = 1.45;
593
594 AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
595
596 // NOMEX
1283eee5 597
37831078 598 amat[0] = 12.011;
599 amat[1] = 1.;
600 amat[2] = 15.999;
601 amat[3] = 14.006;
602
603 zmat[0] = 6.;
604 zmat[1] = 1.;
605 zmat[2] = 8.;
606 zmat[3] = 7.;
607
608 wmat[0] = 14.;
609 wmat[1] = 22.;
610 wmat[2] = 2.;
611 wmat[3] = 2.;
612
613 density = 0.03;
1283eee5 614
fe4da5cc 615
37831078 616 AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
617
618 // Makrolon C16H18O3
619
620 amat[0] = 12.011;
621 amat[1] = 1.;
622 amat[2] = 15.999;
1283eee5 623
37831078 624 zmat[0] = 6.;
625 zmat[1] = 1.;
626 zmat[2] = 8.;
627
628 wmat[0] = 16.;
629 wmat[1] = 18.;
630 wmat[2] = 3.;
631
632 density = 1.2;
633
634 AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
635
1283eee5 636 // Mylar C5H4O2
637
638 amat[0]=12.011;
639 amat[1]=1.;
640 amat[2]=15.9994;
641
642 zmat[0]=6.;
643 zmat[1]=1.;
644 zmat[2]=8.;
645
646 wmat[0]=5.;
647 wmat[1]=4.;
648 wmat[2]=2.;
649
650 density = 1.39;
fe4da5cc 651
37831078 652 AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
1283eee5 653
ba1d05f9 654 // SiO2 - used later for the glass fiber
1283eee5 655
37831078 656 amat[0]=28.086;
657 amat[1]=15.9994;
1283eee5 658
37831078 659 zmat[0]=14.;
660 zmat[1]=8.;
1283eee5 661
37831078 662 wmat[0]=1.;
663 wmat[1]=2.;
1283eee5 664
1283eee5 665
ba1d05f9 666 AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2)
1283eee5 667
37831078 668 // Al
1283eee5 669
37831078 670 amat[0] = 26.98;
671 zmat[0] = 13.;
1283eee5 672
37831078 673 density = 2.7;
1283eee5 674
37831078 675 AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
1283eee5 676
37831078 677 // Si
1283eee5 678
37831078 679 amat[0] = 28.086;
9a3667bd 680 zmat[0] = 14.;
1283eee5 681
37831078 682 density = 2.33;
1283eee5 683
37831078 684 AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
1283eee5 685
37831078 686 // Cu
1283eee5 687
37831078 688 amat[0] = 63.546;
689 zmat[0] = 29.;
1283eee5 690
37831078 691 density = 8.96;
1283eee5 692
37831078 693 AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
1283eee5 694
37831078 695 // Tedlar C2H3F
1283eee5 696
37831078 697 amat[0] = 12.011;
698 amat[1] = 1.;
699 amat[2] = 18.998;
1283eee5 700
37831078 701 zmat[0] = 6.;
702 zmat[1] = 1.;
703 zmat[2] = 9.;
1283eee5 704
37831078 705 wmat[0] = 2.;
706 wmat[1] = 3.;
707 wmat[2] = 1.;
1283eee5 708
37831078 709 density = 1.71;
1283eee5 710
37831078 711 AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
1283eee5 712
1283eee5 713
37831078 714 // Plexiglas C5H8O2
1283eee5 715
37831078 716 amat[0]=12.011;
717 amat[1]=1.;
718 amat[2]=15.9994;
1283eee5 719
37831078 720 zmat[0]=6.;
721 zmat[1]=1.;
722 zmat[2]=8.;
1283eee5 723
37831078 724 wmat[0]=5.;
725 wmat[1]=8.;
726 wmat[2]=2.;
1283eee5 727
37831078 728 density=1.18;
1283eee5 729
37831078 730 AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
1283eee5 731
ba1d05f9 732 // Epoxy - C14 H20 O3
1283eee5 733
5a28a08f 734
735 amat[0]=12.011;
736 amat[1]=1.;
737 amat[2]=15.9994;
738
739 zmat[0]=6.;
740 zmat[1]=1.;
741 zmat[2]=8.;
742
ba1d05f9 743 wmat[0]=14.;
744 wmat[1]=20.;
5a28a08f 745 wmat[2]=3.;
746
ba1d05f9 747 density=1.25;
5a28a08f 748
749 AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat);
37831078 750
ba1d05f9 751 // Carbon
752
753 amat[0]=12.011;
754 zmat[0]=6.;
755 density= 2.265;
756
757 AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.);
758
759 // get epoxy
760
761 gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf);
762
763 // Carbon fiber
764
765 wmat[0]=0.644; // by weight!
766 wmat[1]=0.356;
767
768 density=0.5*(1.25+2.265);
769
770 AliMixture(47,"Cfiber",amat,zmat,density,2,wmat);
771
772 // get SiO2
773
774 gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf);
775
776 wmat[0]=0.725; // by weight!
777 wmat[1]=0.275;
778
779 density=1.7;
780
781 AliMixture(39,"G10",amat,zmat,density,2,wmat);
782
783
784
785
37831078 786 //----------------------------------------------------------
787 // tracking media for gases
788 //----------------------------------------------------------
789
790 AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
791 AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
792 AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
793 AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
794
795 //-----------------------------------------------------------
796 // tracking media for solids
797 //-----------------------------------------------------------
798
799 AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
800 AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
801 AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
802 AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
803 AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
804 AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
805 AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
806 AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
807 AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
808 AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
5a28a08f 809 AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
ba1d05f9 810 AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
1283eee5 811
fe4da5cc 812}
813
fe4da5cc 814
afc42102 815void AliTPC::Digits2Clusters(TFile *of, Int_t eventnumber)
fe4da5cc 816{
3c0f9266 817 //-----------------------------------------------------------------
818 // This is a simple cluster finder.
3c0f9266 819 //-----------------------------------------------------------------
afc42102 820 AliTPCclusterer::Digits2Clusters(fTPCParam,of,eventnumber);
fe4da5cc 821}
822
73042f01 823extern Double_t SigmaY2(Double_t, Double_t, Double_t);
824extern Double_t SigmaZ2(Double_t, Double_t);
fe4da5cc 825//_____________________________________________________________________________
afc42102 826void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
fe4da5cc 827{
8c555625 828 //--------------------------------------------------------
fe4da5cc 829 // TPC simple cluster generator from hits
830 // obtained from the TPC Fast Simulator
8c555625 831 // The point errors are taken from the parametrization
832 //--------------------------------------------------------
833
834 //-----------------------------------------------------------------
835 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
836 //-----------------------------------------------------------------
73042f01 837 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
838 // Jouri.Belikov@cern.ch
839 //----------------------------------------------------------------
840
841 /////////////////////////////////////////////////////////////////////////////
842 //
843 //---------------------------------------------------------------------
844 // ALICE TPC Cluster Parameters
845 //--------------------------------------------------------------------
846
847
848
849 // Cluster width in rphi
850 const Float_t kACrphi=0.18322;
851 const Float_t kBCrphi=0.59551e-3;
852 const Float_t kCCrphi=0.60952e-1;
853 // Cluster width in z
854 const Float_t kACz=0.19081;
855 const Float_t kBCz=0.55938e-3;
856 const Float_t kCCz=0.30428;
857
858 TDirectory *savedir=gDirectory;
859
860 if (!of->IsOpen()) {
861 cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
862 return;
863 }
3c0f9266 864
afc42102 865 //if(fDefaults == 0) SetDefaults();
cc80f89e 866
73042f01 867 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
fe4da5cc 868 //
1578254f 869 TParticle *particle; // pointer to a given particle
fe4da5cc 870 AliTPChit *tpcHit; // pointer to a sigle TPC hit
39c8eb58 871 Int_t sector;
fe4da5cc 872 Int_t ipart;
873 Float_t xyz[5];
874 Float_t pl,pt,tanth,rpad,ratio;
fe4da5cc 875 Float_t cph,sph;
876
877 //---------------------------------------------------------------
878 // Get the access to the tracks
879 //---------------------------------------------------------------
880
73042f01 881 TTree *tH = gAlice->TreeH();
882 Stat_t ntracks = tH->GetEntries();
73042f01 883
884 //Switch to the output file
885 of->cd();
886
afc42102 887 char cname[100];
888
889 sprintf(cname,"TreeC_TPC_%d",eventn);
890
73042f01 891 fTPCParam->Write(fTPCParam->GetTitle());
892 AliTPCClustersArray carray;
893 carray.Setup(fTPCParam);
894 carray.SetClusterType("AliTPCcluster");
895 carray.MakeTree();
896
897 Int_t nclusters=0; //cluster counter
fe4da5cc 898
899 //------------------------------------------------------------
cc80f89e 900 // Loop over all sectors (72 sectors for 20 deg
901 // segmentation for both lower and upper sectors)
902 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
3c0f9266 903 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
fe4da5cc 904 //
3c0f9266 905 // First cluster for sector 0 starts at "0"
fe4da5cc 906 //------------------------------------------------------------
cc80f89e 907
908 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
8c555625 909 //MI change
cc80f89e 910 fTPCParam->AdjustCosSin(isec,cph,sph);
fe4da5cc 911
912 //------------------------------------------------------------
913 // Loop over tracks
914 //------------------------------------------------------------
915
916 for(Int_t track=0;track<ntracks;track++){
917 ResetHits();
73042f01 918 tH->GetEvent(track);
fe4da5cc 919 //
73042f01 920 // Get number of the TPC hits
fe4da5cc 921 //
39c8eb58 922 // nhits=fHits->GetEntriesFast();
fe4da5cc 923 //
39c8eb58 924
925 tpcHit = (AliTPChit*)FirstHit(-1);
926
fe4da5cc 927 // Loop over hits
928 //
39c8eb58 929 // for(Int_t hit=0;hit<nhits;hit++){
930 //tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
931
932 while(tpcHit){
933
934 if (tpcHit->fQ == 0.) {
935 tpcHit = (AliTPChit*) NextHit();
936 continue; //information about track (I.Belikov)
937 }
fe4da5cc 938 sector=tpcHit->fSector; // sector number
39c8eb58 939
940
941 // if(sector != isec) continue; //terminate iteration
942
943 if(sector != isec){
944 tpcHit = (AliTPChit*) NextHit();
945 continue;
946 }
94de3818 947 ipart=tpcHit->Track();
2ab0c725 948 particle=gAlice->Particle(ipart);
1578254f 949 pl=particle->Pz();
950 pt=particle->Pt();
fe4da5cc 951 if(pt < 1.e-9) pt=1.e-9;
952 tanth=pl/pt;
953 tanth = TMath::Abs(tanth);
94de3818 954 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
fe4da5cc 955 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
73042f01 956
fe4da5cc 957 // space-point resolutions
958
73042f01 959 sigmaRphi=SigmaY2(rpad,tanth,pt);
960 sigmaZ =SigmaZ2(rpad,tanth );
fe4da5cc 961
962 // cluster widths
963
73042f01 964 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
965 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
fe4da5cc 966
967 // temporary protection
968
73042f01 969 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
970 if(sigmaZ < 0.) sigmaZ=0.4e-3;
971 if(clRphi < 0.) clRphi=2.5e-3;
972 if(clZ < 0.) clZ=2.5e-5;
fe4da5cc 973
974 //
cc80f89e 975
976 //
977 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
fe4da5cc 978 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
979 //
94de3818 980 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
981 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
73042f01 982 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
983 Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
984 fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
985 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
986 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
94de3818 987 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
988 if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z();
39c8eb58 989 xyz[2]=sigmaRphi; // fSigmaY2
990 xyz[3]=sigmaZ; // fSigmaZ2
991 xyz[4]=tpcHit->fQ; // q
73042f01 992
993 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
994 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
995
94de3818 996 Int_t tracks[3]={tpcHit->Track(), -1, -1};
39c8eb58 997 AliTPCcluster cluster(tracks,xyz);
73042f01 998
999 clrow->InsertCluster(&cluster); nclusters++;
1000
39c8eb58 1001 tpcHit = (AliTPChit*)NextHit();
1002
1003
fe4da5cc 1004 } // end of loop over hits
73042f01 1005
1006 } // end of loop over tracks
1007
1008 Int_t nrows=fTPCParam->GetNRow(isec);
1009 for (Int_t irow=0; irow<nrows; irow++) {
1010 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
1011 if (!clrow) continue;
1012 carray.StoreRow(isec,irow);
1013 carray.ClearRow(isec,irow);
1014 }
1015
cc80f89e 1016 } // end of loop over sectors
73042f01 1017
1018 cerr<<"Number of made clusters : "<<nclusters<<" \n";
1019
6bc3458d 1020 carray.GetTree()->SetName(cname);
1021 carray.GetTree()->Write();
73042f01 1022
1023 savedir->cd(); //switch back to the input file
cc80f89e 1024
1025} // end of function
1026
1027//_________________________________________________________________
1028void AliTPC::Hits2ExactClustersSector(Int_t isec)
1029{
1030 //--------------------------------------------------------
1031 //calculate exact cross point of track and given pad row
1032 //resulting values are expressed in "digit" coordinata
1033 //--------------------------------------------------------
1034
1035 //-----------------------------------------------------------------
1036 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1037 //-----------------------------------------------------------------
1038 //
1039 if (fClustersArray==0){
1040 return;
1041 }
1042 //
1043 TParticle *particle; // pointer to a given particle
1044 AliTPChit *tpcHit; // pointer to a sigle TPC hit
5f16d237 1045 // Int_t sector,nhits;
cc80f89e 1046 Int_t ipart;
73042f01 1047 const Int_t kcmaxhits=30000;
1048 TVector * xxxx = new TVector(kcmaxhits*4);
cc80f89e 1049 TVector & xxx = *xxxx;
73042f01 1050 Int_t maxhits = kcmaxhits;
cc80f89e 1051 //construct array for each padrow
1052 for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++)
1053 fClustersArray->CreateRow(isec,i);
fe4da5cc 1054
cc80f89e 1055 //---------------------------------------------------------------
1056 // Get the access to the tracks
1057 //---------------------------------------------------------------
1058
73042f01 1059 TTree *tH = gAlice->TreeH();
1060 Stat_t ntracks = tH->GetEntries();
2ab0c725 1061 Int_t npart = gAlice->GetNtrack();
5f16d237 1062 TBranch * br = tH->GetBranch("fTrackHitsInfo");
1063 //MI change
1064 TBranch * branch=0;
1065 if (fHitType&2) branch = tH->GetBranch("TPC2");
1066 else branch = tH->GetBranch("TPC");
1067
cc80f89e 1068 //------------------------------------------------------------
1069 // Loop over tracks
1070 //------------------------------------------------------------
fe4da5cc 1071
5f16d237 1072 for(Int_t track=0;track<ntracks;track++){
1073 Bool_t isInSector=kTRUE;
cc80f89e 1074 ResetHits();
5f16d237 1075
1076 if (fHitType&2) {
1077 isInSector=kFALSE;
1078 br->GetEvent(track);
1079 AliObjectArray * ar = fTrackHits->fTrackHitsInfo;
1080 for (UInt_t j=0;j<ar->GetSize();j++){
1081 if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE;
1082 }
1083 }
1084 if (!isInSector) continue;
1085 //MI change
1086 branch->GetEntry(track); // get next track
1087
1088
cc80f89e 1089 //
1090 // Get number of the TPC hits and a pointer
1091 // to the particles
1092 //
5f16d237 1093 //nhits=fHits->GetEntriesFast();
cc80f89e 1094 //
1095 // Loop over hits
1096 //
1097 Int_t currentIndex=0;
1098 Int_t lastrow=-1; //last writen row
5f16d237 1099
1100 //M.I. changes
1101
1102 tpcHit = (AliTPChit*)FirstHit(-1);
1103 while(tpcHit){
1104
1105 Int_t sector=tpcHit->fSector; // sector number
1106 if(sector != isec){
1107 tpcHit = (AliTPChit*) NextHit();
1108 continue;
1109 }
1110
1111 //for(Int_t hit=0;hit<nhits;hit++){
1112 //tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1113 //if (tpcHit==0) continue;
1114 //sector=tpcHit->fSector; // sector number
1115 //if(sector != isec) continue;
94de3818 1116 ipart=tpcHit->Track();
2ab0c725 1117 if (ipart<npart) particle=gAlice->Particle(ipart);
cc80f89e 1118
1119 //find row number
1120
94de3818 1121 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
cc80f89e 1122 Int_t index[3]={1,isec,0};
1123 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
5f16d237 1124 if (currentrow<0) {tpcHit = (AliTPChit*)NextHit(); continue;}
cc80f89e 1125 if (lastrow<0) lastrow=currentrow;
1126 if (currentrow==lastrow){
1127 if ( currentIndex>=maxhits){
73042f01 1128 maxhits+=kcmaxhits;
cc80f89e 1129 xxx.ResizeTo(4*maxhits);
1130 }
1131 xxx(currentIndex*4)=x[0];
1132 xxx(currentIndex*4+1)=x[1];
1133 xxx(currentIndex*4+2)=x[2];
1134 xxx(currentIndex*4+3)=tpcHit->fQ;
1135 currentIndex++;
1136 }
1137 else
1138 if (currentIndex>2){
1139 Float_t sumx=0;
1140 Float_t sumx2=0;
1141 Float_t sumx3=0;
1142 Float_t sumx4=0;
1143 Float_t sumy=0;
1144 Float_t sumxy=0;
1145 Float_t sumx2y=0;
1146 Float_t sumz=0;
1147 Float_t sumxz=0;
1148 Float_t sumx2z=0;
1149 Float_t sumq=0;
1150 for (Int_t index=0;index<currentIndex;index++){
1151 Float_t x,x2,x3,x4;
1152 x=x2=x3=x4=xxx(index*4);
1153 x2*=x;
1154 x3*=x2;
1155 x4*=x3;
1156 sumx+=x;
1157 sumx2+=x2;
1158 sumx3+=x3;
1159 sumx4+=x4;
1160 sumy+=xxx(index*4+1);
1161 sumxy+=xxx(index*4+1)*x;
1162 sumx2y+=xxx(index*4+1)*x2;
1163 sumz+=xxx(index*4+2);
1164 sumxz+=xxx(index*4+2)*x;
1165 sumx2z+=xxx(index*4+2)*x2;
1166 sumq+=xxx(index*4+3);
1167 }
73042f01 1168 Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
cc80f89e 1169 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
1170 sumx2*(sumx*sumx3-sumx2*sumx2);
1171
1172 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
1173 sumx2*(sumxy*sumx3-sumx2y*sumx2);
1174 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
1175 sumx2*(sumxz*sumx3-sumx2z*sumx2);
1176
1177 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
1178 sumx2*(sumx*sumx2y-sumx2*sumxy);
1179 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
1180 sumx2*(sumx*sumx2z-sumx2*sumxz);
1181
73042f01 1182 Float_t y=detay/det+centralPad;
cc80f89e 1183 Float_t z=detaz/det;
1184 Float_t by=detby/det; //y angle
1185 Float_t bz=detbz/det; //z angle
1186 sumy/=Float_t(currentIndex);
1187 sumz/=Float_t(currentIndex);
5f16d237 1188
cc80f89e 1189 AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
5f16d237 1190 if (row!=0) {
1191 AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
1192 // AliComplexCluster cl;
1193 cl->fX=z;
1194 cl->fY=y;
1195 cl->fQ=sumq;
1196 cl->fSigmaX2=bz;
1197 cl->fSigmaY2=by;
1198 cl->fTracks[0]=ipart;
1199 }
cc80f89e 1200 currentIndex=0;
1201 lastrow=currentrow;
1202 } //end of calculating cluster for given row
1203
1204
5f16d237 1205 tpcHit = (AliTPChit*)NextHit();
cc80f89e 1206 } // end of loop over hits
1207 } // end of loop over tracks
1208 //write padrows to tree
1209 for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1210 fClustersArray->StoreRow(isec,ii);
1211 fClustersArray->ClearRow(isec,ii);
1212 }
1213 xxxx->Delete();
1214
fe4da5cc 1215}
0a61bf9d 1216
1217
1218
1219//__
1220void AliTPC::SDigits2Digits2(Int_t eventnumber)
1221{
1222 //create digits from summable digits
1223
1224 char sname[100];
1225 char dname[100];
1226 sprintf(sname,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
1227 sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1228
1229 //conect tree with sSDigits
1230 TTree *t = (TTree *)gDirectory->Get(sname);
1231 AliSimDigits digarr, *dummy=&digarr;
1232 t->GetBranch("Segment")->SetAddress(&dummy);
1233 Stat_t nentries = t->GetEntries();
1234
5f16d237 1235 // set zero suppression
1236
1237 fTPCParam->SetZeroSup(2);
1238
1239 // get zero suppression
1240
1241 Int_t zerosup = fTPCParam->GetZeroSup();
1242
1243 //make tree with digits
1244
0a61bf9d 1245 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1246 arr->SetClass("AliSimDigits");
1247 arr->Setup(fTPCParam);
1248 arr->MakeTree(fDigitsFile);
1249
1250 AliTPCParam * par =fTPCParam;
5f16d237 1251
0a61bf9d 1252 //Loop over segments of the TPC
5f16d237 1253
0a61bf9d 1254 for (Int_t n=0; n<nentries; n++) {
1255 t->GetEvent(n);
1256 Int_t sec, row;
1257 if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1258 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
1259 continue;
1260 }
1261 AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1262 Int_t nrows = digrow->GetNRows();
1263 Int_t ncols = digrow->GetNCols();
1264
1265 digrow->ExpandBuffer();
1266 digarr.ExpandBuffer();
1267 digrow->ExpandTrackBuffer();
1268 digarr.ExpandTrackBuffer();
1269
5f16d237 1270
1271 for (Int_t rows=0;rows<nrows; rows++){
0a61bf9d 1272 for (Int_t col=0;col<ncols; col++){
5f16d237 1273 Float_t q = digarr.GetDigitFast(rows,col);
0a61bf9d 1274 //add noise
5f16d237 1275 q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()*16);
1276 q/=16; //conversion faktor
1277 q=(Short_t)q; // digits are short integers
1278
1279 if (q > zerosup){
1280 if(q > fTPCParam->GetADCSat()) q = (Short_t)(fTPCParam->GetADCSat());
1281 digrow->SetDigitFast(q,rows,col);
1282 for (Int_t tr=0;tr<3;tr++)
1283
1284((AliSimDigits*)digrow)->SetTrackIDFast(digarr.GetTrackIDFast(rows,col,tr)-2,rows,col,tr);
0a61bf9d 1285 }
5f16d237 1286 }
1287 }
1288 arr->StoreRow(sec,row);
1289 arr->ClearRow(sec,row);
9fe90376 1290 // cerr<<sec<<"\t"<<row<<"\n";
5f16d237 1291 }
0a61bf9d 1292
0a61bf9d 1293
5f16d237 1294 //write results
0a61bf9d 1295
5f16d237 1296
1297 arr->GetTree()->SetName(dname);
1298 arr->GetTree()->Write();
1299 delete arr;
1300}
1301//_________________________________________
1302void AliTPC::Merge(TTree * intree, Int_t *mask, Int_t nin, Int_t outid)
1303{
1304
1305 //intree - pointer to array of trees with s digits
1306 //mask - mask for each
1307 //nin - number of inputs
1308 //outid - event number of the output event
1309
1310 //create digits from summable digits
1311 //conect tree with sSDigits
1312
1313
1314 AliSimDigits ** digarr = new AliSimDigits*[nin];
1315 for (Int_t i1=0;i1<nin; i1++){
1316 digarr[i1]=0;
1317 intree[i1].GetBranch("Segment")->SetAddress(&digarr[i1]);
0a61bf9d 1318 }
5f16d237 1319 Stat_t nentries = intree[0].GetEntries();
1320
1321 //make tree with digits
1322 char dname[100];
1323 sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),outid);
1324 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1325 arr->SetClass("AliSimDigits");
1326 arr->Setup(fTPCParam);
1327 arr->MakeTree(fDigitsFile);
0a61bf9d 1328
5f16d237 1329 // set zero suppression
0a61bf9d 1330
5f16d237 1331 fTPCParam->SetZeroSup(2);
0a61bf9d 1332
5f16d237 1333 // get zero suppression
1334
1335 Int_t zerosup = fTPCParam->GetZeroSup();
1336
1337
1338 AliTPCParam * par =fTPCParam;
1339
1340 //Loop over segments of the TPC
1341 for (Int_t n=0; n<nentries; n++) {
0a61bf9d 1342
5f16d237 1343 for (Int_t i=0;i<nin; i++){
1344 //connect proper digits
1345 intree[i].GetEvent(n);
1346 digarr[i]->ExpandBuffer();
1347 digarr[i]->ExpandTrackBuffer();
1348 }
1349 Int_t sec, row;
1350 if (!par->AdjustSectorRow(digarr[0]->GetID(),sec,row)) {
1351 cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
1352 continue;
1353 }
0a61bf9d 1354
5f16d237 1355 AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1356 Int_t nrows = digrow->GetNRows();
1357 Int_t ncols = digrow->GetNCols();
0a61bf9d 1358
5f16d237 1359 digrow->ExpandBuffer();
1360 digrow->ExpandTrackBuffer();
1361
1362 for (Int_t rows=0;rows<nrows; rows++){
1363 for (Int_t col=0;col<ncols; col++){
1364 Float_t q=0;
1365 Int_t label[1000]; // i hope no more than 300 events merged
1366 Int_t labptr = 0;
1367 // looop over digits
1368 for (Int_t i=0;i<nin; i++){
1369 q += digarr[i]->GetDigitFast(rows,col);
1370 for (Int_t tr=0;tr<3;tr++) {
1371 Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
1372 if ( lab > 1) {
1373 label[labptr]=lab+mask[i]-2;
1374 labptr++;
1375 }
1376 }
1377 }
1378 //add noise
1379 q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()*16);
1380
1381 q/=16; //conversion faktor
1382 q=(Short_t)q;
1383
1384 if (q> zerosup){
1385
1386 if(q > fTPCParam->GetADCSat()) q = (Short_t)(fTPCParam->GetADCSat());
1387 digrow->SetDigitFast(q,rows,col);
1388 for (Int_t tr=0;tr<3;tr++){
1389 if (tr<labptr)
1390 ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],rows,col,tr);
1391 else
1392 ((AliSimDigits*)digrow)->SetTrackIDFast(-1,rows,col,tr);
1393 }
1394 }
1395 }
1396 }
1397 arr->StoreRow(sec,row);
1398
1399 arr->ClearRow(sec,row);
9fe90376 1400 // cerr<<sec<<"\t"<<row<<"\n";
5f16d237 1401 }
1402
1403 delete digarr;
1404 arr->GetTree()->SetName(dname);
1405 arr->GetTree()->Write();
1406
1407 delete arr;
0a61bf9d 1408
0a61bf9d 1409}
1410
0a61bf9d 1411/*_________________________________________
afc42102 1412void AliTPC::SDigits2Digits(Int_t eventnumber)
2ab0c725 1413{
afc42102 1414
1415
1416 cerr<<"Digitizing TPC...\n";
1417
1418 Hits2Digits(eventnumber);
1419
1420
1421 //write results
1422
1423 // char treeName[100];
1424
1425 // sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1426
1427 // GetDigitsArray()->GetTree()->Write(treeName);
1428}
0a61bf9d 1429*/
afc42102 1430//__________________________________________________________________
1431void AliTPC::SetDefaults(){
1432
1433
1434 cerr<<"Setting default parameters...\n";
1435
1436 // Set response functions
1437
2ab0c725 1438 AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1439 AliTPCPRF2D * prfinner = new AliTPCPRF2D;
1440 AliTPCPRF2D * prfouter = new AliTPCPRF2D;
1441 AliTPCRF1D * rf = new AliTPCRF1D(kTRUE);
2ab0c725 1442 rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1443 rf->SetOffset(3*param->GetZSigma());
1444 rf->Update();
afc42102 1445
1446 TDirectory *savedir=gDirectory;
2ab0c725 1447 TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1448 if (!f->IsOpen()) {
afc42102 1449 cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ;
2ab0c725 1450 exit(3);
1451 }
1452 prfinner->Read("prf_07504_Gati_056068_d02");
1453 prfouter->Read("prf_10006_Gati_047051_d03");
1454 f->Close();
afc42102 1455 savedir->cd();
1456
2ab0c725 1457 param->SetInnerPRF(prfinner);
1458 param->SetOuterPRF(prfouter);
1459 param->SetTimeRF(rf);
1460
afc42102 1461 // set fTPCParam
1462
2ab0c725 1463 SetParam(param);
1464
afc42102 1465
1466 fDefaults = 1;
1467
1468}
1469//__________________________________________________________________
1470void AliTPC::Hits2Digits(Int_t eventnumber)
1471{
1472 //----------------------------------------------------
1473 // Loop over all sectors for a single event
1474 //----------------------------------------------------
1475
1476
1477 if(fDefaults == 0) SetDefaults(); // check if the parameters are set
2ab0c725 1478
1479 //setup TPCDigitsArray
afc42102 1480
1481 if(GetDigitsArray()) delete GetDigitsArray();
1482
2ab0c725 1483 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1484 arr->SetClass("AliSimDigits");
afc42102 1485 arr->Setup(fTPCParam);
2ab0c725 1486 arr->MakeTree(fDigitsFile);
2ab0c725 1487 SetDigitsArray(arr);
1488
afc42102 1489 fDigitsSwitch=0; // standard digits
2ab0c725 1490
afc42102 1491 cerr<<"Digitizing TPC...\n";
fe4da5cc 1492
afc42102 1493 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
cc80f89e 1494
afc42102 1495 // write results
1496
1497 char treeName[100];
1498
1499 sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
f8cf550c 1500
6bc3458d 1501 GetDigitsArray()->GetTree()->SetName(treeName);
1502 GetDigitsArray()->GetTree()->Write();
cc80f89e 1503
cc80f89e 1504
8c555625 1505}
1506
afc42102 1507
1508
f8cf550c 1509//__________________________________________________________________
0a61bf9d 1510void AliTPC::Hits2SDigits2(Int_t eventnumber)
f8cf550c 1511{
1512
1513 //-----------------------------------------------------------
1514 // summable digits - 16 bit "ADC", no noise, no saturation
1515 //-----------------------------------------------------------
1516
1517 //----------------------------------------------------
afc42102 1518 // Loop over all sectors for a single event
f8cf550c 1519 //----------------------------------------------------
1520
afc42102 1521
1522 if(fDefaults == 0) SetDefaults();
1523
1524 //setup TPCDigitsArray
1525
1526 if(GetDigitsArray()) delete GetDigitsArray();
1527
1528 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1529 arr->SetClass("AliSimDigits");
1530 arr->Setup(fTPCParam);
1531 arr->MakeTree(fDigitsFile);
1532 SetDigitsArray(arr);
1533
1534 cerr<<"Digitizing TPC...\n";
1535
1536 fDigitsSwitch=1; // summable digits
5f16d237 1537
1538 // set zero suppression to "0"
1539
1540 fTPCParam->SetZeroSup(0);
f8cf550c 1541
1542 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1543
afc42102 1544
1545 // write results
1546
1547 char treeName[100];
1548
1549 sprintf(treeName,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
1550
6bc3458d 1551 GetDigitsArray()->GetTree()->SetName(treeName);
1552 GetDigitsArray()->GetTree()->Write();
afc42102 1553
f8cf550c 1554}
8c555625 1555
afc42102 1556
1557
0a61bf9d 1558//__________________________________________________________________
1559void AliTPC::Hits2SDigits()
1560{
1561
1562 //-----------------------------------------------------------
1563 // summable digits - 16 bit "ADC", no noise, no saturation
1564 //-----------------------------------------------------------
1565
1566 //----------------------------------------------------
1567 // Loop over all sectors for a single event
1568 //----------------------------------------------------
1569 //MI change - for pp run
1570 Int_t eventnumber = gAlice->GetEvNumber();
1571
1572 if(fDefaults == 0) SetDefaults();
1573
1574 //setup TPCDigitsArray
1575
1576 if(GetDigitsArray()) delete GetDigitsArray();
1577
1578 AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1579 arr->SetClass("AliSimDigits");
1580 arr->Setup(fTPCParam);
1581 arr->MakeTree(fDigitsFile);
1582 SetDigitsArray(arr);
1583
1584 cerr<<"Digitizing TPC...\n";
1585
1586 // fDigitsSwitch=1; // summable digits -for the moment direct
1587
1588 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1589
1590
1591 // write results
1592 char treeName[100];
1593
1594 sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1595
6bc3458d 1596 GetDigitsArray()->GetTree()->SetName(treeName);
1597 GetDigitsArray()->GetTree()->Write();
0a61bf9d 1598
1599}
1600
1601
fe4da5cc 1602//_____________________________________________________________________________
8c555625 1603void AliTPC::Hits2DigitsSector(Int_t isec)
fe4da5cc 1604{
8c555625 1605 //-------------------------------------------------------------------
fe4da5cc 1606 // TPC conversion from hits to digits.
8c555625 1607 //-------------------------------------------------------------------
1608
1609 //-----------------------------------------------------------------
1610 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1611 //-----------------------------------------------------------------
1612
fe4da5cc 1613 //-------------------------------------------------------
8c555625 1614 // Get the access to the track hits
fe4da5cc 1615 //-------------------------------------------------------
8c555625 1616
afc42102 1617 // check if the parameters are set - important if one calls this method
1618 // directly, not from the Hits2Digits
1619
1620 if(fDefaults == 0) SetDefaults();
cc80f89e 1621
73042f01 1622 TTree *tH = gAlice->TreeH(); // pointer to the hits tree
1623 Stat_t ntracks = tH->GetEntries();
8c555625 1624
1625 if( ntracks > 0){
1626
1627 //-------------------------------------------
1628 // Only if there are any tracks...
1629 //-------------------------------------------
1630
8c555625 1631 TObjArray **row;
fe4da5cc 1632
a62d28a2 1633 //printf("*** Processing sector number %d ***\n",isec);
8c555625 1634
1635 Int_t nrows =fTPCParam->GetNRow(isec);
1636
1637 row= new TObjArray* [nrows];
fe4da5cc 1638
73042f01 1639 MakeSector(isec,nrows,tH,ntracks,row);
8c555625 1640
1641 //--------------------------------------------------------
1642 // Digitize this sector, row by row
1643 // row[i] is the pointer to the TObjArray of TVectors,
1644 // each one containing electrons accepted on this
1645 // row, assigned into tracks
1646 //--------------------------------------------------------
1647
1648 Int_t i;
1649
2ab0c725 1650 if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree(fDigitsFile);
8c555625 1651
cc80f89e 1652 for (i=0;i<nrows;i++){
8c555625 1653
afc42102 1654 AliDigits * dig = fDigitsArray->CreateRow(isec,i);
8c555625 1655
cc80f89e 1656 DigitizeRow(i,isec,row);
8c555625 1657
cc80f89e 1658 fDigitsArray->StoreRow(isec,i);
8c555625 1659
9e1a0ddb 1660 Int_t ndig = dig->GetDigitSize();
5a28a08f 1661
cc80f89e 1662
a62d28a2 1663 //printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
cc80f89e 1664
1665 fDigitsArray->ClearRow(isec,i);
8c555625 1666
cc80f89e 1667
8c555625 1668 } // end of the sector digitization
8c555625 1669
cc80f89e 1670 for(i=0;i<nrows;i++){
39c8eb58 1671 row[i]->Delete();
1672 delete row[i];
cc80f89e 1673 }
1674
8c555625 1675 delete [] row; // delete the array of pointers to TObjArray-s
1676
1677 } // ntracks >0
8c555625 1678
cc80f89e 1679} // end of Hits2DigitsSector
8c555625 1680
8c555625 1681
8c555625 1682//_____________________________________________________________________________
cc80f89e 1683void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
8c555625 1684{
1685 //-----------------------------------------------------------
1686 // Single row digitization, coupling from the neighbouring
1687 // rows taken into account
1688 //-----------------------------------------------------------
1689
1690 //-----------------------------------------------------------------
1691 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
cc80f89e 1692 // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
8c555625 1693 //-----------------------------------------------------------------
1694
1695
8c555625 1696 Float_t zerosup = fTPCParam->GetZeroSup();
1697 Int_t nrows =fTPCParam->GetNRow(isec);
cc80f89e 1698 fCurrentIndex[1]= isec;
8c555625 1699
8c555625 1700
73042f01 1701 Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1702 Int_t nofTbins = fTPCParam->GetMaxTBin();
1703 Int_t indexRange[4];
8c555625 1704 //
1705 // Integrated signal for this row
1706 // and a single track signal
cc80f89e 1707 //
73042f01 1708 TMatrix *m1 = new TMatrix(0,nofPads,0,nofTbins); // integrated
1709 TMatrix *m2 = new TMatrix(0,nofPads,0,nofTbins); // single
8c555625 1710 //
73042f01 1711 TMatrix &total = *m1;
8c555625 1712
1713 // Array of pointers to the label-signal list
1714
73042f01 1715 Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1716 Float_t **pList = new Float_t* [nofDigits];
8c555625 1717
1718 Int_t lp;
cc80f89e 1719 Int_t i1;
73042f01 1720 for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
8c555625 1721 //
cc80f89e 1722 //calculate signal
1723 //
1724 Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1725 Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1726 for (Int_t row= row1;row<=row2;row++){
1727 Int_t nTracks= rows[row]->GetEntries();
1728 for (i1=0;i1<nTracks;i1++){
1729 fCurrentIndex[2]= row;
1730 fCurrentIndex[3]=irow;
1731 if (row==irow){
1732 m2->Zero(); // clear single track signal matrix
73042f01 1733 Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1734 GetList(trackLabel,nofPads,m2,indexRange,pList);
cc80f89e 1735 }
73042f01 1736 else GetSignal(rows[row],i1,0,m1,indexRange);
cc80f89e 1737 }
8c555625 1738 }
cc80f89e 1739
8c555625 1740 Int_t tracks[3];
8c555625 1741
cc80f89e 1742 AliDigits *dig = fDigitsArray->GetRow(isec,irow);
73042f01 1743 for(Int_t ip=0;ip<nofPads;ip++){
1744 for(Int_t it=0;it<nofTbins;it++){
8c555625 1745
73042f01 1746 Float_t q = total(ip,it);
8c555625 1747
73042f01 1748 Int_t gi =it*nofPads+ip; // global index
8c555625 1749
f8cf550c 1750 if(fDigitsSwitch == 0){
1751
1752 q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
1753
68771f83 1754 q = TMath::Nint(q);
8c555625 1755
f8cf550c 1756 if(q <=zerosup) continue; // do not fill zeros
1757 if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat(); // saturation
cc80f89e 1758
f8cf550c 1759 }
1760
1761 else {
5f16d237 1762 if(q <= 0.) continue; // do not fill zeros
68771f83 1763 if(q>2000.) q=2000.;
f8cf550c 1764 q *= 16.;
68771f83 1765 q = TMath::Nint(q);
f8cf550c 1766 }
8c555625 1767
1768 //
1769 // "real" signal or electronic noise (list = -1)?
1770 //
1771
1772 for(Int_t j1=0;j1<3;j1++){
1773 tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1774 }
1775
cc80f89e 1776//Begin_Html
1777/*
1778 <A NAME="AliDigits"></A>
1779 using of AliDigits object
1780*/
1781//End_Html
1782 dig->SetDigitFast((Short_t)q,it,ip);
1783 if (fDigitsArray->IsSimulated())
1784 {
1785 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1786 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1787 ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1788 }
1789
8c555625 1790
1791 } // end of loop over time buckets
1792 } // end of lop over pads
1793
1794 //
1795 // This row has been digitized, delete nonused stuff
1796 //
1797
73042f01 1798 for(lp=0;lp<nofDigits;lp++){
8c555625 1799 if(pList[lp]) delete [] pList[lp];
1800 }
1801
1802 delete [] pList;
1803
1804 delete m1;
1805 delete m2;
cc80f89e 1806 // delete m3;
8c555625 1807
1808} // end of DigitizeRow
cc80f89e 1809
8c555625 1810//_____________________________________________________________________________
cc80f89e 1811
1812Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
73042f01 1813 Int_t *indexRange)
8c555625 1814{
1815
1816 //---------------------------------------------------------------
1817 // Calculates 2-D signal (pad,time) for a single track,
1818 // returns a pointer to the signal matrix and the track label
1819 // No digitization is performed at this level!!!
1820 //---------------------------------------------------------------
1821
1822 //-----------------------------------------------------------------
1823 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
cc80f89e 1824 // Modified: Marian Ivanov
8c555625 1825 //-----------------------------------------------------------------
1826
1827 TVector *tv;
8c555625 1828
8c555625 1829 tv = (TVector*)p1->At(ntr); // pointer to a track
1830 TVector &v = *tv;
1831
1832 Float_t label = v(0);
73042f01 1833 Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
8c555625 1834
8c555625 1835 Int_t nElectrons = (tv->GetNrows()-1)/4;
73042f01 1836 indexRange[0]=9999; // min pad
1837 indexRange[1]=-1; // max pad
1838 indexRange[2]=9999; //min time
1839 indexRange[3]=-1; // max time
8c555625 1840
cc80f89e 1841 // Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1842
1843 TMatrix &signal = *m1;
1844 TMatrix &total = *m2;
8c555625 1845 //
1846 // Loop over all electrons
1847 //
8c555625 1848 for(Int_t nel=0; nel<nElectrons; nel++){
cc80f89e 1849 Int_t idx=nel*4;
1850 Float_t aval = v(idx+4);
1851 Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1852 Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1853 Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
8c555625 1854
cc80f89e 1855 if (n>0) for (Int_t i =0; i<n; i++){
1856 Int_t *index = fTPCParam->GetResBin(i);
73042f01 1857 Int_t pad=index[1]+centralPad; //in digit coordinates central pad has coordinate 0
cc80f89e 1858 if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
1859 Int_t time=index[2];
1860 Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
1861 weight *= eltoadcfac;
1862
1863 if (m1!=0) signal(pad,time)+=weight;
1864 total(pad,time)+=weight;
73042f01 1865 indexRange[0]=TMath::Min(indexRange[0],pad);
1866 indexRange[1]=TMath::Max(indexRange[1],pad);
1867 indexRange[2]=TMath::Min(indexRange[2],time);
1868 indexRange[3]=TMath::Max(indexRange[3],time);
cc80f89e 1869 }
1870 }
8c555625 1871 } // end of loop over electrons
cc80f89e 1872
8c555625 1873 return label; // returns track label when finished
1874}
1875
1876//_____________________________________________________________________________
73042f01 1877void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
8c555625 1878 Float_t **pList)
1879{
1880 //----------------------------------------------------------------------
1881 // Updates the list of tracks contributing to digits for a given row
1882 //----------------------------------------------------------------------
1883
1884 //-----------------------------------------------------------------
1885 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1886 //-----------------------------------------------------------------
1887
1888 TMatrix &signal = *m;
1889
1890 // lop over nonzero digits
1891
73042f01 1892 for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1893 for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
8c555625 1894
1895
921bf71a 1896 // accept only the contribution larger than 500 electrons (1/2 s_noise)
1897
cc80f89e 1898 if(signal(ip,it)<0.5) continue;
921bf71a 1899
1900
73042f01 1901 Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
8c555625 1902
73042f01 1903 if(!pList[globalIndex]){
8c555625 1904
1905 //
1906 // Create new list (6 elements - 3 signals and 3 labels),
8c555625 1907 //
1908
73042f01 1909 pList[globalIndex] = new Float_t [6];
8c555625 1910
1911 // set list to -1
1912
73042f01 1913 *pList[globalIndex] = -1.;
1914 *(pList[globalIndex]+1) = -1.;
1915 *(pList[globalIndex]+2) = -1.;
1916 *(pList[globalIndex]+3) = -1.;
1917 *(pList[globalIndex]+4) = -1.;
1918 *(pList[globalIndex]+5) = -1.;
8c555625 1919
1920
73042f01 1921 *pList[globalIndex] = label;
1922 *(pList[globalIndex]+3) = signal(ip,it);
8c555625 1923 }
1924 else{
1925
1926 // check the signal magnitude
1927
73042f01 1928 Float_t highest = *(pList[globalIndex]+3);
1929 Float_t middle = *(pList[globalIndex]+4);
1930 Float_t lowest = *(pList[globalIndex]+5);
8c555625 1931
1932 //
1933 // compare the new signal with already existing list
1934 //
1935
1936 if(signal(ip,it)<lowest) continue; // neglect this track
1937
1938 //
1939
1940 if (signal(ip,it)>highest){
73042f01 1941 *(pList[globalIndex]+5) = middle;
1942 *(pList[globalIndex]+4) = highest;
1943 *(pList[globalIndex]+3) = signal(ip,it);
8c555625 1944
73042f01 1945 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1946 *(pList[globalIndex]+1) = *pList[globalIndex];
1947 *pList[globalIndex] = label;
8c555625 1948 }
1949 else if (signal(ip,it)>middle){
73042f01 1950 *(pList[globalIndex]+5) = middle;
1951 *(pList[globalIndex]+4) = signal(ip,it);
8c555625 1952
73042f01 1953 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1954 *(pList[globalIndex]+1) = label;
8c555625 1955 }
1956 else{
73042f01 1957 *(pList[globalIndex]+5) = signal(ip,it);
1958 *(pList[globalIndex]+2) = label;
8c555625 1959 }
1960 }
1961
1962 } // end of loop over pads
1963 } // end of loop over time bins
1964
1965
1966
8c555625 1967}//end of GetList
1968//___________________________________________________________________
1969void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1970 Stat_t ntracks,TObjArray **row)
1971{
1972
1973 //-----------------------------------------------------------------
1974 // Prepares the sector digitization, creates the vectors of
1975 // tracks for each row of this sector. The track vector
1976 // contains the track label and the position of electrons.
1977 //-----------------------------------------------------------------
1978
1979 //-----------------------------------------------------------------
1980 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1981 //-----------------------------------------------------------------
1982
cc80f89e 1983 Float_t gasgain = fTPCParam->GetGasGain();
8c555625 1984 Int_t i;
cc80f89e 1985 Float_t xyz[4];
8c555625 1986
1987 AliTPChit *tpcHit; // pointer to a sigle TPC hit
39c8eb58 1988 //MI change
1989 TBranch * branch=0;
1990 if (fHitType&2) branch = TH->GetBranch("TPC2");
1991 else branch = TH->GetBranch("TPC");
1992
8c555625 1993
1994 //----------------------------------------------
1995 // Create TObjArray-s, one for each row,
1996 // each TObjArray will store the TVectors
1997 // of electrons, one TVector per each track.
1998 //----------------------------------------------
1999
f74bb6f5 2000 Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
2001 TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
8c555625 2002 for(i=0; i<nrows; i++){
2003 row[i] = new TObjArray;
f74bb6f5 2004 nofElectrons[i]=0;
2005 tracks[i]=0;
8c555625 2006 }
8c555625 2007
37831078 2008
2009
8c555625 2010 //--------------------------------------------------------------------
2011 // Loop over tracks, the "track" contains the full history
2012 //--------------------------------------------------------------------
2013
2014 Int_t previousTrack,currentTrack;
2015 previousTrack = -1; // nothing to store so far!
2016
2017 for(Int_t track=0;track<ntracks;track++){
39c8eb58 2018 Bool_t isInSector=kTRUE;
8c555625 2019 ResetHits();
2020
39c8eb58 2021 if (fHitType&2) {
2022 isInSector=kFALSE;
2023 TBranch * br = TH->GetBranch("fTrackHitsInfo");
2024 br->GetEvent(track);
2025 AliObjectArray * ar = fTrackHits->fTrackHitsInfo;
2026 for (UInt_t j=0;j<ar->GetSize();j++){
2027 if ( ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE;
2028 }
2029 }
2030 if (!isInSector) continue;
2031 //MI change
2032 branch->GetEntry(track); // get next track
2033
2034 //M.I. changes
8c555625 2035
39c8eb58 2036 tpcHit = (AliTPChit*)FirstHit(-1);
8c555625 2037
2038 //--------------------------------------------------------------
2039 // Loop over hits
2040 //--------------------------------------------------------------
2041
8c555625 2042
39c8eb58 2043 while(tpcHit){
8c555625 2044
2045 Int_t sector=tpcHit->fSector; // sector number
39c8eb58 2046 // if(sector != isec) continue;
2047 if(sector != isec){
2048 tpcHit = (AliTPChit*) NextHit();
2049 continue;
2050 }
8c555625 2051
94de3818 2052 currentTrack = tpcHit->Track(); // track number
39c8eb58 2053
2054
8c555625 2055 if(currentTrack != previousTrack){
2056
2057 // store already filled fTrack
2058
2059 for(i=0;i<nrows;i++){
2060 if(previousTrack != -1){
73042f01 2061 if(nofElectrons[i]>0){
cc80f89e 2062 TVector &v = *tracks[i];
8c555625 2063 v(0) = previousTrack;
73042f01 2064 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
cc80f89e 2065 row[i]->Add(tracks[i]);
8c555625 2066 }
2067 else{
cc80f89e 2068 delete tracks[i]; // delete empty TVector
2069 tracks[i]=0;
8c555625 2070 }
2071 }
2072
73042f01 2073 nofElectrons[i]=0;
cc80f89e 2074 tracks[i] = new TVector(481); // TVectors for the next fTrack
8c555625 2075
2076 } // end of loop over rows
2077
2078 previousTrack=currentTrack; // update track label
2079 }
2080
73042f01 2081 Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
8c555625 2082
2083 //---------------------------------------------------
2084 // Calculate the electron attachment probability
2085 //---------------------------------------------------
2086
cc80f89e 2087
94de3818 2088 Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
cc80f89e 2089 /fTPCParam->GetDriftV();
8c555625 2090 // in microseconds!
73042f01 2091 Float_t attProb = fTPCParam->GetAttCoef()*
8c555625 2092 fTPCParam->GetOxyCont()*time; // fraction!
2093
2094 //-----------------------------------------------
2095 // Loop over electrons
2096 //-----------------------------------------------
cc80f89e 2097 Int_t index[3];
2098 index[1]=isec;
73042f01 2099 for(Int_t nel=0;nel<qI;nel++){
8c555625 2100 // skip if electron lost due to the attachment
73042f01 2101 if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
94de3818 2102 xyz[0]=tpcHit->X();
2103 xyz[1]=tpcHit->Y();
1705494b 2104 xyz[2]=tpcHit->Z();
2105 //
2106 // protection for the nonphysical avalanche size (10**6 maximum)
2107 //
710dacf5 2108 Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2109 xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
cc80f89e 2110 index[0]=1;
2111
2112 TransportElectron(xyz,index); //MI change -august
73042f01 2113 Int_t rowNumber;
cc80f89e 2114 fTPCParam->GetPadRow(xyz,index); //MI change august
73042f01 2115 rowNumber = index[2];
cc80f89e 2116 //transform position to local digit coordinates
2117 //relative to nearest pad row
73042f01 2118 if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue;
2119 nofElectrons[rowNumber]++;
8c555625 2120 //----------------------------------
2121 // Expand vector if necessary
2122 //----------------------------------
73042f01 2123 if(nofElectrons[rowNumber]>120){
2124 Int_t range = tracks[rowNumber]->GetNrows();
2125 if((nofElectrons[rowNumber])>(range-1)/4){
cc80f89e 2126
73042f01 2127 tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
fe4da5cc 2128 }
2129 }
2130
73042f01 2131 TVector &v = *tracks[rowNumber];
2132 Int_t idx = 4*nofElectrons[rowNumber]-3;
8c555625 2133
cc80f89e 2134 v(idx)= xyz[0]; // X - pad row coordinate
2135 v(idx+1)=xyz[1]; // Y - pad coordinate (along the pad-row)
2136 v(idx+2)=xyz[2]; // Z - time bin coordinate
2137 v(idx+3)=xyz[3]; // avalanche size
8c555625 2138 } // end of loop over electrons
39c8eb58 2139
2140 tpcHit = (AliTPChit*)NextHit();
8c555625 2141
2142 } // end of loop over hits
2143 } // end of loop over tracks
2144
2145 //
2146 // store remaining track (the last one) if not empty
2147 //
2148
2149 for(i=0;i<nrows;i++){
73042f01 2150 if(nofElectrons[i]>0){
cc80f89e 2151 TVector &v = *tracks[i];
8c555625 2152 v(0) = previousTrack;
73042f01 2153 tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
cc80f89e 2154 row[i]->Add(tracks[i]);
fe4da5cc 2155 }
2156 else{
cc80f89e 2157 delete tracks[i];
2158 tracks[i]=0;
8c555625 2159 }
2160 }
2161
cc80f89e 2162 delete [] tracks;
73042f01 2163 delete [] nofElectrons;
8c555625 2164
8c555625 2165
cc80f89e 2166} // end of MakeSector
8c555625 2167
fe4da5cc 2168
2169//_____________________________________________________________________________
2170void AliTPC::Init()
2171{
2172 //
2173 // Initialise TPC detector after definition of geometry
2174 //
2175 Int_t i;
2176 //
9e1a0ddb 2177 if(fDebug) {
2178 printf("\n%s: ",ClassName());
2179 for(i=0;i<35;i++) printf("*");
2180 printf(" TPC_INIT ");
2181 for(i=0;i<35;i++) printf("*");
2182 printf("\n%s: ",ClassName());
2183 //
2184 for(i=0;i<80;i++) printf("*");
2185 printf("\n");
2186 }
fe4da5cc 2187}
2188
2189//_____________________________________________________________________________
9e1a0ddb 2190void AliTPC::MakeBranch(Option_t* option, const char *file)
fe4da5cc 2191{
2192 //
2193 // Create Tree branches for the TPC.
2194 //
2195 Int_t buffersize = 4000;
2196 char branchname[10];
2197 sprintf(branchname,"%s",GetName());
2198
2ab0c725 2199 AliDetector::MakeBranch(option,file);
fe4da5cc 2200
5cf7bbad 2201 const char *d = strstr(option,"D");
fe4da5cc 2202
73042f01 2203 if (fDigits && gAlice->TreeD() && d) {
9e1a0ddb 2204 MakeBranchInTree(gAlice->TreeD(),
2205 branchname, &fDigits, buffersize, file);
fe4da5cc 2206 }
39c8eb58 2207
2ab0c725 2208 if (fHitType&2) MakeBranch2(option,file); // MI change 14.09.2000
fe4da5cc 2209}
2210
2211//_____________________________________________________________________________
2212void AliTPC::ResetDigits()
2213{
2214 //
2215 // Reset number of digits and the digits array for this detector
fe4da5cc 2216 //
2217 fNdigits = 0;
cc80f89e 2218 if (fDigits) fDigits->Clear();
fe4da5cc 2219}
2220
2221//_____________________________________________________________________________
2222void AliTPC::SetSecAL(Int_t sec)
2223{
8c555625 2224 //---------------------------------------------------
fe4da5cc 2225 // Activate/deactivate selection for lower sectors
8c555625 2226 //---------------------------------------------------
2227
2228 //-----------------------------------------------------------------
2229 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2230 //-----------------------------------------------------------------
2231
fe4da5cc 2232 fSecAL = sec;
2233}
2234
2235//_____________________________________________________________________________
2236void AliTPC::SetSecAU(Int_t sec)
2237{
8c555625 2238 //----------------------------------------------------
fe4da5cc 2239 // Activate/deactivate selection for upper sectors
8c555625 2240 //---------------------------------------------------
2241
2242 //-----------------------------------------------------------------
2243 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2244 //-----------------------------------------------------------------
2245
fe4da5cc 2246 fSecAU = sec;
2247}
2248
2249//_____________________________________________________________________________
2250void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2251{
8c555625 2252 //----------------------------------------
fe4da5cc 2253 // Select active lower sectors
8c555625 2254 //----------------------------------------
2255
2256 //-----------------------------------------------------------------
2257 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2258 //-----------------------------------------------------------------
2259
fe4da5cc 2260 fSecLows[0] = s1;
2261 fSecLows[1] = s2;
2262 fSecLows[2] = s3;
2263 fSecLows[3] = s4;
2264 fSecLows[4] = s5;
2265 fSecLows[5] = s6;
2266}
2267
2268//_____________________________________________________________________________
2269void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2270 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2271 Int_t s11 , Int_t s12)
2272{
8c555625 2273 //--------------------------------
fe4da5cc 2274 // Select active upper sectors
8c555625 2275 //--------------------------------
2276
2277 //-----------------------------------------------------------------
2278 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2279 //-----------------------------------------------------------------
2280
fe4da5cc 2281 fSecUps[0] = s1;
2282 fSecUps[1] = s2;
2283 fSecUps[2] = s3;
2284 fSecUps[3] = s4;
2285 fSecUps[4] = s5;
2286 fSecUps[5] = s6;
2287 fSecUps[6] = s7;
2288 fSecUps[7] = s8;
2289 fSecUps[8] = s9;
2290 fSecUps[9] = s10;
2291 fSecUps[10] = s11;
2292 fSecUps[11] = s12;
2293}
2294
2295//_____________________________________________________________________________
2296void AliTPC::SetSens(Int_t sens)
2297{
8c555625 2298
2299 //-------------------------------------------------------------
2300 // Activates/deactivates the sensitive strips at the center of
2301 // the pad row -- this is for the space-point resolution calculations
2302 //-------------------------------------------------------------
2303
2304 //-----------------------------------------------------------------
2305 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2306 //-----------------------------------------------------------------
2307
fe4da5cc 2308 fSens = sens;
2309}
2b06d5c3 2310
4b0fdcad 2311
73042f01 2312void AliTPC::SetSide(Float_t side=0.)
4b0fdcad 2313{
73042f01 2314 // choice of the TPC side
2315
4b0fdcad 2316 fSide = side;
2317
2318}
1283eee5 2319//____________________________________________________________________________
2320void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2321 Float_t p2,Float_t p3)
2322{
fe4da5cc 2323
73042f01 2324 // gax mixture definition
2325
1283eee5 2326 fNoComp = nc;
2327
2328 fMixtComp[0]=c1;
2329 fMixtComp[1]=c2;
2330 fMixtComp[2]=c3;
2331
2332 fMixtProp[0]=p1;
2333 fMixtProp[1]=p2;
2334 fMixtProp[2]=p3;
2335
2336
cc80f89e 2337}
2338//_____________________________________________________________________________
2339
2340void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2341{
2342 //
2343 // electron transport taking into account:
2344 // 1. diffusion,
2345 // 2.ExB at the wires
2346 // 3. nonisochronity
2347 //
2348 // xyz and index must be already transformed to system 1
2349 //
2350
2351 fTPCParam->Transform1to2(xyz,index);
2352
2353 //add diffusion
2354 Float_t driftl=xyz[2];
2355 if(driftl<0.01) driftl=0.01;
2356 driftl=TMath::Sqrt(driftl);
73042f01 2357 Float_t sigT = driftl*(fTPCParam->GetDiffT());
2358 Float_t sigL = driftl*(fTPCParam->GetDiffL());
2359 xyz[0]=gRandom->Gaus(xyz[0],sigT);
2360 xyz[1]=gRandom->Gaus(xyz[1],sigT);
2361 xyz[2]=gRandom->Gaus(xyz[2],sigL);
cc80f89e 2362
2363 // ExB
2364
2365 if (fTPCParam->GetMWPCReadout()==kTRUE){
2366 Float_t x1=xyz[0];
2367 fTPCParam->Transform2to2NearestWire(xyz,index);
2368 Float_t dx=xyz[0]-x1;
2369 xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2370 }
2371 //add nonisochronity (not implemented yet)
2372
1283eee5 2373}
fe4da5cc 2374
fe4da5cc 2375ClassImp(AliTPCdigit)
2376
2377//_____________________________________________________________________________
2378AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2379 AliDigit(tracks)
2380{
2381 //
2382 // Creates a TPC digit object
2383 //
2384 fSector = digits[0];
2385 fPadRow = digits[1];
2386 fPad = digits[2];
2387 fTime = digits[3];
2388 fSignal = digits[4];
2389}
2390
2391
2392ClassImp(AliTPChit)
2393
2394//_____________________________________________________________________________
2395AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2396AliHit(shunt,track)
2397{
2398 //
2399 // Creates a TPC hit object
2400 //
2401 fSector = vol[0];
2402 fPadRow = vol[1];
2403 fX = hits[0];
2404 fY = hits[1];
2405 fZ = hits[2];
2406 fQ = hits[3];
2407}
2408
8c555625 2409
39c8eb58 2410//________________________________________________________________________
2411// Additional code because of the AliTPCTrackHits
2412
9e1a0ddb 2413void AliTPC::MakeBranch2(Option_t *option,const char *file)
39c8eb58 2414{
2415 //
2416 // Create a new branch in the current Root Tree
2417 // The branch of fHits is automatically split
2418 // MI change 14.09.2000
2419 if (fHitType&2==0) return;
2420 char branchname[10];
2421 sprintf(branchname,"%s2",GetName());
2422 //
2423 // Get the pointer to the header
5cf7bbad 2424 const char *cH = strstr(option,"H");
39c8eb58 2425 //
2426 if (fTrackHits && gAlice->TreeH() && cH) {
39c8eb58 2427 AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHits,
8f90302f 2428 gAlice->TreeH(),fBufferSize,99);
39c8eb58 2429 gAlice->TreeH()->GetListOfBranches()->Add(branch);
9e1a0ddb 2430 if (GetDebug()>1)
2431 printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname);
2432 const char folder [] = "RunMC/Event/Data";
2433 if (GetDebug())
2434 printf("%15s: Publishing %s to %s\n",ClassName(),branchname,folder);
2435 Publish(folder,&fTrackHits,branchname);
2ab0c725 2436 if (file) {
2437 TBranch *b = gAlice->TreeH()->GetBranch(branchname);
2438 TDirectory *wd = gDirectory;
2439 b->SetFile(file);
2440 TIter next( b->GetListOfBranches());
2441 while ((b=(TBranch*)next())) {
9e1a0ddb 2442 b->SetFile(file);
2ab0c725 2443 }
2444 wd->cd();
9e1a0ddb 2445 if (GetDebug()>1)
2446 cout << "Diverting branch " << branchname << " to file " << file << endl;
2ab0c725 2447 }
39c8eb58 2448 }
2449}
2450
2451void AliTPC::SetTreeAddress()
2452{
2453 if (fHitType&1) AliDetector::SetTreeAddress();
2454 if (fHitType&2) SetTreeAddress2();
2455}
2456
2457void AliTPC::SetTreeAddress2()
2458{
2459 //
2460 // Set branch address for the TrackHits Tree
2461 //
2462 TBranch *branch;
2463 char branchname[20];
2464 sprintf(branchname,"%s2",GetName());
2465 //
2466 // Branch address for hit tree
2467 TTree *treeH = gAlice->TreeH();
2468 if (treeH) {
2469 branch = treeH->GetBranch(branchname);
2470 if (branch) branch->SetAddress(&fTrackHits);
2471 }
2472}
2473
2474void AliTPC::FinishPrimary()
2475{
2476 if (fTrackHits) fTrackHits->FlushHitStack();
2477}
2478
2479
2480void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2481{
2482 //
2483 // add hit to the list
39c8eb58 2484 Int_t rtrack;
2485 if (fIshunt) {
2486 int primary = gAlice->GetPrimary(track);
2ab0c725 2487 gAlice->Particle(primary)->SetBit(kKeepBit);
39c8eb58 2488 rtrack=primary;
2489 } else {
2490 rtrack=track;
2491 gAlice->FlagTrack(track);
2492 }
2493 //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
2494 //if (hit->fTrack!=rtrack)
2495 // cout<<"bad track number\n";
2496 if (fTrackHits)
2497 fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2498 hits[1],hits[2],(Int_t)hits[3]);
2499}
2500
2501void AliTPC::ResetHits()
2502{
2503 if (fHitType&1) AliDetector::ResetHits();
2504 if (fHitType&2) ResetHits2();
2505}
2506
2507void AliTPC::ResetHits2()
2508{
2509 //
2510 //reset hits
2511 if (fTrackHits) fTrackHits->Clear();
2512}
2513
2514AliHit* AliTPC::FirstHit(Int_t track)
2515{
2516 if (fHitType&2) return FirstHit2(track);
2517 return AliDetector::FirstHit(track);
2518}
2519AliHit* AliTPC::NextHit()
2520{
2521 if (fHitType&2) return NextHit2();
2522 return AliDetector::NextHit();
2523}
2524
2525AliHit* AliTPC::FirstHit2(Int_t track)
2526{
2527 //
2528 // Initialise the hit iterator
2529 // Return the address of the first hit for track
2530 // If track>=0 the track is read from disk
2531 // while if track<0 the first hit of the current
2532 // track is returned
2533 //
2534 if(track>=0) {
2535 gAlice->ResetHits();
2536 gAlice->TreeH()->GetEvent(track);
2537 }
2538 //
2539 if (fTrackHits) {
2540 fTrackHits->First();
2541 return fTrackHits->GetHit();
2542 }
2543 else return 0;
2544}
2545
2546AliHit* AliTPC::NextHit2()
2547{
2548 //
2549 //Return the next hit for the current track
2550
2551 if (fTrackHits) {
2552 fTrackHits->Next();
2553 return fTrackHits->GetHit();
2554 }
2555 else
2556 return 0;
2557}
2558
2559void AliTPC::LoadPoints(Int_t)
2560{
2561 //
2562 Int_t a = 0;
f2e8b846 2563 /* if(fHitType==1) return AliDetector::LoadPoints(a);
39c8eb58 2564 LoadPoints2(a);
f2e8b846 2565 */
2566 if(fHitType==1) AliDetector::LoadPoints(a);
2567 else LoadPoints2(a);
2568
39c8eb58 2569 // LoadPoints3(a);
2570
2571}
2572
2573
2574void AliTPC::RemapTrackHitIDs(Int_t *map)
2575{
2576 if (!fTrackHits) return;
2577 AliObjectArray * arr = fTrackHits->fTrackHitsInfo;
2578 for (UInt_t i=0;i<arr->GetSize();i++){
2579 AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2580 info->fTrackID = map[info->fTrackID];
2581 }
2582
2583}
2584
2585
2586//_____________________________________________________________________________
2587void AliTPC::LoadPoints2(Int_t)
2588{
2589 //
2590 // Store x, y, z of all hits in memory
2591 //
2592 if (fTrackHits == 0) return;
2593 //
2594 Int_t nhits = fTrackHits->GetEntriesFast();
2595 if (nhits == 0) return;
2596 Int_t tracks = gAlice->GetNtrack();
2597 if (fPoints == 0) fPoints = new TObjArray(tracks);
2598 AliHit *ahit;
2599 //
2600 Int_t *ntrk=new Int_t[tracks];
2601 Int_t *limi=new Int_t[tracks];
2602 Float_t **coor=new Float_t*[tracks];
2603 for(Int_t i=0;i<tracks;i++) {
2604 ntrk[i]=0;
2605 coor[i]=0;
2606 limi[i]=0;
2607 }
2608 //
2609 AliPoints *points = 0;
2610 Float_t *fp=0;
2611 Int_t trk;
2612 Int_t chunk=nhits/4+1;
2613 //
2614 // Loop over all the hits and store their position
2615 //
2616 ahit = FirstHit2(-1);
2617 //for (Int_t hit=0;hit<nhits;hit++) {
2618 while (ahit){
2619 // ahit = (AliHit*)fHits->UncheckedAt(hit);
2620 trk=ahit->GetTrack();
2621 if(ntrk[trk]==limi[trk]) {
2622 //
2623 // Initialise a new track
2624 fp=new Float_t[3*(limi[trk]+chunk)];
2625 if(coor[trk]) {
2626 memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2627 delete [] coor[trk];
2628 }
2629 limi[trk]+=chunk;
2630 coor[trk] = fp;
2631 } else {
2632 fp = coor[trk];
2633 }
2634 fp[3*ntrk[trk] ] = ahit->X();
2635 fp[3*ntrk[trk]+1] = ahit->Y();
2636 fp[3*ntrk[trk]+2] = ahit->Z();
2637 ntrk[trk]++;
2638 ahit = NextHit2();
2639 }
2640 //
2641 for(trk=0; trk<tracks; ++trk) {
2642 if(ntrk[trk]) {
2643 points = new AliPoints();
2644 points->SetMarkerColor(GetMarkerColor());
2645 points->SetMarkerSize(GetMarkerSize());
2646 points->SetDetector(this);
2647 points->SetParticle(trk);
2648 points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
2649 fPoints->AddAt(points,trk);
2650 delete [] coor[trk];
2651 coor[trk]=0;
2652 }
2653 }
2654 delete [] coor;
2655 delete [] ntrk;
2656 delete [] limi;
2657}
2658
2659
2660//_____________________________________________________________________________
2661void AliTPC::LoadPoints3(Int_t)
2662{
2663 //
2664 // Store x, y, z of all hits in memory
2665 // - only intersection point with pad row
2666 if (fTrackHits == 0) return;
2667 //
2668 Int_t nhits = fTrackHits->GetEntriesFast();
2669 if (nhits == 0) return;
2670 Int_t tracks = gAlice->GetNtrack();
2671 if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2672 fPoints->Expand(2*tracks);
2673 AliHit *ahit;
2674 //
2675 Int_t *ntrk=new Int_t[tracks];
2676 Int_t *limi=new Int_t[tracks];
2677 Float_t **coor=new Float_t*[tracks];
2678 for(Int_t i=0;i<tracks;i++) {
2679 ntrk[i]=0;
2680 coor[i]=0;
2681 limi[i]=0;
2682 }
2683 //
2684 AliPoints *points = 0;
2685 Float_t *fp=0;
2686 Int_t trk;
2687 Int_t chunk=nhits/4+1;
2688 //
2689 // Loop over all the hits and store their position
2690 //
2691 ahit = FirstHit2(-1);
2692 //for (Int_t hit=0;hit<nhits;hit++) {
2693
2694 Int_t lastrow = -1;
2695 while (ahit){
2696 // ahit = (AliHit*)fHits->UncheckedAt(hit);
2697 trk=ahit->GetTrack();
2698 Float_t x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2699 Int_t index[3]={1,((AliTPChit*)ahit)->fSector,0};
2700 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
2701 if (currentrow!=lastrow){
2702 lastrow = currentrow;
2703 //later calculate intersection point
2704 if(ntrk[trk]==limi[trk]) {
2705 //
2706 // Initialise a new track
2707 fp=new Float_t[3*(limi[trk]+chunk)];
2708 if(coor[trk]) {
2709 memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2710 delete [] coor[trk];
2711 }
2712 limi[trk]+=chunk;
2713 coor[trk] = fp;
2714 } else {
2715 fp = coor[trk];
2716 }
2717 fp[3*ntrk[trk] ] = ahit->X();
2718 fp[3*ntrk[trk]+1] = ahit->Y();
2719 fp[3*ntrk[trk]+2] = ahit->Z();
2720 ntrk[trk]++;
2721 }
2722 ahit = NextHit2();
2723 }
2724
2725 //
2726 for(trk=0; trk<tracks; ++trk) {
2727 if(ntrk[trk]) {
2728 points = new AliPoints();
2729 points->SetMarkerColor(GetMarkerColor()+1);
2730 points->SetMarkerStyle(5);
2731 points->SetMarkerSize(0.2);
2732 points->SetDetector(this);
2733 points->SetParticle(trk);
2734 // points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20);
2735 points->SetPolyMarker(ntrk[trk],coor[trk],30);
2736 fPoints->AddAt(points,tracks+trk);
2737 delete [] coor[trk];
2738 coor[trk]=0;
2739 }
2740 }
2741 delete [] coor;
2742 delete [] ntrk;
2743 delete [] limi;
2744}
2745
2746
2747
2748void AliTPC::FindTrackHitsIntersection(TClonesArray * arr)
2749{
2750
2751 //
2752 //fill clones array with intersection of current point with the
2753 //middle of the row
2754 Int_t sector;
2755 Int_t ipart;
2756
2757 const Int_t kcmaxhits=30000;
2758 TVector * xxxx = new TVector(kcmaxhits*4);
2759 TVector & xxx = *xxxx;
2760 Int_t maxhits = kcmaxhits;
2761
2762 //
2763 AliTPChit * tpcHit=0;
2764 tpcHit = (AliTPChit*)FirstHit2(-1);
2765 Int_t currentIndex=0;
2766 Int_t lastrow=-1; //last writen row
2767
2768 while (tpcHit){
2769 if (tpcHit==0) continue;
2770 sector=tpcHit->fSector; // sector number
2771 ipart=tpcHit->Track();
2772
2773 //find row number
2774
2775 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
2776 Int_t index[3]={1,sector,0};
2777 Int_t currentrow = fTPCParam->GetPadRow(x,index) ;
2778 if (currentrow<0) continue;
2779 if (lastrow<0) lastrow=currentrow;
2780 if (currentrow==lastrow){
2781 if ( currentIndex>=maxhits){
2782 maxhits+=kcmaxhits;
2783 xxx.ResizeTo(4*maxhits);
2784 }
2785 xxx(currentIndex*4)=x[0];
2786 xxx(currentIndex*4+1)=x[1];
2787 xxx(currentIndex*4+2)=x[2];
2788 xxx(currentIndex*4+3)=tpcHit->fQ;
2789 currentIndex++;
2790 }
2791 else
2792 if (currentIndex>2){
2793 Float_t sumx=0;
2794 Float_t sumx2=0;
2795 Float_t sumx3=0;
2796 Float_t sumx4=0;
2797 Float_t sumy=0;
2798 Float_t sumxy=0;
2799 Float_t sumx2y=0;
2800 Float_t sumz=0;
2801 Float_t sumxz=0;
2802 Float_t sumx2z=0;
2803 Float_t sumq=0;
2804 for (Int_t index=0;index<currentIndex;index++){
2805 Float_t x,x2,x3,x4;
2806 x=x2=x3=x4=xxx(index*4);
2807 x2*=x;
2808 x3*=x2;
2809 x4*=x3;
2810 sumx+=x;
2811 sumx2+=x2;
2812 sumx3+=x3;
2813 sumx4+=x4;
2814 sumy+=xxx(index*4+1);
2815 sumxy+=xxx(index*4+1)*x;
2816 sumx2y+=xxx(index*4+1)*x2;
2817 sumz+=xxx(index*4+2);
2818 sumxz+=xxx(index*4+2)*x;
2819 sumx2z+=xxx(index*4+2)*x2;
2820 sumq+=xxx(index*4+3);
2821 }
2822 Float_t centralPad = (fTPCParam->GetNPads(sector,lastrow)-1)/2;
2823 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
2824 sumx2*(sumx*sumx3-sumx2*sumx2);
2825
2826 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
2827 sumx2*(sumxy*sumx3-sumx2y*sumx2);
2828 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
2829 sumx2*(sumxz*sumx3-sumx2z*sumx2);
2830
2831 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
2832 sumx2*(sumx*sumx2y-sumx2*sumxy);
2833 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
2834 sumx2*(sumx*sumx2z-sumx2*sumxz);
2835
2836 Float_t y=detay/det+centralPad;
2837 Float_t z=detaz/det;
2838 Float_t by=detby/det; //y angle
2839 Float_t bz=detbz/det; //z angle
2840 sumy/=Float_t(currentIndex);
2841 sumz/=Float_t(currentIndex);
2842
2843 AliComplexCluster cl;
2844 cl.fX=z;
2845 cl.fY=y;
2846 cl.fQ=sumq;
2847 cl.fSigmaX2=bz;
2848 cl.fSigmaY2=by;
2849 cl.fTracks[0]=ipart;
2850
2851 AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow));
2852 if (row!=0) row->InsertCluster(&cl);
2853 currentIndex=0;
2854 lastrow=currentrow;
2855 } //end of calculating cluster for given row
2856
2857 } // end of loop over hits
2858 xxxx->Delete();
39c8eb58 2859
2860}
afc42102 2861//_______________________________________________________________________________
a62d28a2 2862void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
afc42102 2863{
a62d28a2 2864 // produces rec points from digits and writes them on the root file
2865 // AliTPCclusters.root
2866
2867 TDirectory *cwd = gDirectory;
2868
2869
2870 AliTPCParam *dig=(AliTPCParam *)gDirectory->Get("75x40_100x60");
2871 SetParam(dig);
2872 cout<<"AliTPC::Digits2Reco: TPC parameteres have been set"<<endl;
2873 TFile *out;
2874 if(!gSystem->Getenv("CONFIG_FILE")){
2875 out=TFile::Open("AliTPCclusters.root","recreate");
2876 }
2877 else {
2878 const char *tmp1;
2879 const char *tmp2;
2880 char tmp3[80];
2881 tmp1=gSystem->Getenv("CONFIG_FILE_PREFIX");
2882 tmp2=gSystem->Getenv("CONFIG_OUTDIR");
2883 sprintf(tmp3,"%s%s/AliTPCclusters.root",tmp1,tmp2);
2884 out=TFile::Open(tmp3,"recreate");
2885 }
2886
2887 TStopwatch timer;
2888 cout<<"AliTPC::Digits2Reco - determination of rec points begins"<<endl;
2889 timer.Start();
2890 cwd->cd();
2891 for(Int_t iev=firstevent;iev<lastevent+1;iev++){
2892
2893 printf("Processing event %d\n",iev);
2894 Digits2Clusters(out,iev);
2895 }
2896 cout<<"AliTPC::Digits2Reco - determination of rec points ended"<<endl;
2897 timer.Stop();
2898 timer.Print();
2899 out->Close();
2900 cwd->cd();
2901
2902
afc42102 2903}