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