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