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