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