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