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