15 minutes interval for calculation of drift correction
[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 "AliTPCRawStreamV3.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   //
840   AliSimDigits digarr;
841   AliSimDigits* digrow = &digarr;
842   digits->GetBranch("Segment")->SetAddress(&digrow);
843
844   const char* fileName = "AliTPCDDL.dat";
845   AliTPCBuffer* buffer  = new AliTPCBuffer(fileName);
846   //Verbose level
847   // 0: Silent
848   // 1: cout messages
849   // 2: txt files with digits 
850   //BE CAREFUL, verbose level 2 MUST be used only for debugging and
851   //it is highly suggested to use this mode only for debugging digits files
852   //reasonably small, because otherwise the size of the txt files can reach
853   //quickly several MB wasting time and disk space.
854   buffer->SetVerbose(0);
855
856   Int_t nEntries = Int_t(digits->GetEntries());
857   Int_t previousSector = -1;
858   Int_t subSector = 0;
859   for (Int_t i = 0; i < nEntries; i++) {
860     digits->GetEntry(i);
861     Int_t sector, row;
862     fTPCParam->AdjustSectorRow(digarr.GetID(), sector, row);
863     if(previousSector != sector) {
864       subSector = 0;
865       previousSector = sector;
866     }
867
868     if (sector < 36) { //inner sector [0;35]
869       if (row != 30) {
870         //the whole row is written into the output file
871         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0, 
872                                sector, subSector, row);
873       } else {
874         //only the pads in the range [37;48] are written into the output file
875         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 1, 
876                                sector, subSector, row);
877         subSector = 1;
878         //only the pads outside the range [37;48] are written into the output file
879         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 2, 
880                                sector, subSector, row);
881       }//end else
882
883     } else { //outer sector [36;71]
884       if (row == 54) subSector = 2;
885       if ((row != 27) && (row != 76)) {
886         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
887                                sector, subSector, row);
888       } else if (row == 27) {
889         //only the pads outside the range [43;46] are written into the output file
890         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 2,
891                                  sector, subSector, row);
892         subSector = 1;
893         //only the pads in the range [43;46] are written into the output file
894         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 1,
895                                  sector, subSector, row);
896       } else if (row == 76) {
897         //only the pads outside the range [33;88] are written into the output file
898         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 2,
899                                sector, subSector, row);
900         subSector = 3;
901         //only the pads in the range [33;88] are written into the output file
902         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 1,
903                                  sector, subSector, row);
904       }
905     }//end else
906   }//end for
907
908   delete buffer;
909   fLoader->UnloadDigits();
910
911   AliTPCDDLRawData rawWriter;
912   rawWriter.SetVerbose(0);
913
914   rawWriter.RawData(fileName);
915   gSystem->Unlink(fileName);
916
917 }
918
919
920 //_____________________________________________________________________________
921 Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
922   // Converts the TPC raw data into summable digits
923   // The method is used for merging simulated and
924   // real data events
925   if (fLoader->TreeS() == 0x0 ) {
926     fLoader->MakeTree("S");
927   }
928
929   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
930
931   //setup TPCDigitsArray 
932   if(GetDigitsArray()) delete GetDigitsArray();
933
934   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
935   arr->SetClass("AliSimDigits");
936   arr->Setup(fTPCParam);
937   arr->MakeTree(fLoader->TreeS());
938
939   SetDigitsArray(arr);
940
941   // set zero suppression to "0"
942   fTPCParam->SetZeroSup(0);
943
944   // Loop over sectors
945   const Int_t kmaxTime = fTPCParam->GetMaxTBin();
946   const Int_t kNIS = fTPCParam->GetNInnerSector();
947   const Int_t kNOS = fTPCParam->GetNOuterSector();
948   const Int_t kNS = kNIS + kNOS;
949
950   // Setup storage
951   AliTPCROC * roc = AliTPCROC::Instance();
952   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
953   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
954   Short_t** allBins = new Short_t*[nRowsMax];
955   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
956     Int_t maxBin = kmaxTime*nPadsMax;
957     allBins[iRow] = new Short_t[maxBin];
958     memset(allBins[iRow],0,sizeof(Short_t)*maxBin);
959   }
960
961   for(Int_t iSector = 0; iSector < kNS; iSector++) {
962     
963     Int_t nRows = fTPCParam->GetNRow(iSector);
964     Int_t nDDLs = 0, indexDDL = 0;
965     if (iSector < kNIS) {
966       nDDLs = 2;
967       indexDDL = iSector * 2;
968     }
969     else {
970       nDDLs = 4;
971       indexDDL = (iSector-kNIS) * 4 + kNIS * 2;
972     }
973
974     // Load the raw data for corresponding DDLs
975     rawReader->Reset();
976
977     AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
978     AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping);
979     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
980
981     // Clean storage
982     for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
983       Int_t maxBin = kmaxTime*nPadsMax;
984       memset(allBins[iRow],0,sizeof(Short_t)*maxBin);
985     }
986
987     // Begin loop over altro data
988     while (input.NextDDL()) {
989
990       if (input.GetSector() != iSector)
991         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSector,input.GetSector()));
992
993       //loop over pads
994       while ( input.NextChannel() ) {
995
996         Int_t iRow = input.GetRow();
997         if (iRow < 0 || iRow >= nRows)
998           AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
999                         iRow, 0, nRows -1));
1000         Int_t iPad = input.GetPad();
1001
1002         Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
1003
1004         if (iPad < 0 || iPad >= maxPad)
1005           AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
1006                         iPad, 0, maxPad -1));
1007
1008         //loop over bunches
1009         while ( input.NextBunch() ){
1010           Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
1011           Int_t  bunchlength  = (Int_t)input.GetBunchLength();
1012           const UShort_t *sig = input.GetSignals();
1013           for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1014             Int_t iTimeBin=startTbin-iTime;
1015             if ( iTimeBin < 0 || iTimeBin >= kmaxTime) {
1016               continue;
1017               //AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1018               //               iTimeBin, 0, kmaxTime -1));
1019             }
1020
1021             Int_t maxBin = kmaxTime*maxPad;
1022             if (((iPad*kmaxTime+iTimeBin) >= maxBin) ||
1023                 ((iPad*kmaxTime+iTimeBin) < 0))
1024               AliFatal(Form("Index outside the allowed range"
1025                             " Sector=%d Row=%d Pad=%d Timebin=%d"
1026                             " (Max.index=%d)",iSector,iRow,iPad,iTimeBin,maxBin));
1027             allBins[iRow][iPad*kmaxTime+iTimeBin] = sig[iTime];
1028           }
1029         }
1030       } // End loop over altro data
1031     }
1032
1033     // Now fill the digits array
1034     if (fDigitsArray->GetTree()==0) {
1035       AliFatal("Tree not set in fDigitsArray");
1036     }
1037
1038     for (Int_t iRow = 0; iRow < nRows; iRow++) {
1039       AliDigits * dig = fDigitsArray->CreateRow(iSector,iRow);
1040       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
1041       for(Int_t iPad = 0; iPad < maxPad; iPad++) {
1042         for(Int_t iTimeBin = 0; iTimeBin < kmaxTime; iTimeBin++) {
1043           Short_t q = allBins[iRow][iPad*kmaxTime + iTimeBin];
1044           if (q <= 0) continue;
1045           q *= 16;
1046           dig->SetDigitFast((Short_t)q,iTimeBin,iPad);
1047         }
1048       }
1049       fDigitsArray->StoreRow(iSector,iRow);
1050       Int_t ndig = dig->GetDigitSize(); 
1051         
1052       AliDebug(10,
1053                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1054                     iSector,iRow,ndig));        
1055         
1056       fDigitsArray->ClearRow(iSector,iRow);  
1057
1058     } // end of the sector digitization
1059   }
1060   // get LHC clock phase from the digits tree
1061
1062   TParameter<float> *ph; 
1063   Float_t phase;
1064   TTree *digtree = fLoader->TreeD();
1065   //
1066   if(digtree){ // if TreeD exists
1067     ph = (TParameter<float>*)digtree->GetUserInfo()->FindObject("lhcphase0");
1068     phase = ph->GetVal();
1069   }
1070   else{ //TreeD does not exist
1071     phase = 0.; 
1072   }
1073     //
1074     // store lhc clock phase in S-digits tree
1075     //
1076     fLoader->TreeS()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",phase));
1077    //
1078    fLoader->WriteSDigits("OVERWRITE");
1079
1080   if(GetDigitsArray()) delete GetDigitsArray();
1081   SetDigitsArray(0x0);
1082
1083   // cleanup storage
1084   for (Int_t iRow = 0; iRow < nRowsMax; iRow++)
1085     delete [] allBins[iRow];
1086   delete [] allBins;
1087
1088   return kTRUE;
1089 }
1090
1091 //______________________________________________________________________
1092 AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
1093 {
1094   return new AliTPCDigitizer(manager);
1095 }
1096 //__
1097 void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)  
1098 {
1099   //create digits from summable digits
1100   GenerNoise(500000); //create teble with noise
1101
1102   //conect tree with sSDigits
1103   TTree *t = fLoader->TreeS();
1104
1105   if (t == 0x0) {
1106     fLoader->LoadSDigits("READ");
1107     t = fLoader->TreeS();
1108     if (t == 0x0) {
1109       AliError("Can not get input TreeS");
1110       return;
1111     }
1112   }
1113   
1114   if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D");
1115   
1116   AliSimDigits digarr, *dummy=&digarr;
1117   TBranch* sdb = t->GetBranch("Segment");
1118   if (sdb == 0x0) {
1119     AliError("Can not find branch with segments in TreeS.");
1120     return;
1121   }  
1122
1123   sdb->SetAddress(&dummy);
1124       
1125   Stat_t nentries = t->GetEntries();
1126
1127   // set zero suppression
1128
1129   fTPCParam->SetZeroSup(2);
1130
1131   // get zero suppression
1132
1133   Int_t zerosup = fTPCParam->GetZeroSup();
1134
1135   //make tree with digits 
1136   
1137   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1138   arr->SetClass("AliSimDigits");
1139   arr->Setup(fTPCParam);
1140   arr->MakeTree(fLoader->TreeD());
1141   
1142   AliTPCParam * par = fTPCParam;
1143
1144   //Loop over segments of the TPC
1145
1146   for (Int_t n=0; n<nentries; n++) {
1147     t->GetEvent(n);
1148     Int_t sec, row;
1149     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1150       AliWarning(Form("Invalid segment ID ! %d",digarr.GetID()));
1151       continue;
1152     }
1153     if (!IsSectorActive(sec)) continue;
1154     
1155     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1156     Int_t nrows = digrow->GetNRows();
1157     Int_t ncols = digrow->GetNCols();
1158
1159     digrow->ExpandBuffer();
1160     digarr.ExpandBuffer();
1161     digrow->ExpandTrackBuffer();
1162     digarr.ExpandTrackBuffer();
1163
1164     
1165     Short_t * pamp0 = digarr.GetDigits();
1166     Int_t   * ptracks0 = digarr.GetTracks();
1167     Short_t * pamp1 = digrow->GetDigits();
1168     Int_t   * ptracks1 = digrow->GetTracks();
1169     Int_t  nelems =nrows*ncols;
1170     Int_t saturation = fTPCParam->GetADCSat() - 1;
1171     //use internal structure of the AliDigits - for speed reason
1172     //if you cahnge implementation
1173     //of the Alidigits - it must be rewriten -
1174     for (Int_t i= 0; i<nelems; i++){
1175       Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
1176       if (q>zerosup){
1177         if (q>saturation) q=saturation;      
1178         *pamp1=(Short_t)q;
1179
1180         ptracks1[0]=ptracks0[0];        
1181         ptracks1[nelems]=ptracks0[nelems];
1182         ptracks1[2*nelems]=ptracks0[2*nelems];
1183       }
1184       pamp0++;
1185       pamp1++;
1186       ptracks0++;
1187       ptracks1++;        
1188     }
1189
1190     arr->StoreRow(sec,row);
1191     arr->ClearRow(sec,row);   
1192   }  
1193
1194     
1195   //write results
1196   fLoader->WriteDigits("OVERWRITE");
1197    
1198   delete arr;
1199 }
1200 //__________________________________________________________________
1201 void AliTPC::SetDefaults(){
1202   //
1203   // setting the defaults
1204   //
1205    
1206   // Set response functions
1207
1208   //
1209   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1210   rl->CdGAFile();
1211   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1212   // if(param){
1213 //     AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
1214 //     delete param;
1215 //     param = new AliTPCParamSR();
1216 //   }
1217 //   else {
1218 //     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1219 //   }
1220   param = (AliTPCParamSR*)AliTPCcalibDB::Instance()->GetParameters();
1221   if (!param->IsGeoRead()){
1222       //
1223       // read transformation matrices for gGeoManager
1224       //
1225       param->ReadGeoMatrices();
1226     }
1227   if(!param){
1228     AliFatal("No TPC parameters found");
1229   }
1230
1231
1232   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1233   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1234   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;  
1235
1236   
1237   //AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1238   //rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1239   //rf->SetOffset(3*param->GetZSigma());
1240   //rf->Update();
1241   //
1242   // Use gamma 4
1243   //
1244   char  strgamma4[1000];
1245   sprintf(strgamma4,"AliTPCRF1D::Gamma4((x-0.135+%f)*%f,55,160)",3*param->GetZSigma(), 1000000000*param->GetTSample()/param->GetZWidth());
1246   
1247   TF1 * fgamma4 = new TF1("fgamma4",strgamma4, -1,1);
1248   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE,1000);
1249   rf->SetParam(fgamma4,param->GetZWidth(), 1,0.2);
1250   rf->SetOffset(3*param->GetZSigma()); 
1251   rf->Update();
1252
1253   TDirectory *savedir=gDirectory;
1254   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1255   if (!f->IsOpen()) 
1256     AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !");
1257
1258   TString s;
1259   prfinner->Read("prf_07504_Gati_056068_d02");
1260   //PH Set different names
1261   s=prfinner->GetGRF()->GetName();
1262   s+="in";
1263   prfinner->GetGRF()->SetName(s.Data());
1264
1265   prfouter1->Read("prf_10006_Gati_047051_d03");
1266   s=prfouter1->GetGRF()->GetName();
1267   s+="out1";
1268   prfouter1->GetGRF()->SetName(s.Data());
1269
1270   prfouter2->Read("prf_15006_Gati_047051_d03");  
1271   s=prfouter2->GetGRF()->GetName();
1272   s+="out2";
1273   prfouter2->GetGRF()->SetName(s.Data());
1274
1275   f->Close();
1276   savedir->cd();
1277
1278   param->SetInnerPRF(prfinner);
1279   param->SetOuter1PRF(prfouter1); 
1280   param->SetOuter2PRF(prfouter2);
1281   param->SetTimeRF(rf);
1282
1283   // set fTPCParam
1284
1285   SetParam(param);
1286
1287
1288   fDefaults = 1;
1289
1290 }
1291 //__________________________________________________________________  
1292 void AliTPC::Hits2Digits()  
1293 {
1294   //
1295   // creates digits from hits
1296   //
1297   if (!fTPCParam->IsGeoRead()){
1298     //
1299     // read transformation matrices for gGeoManager
1300     //
1301     fTPCParam->ReadGeoMatrices();
1302   }
1303
1304   fLoader->LoadHits("read");
1305   fLoader->LoadDigits("recreate");
1306   AliRunLoader* runLoader = fLoader->GetRunLoader(); 
1307
1308   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1309     //PH    runLoader->GetEvent(iEvent);
1310     Hits2Digits(iEvent);
1311   }
1312
1313   fLoader->UnloadHits();
1314   fLoader->UnloadDigits();
1315
1316 //__________________________________________________________________  
1317 void AliTPC::Hits2Digits(Int_t eventnumber)  
1318
1319  //----------------------------------------------------
1320  // Loop over all sectors for a single event
1321  //----------------------------------------------------
1322   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1323   rl->GetEvent(eventnumber);
1324   SetActiveSectors();   
1325   if (fLoader->TreeH() == 0x0) {
1326     if(fLoader->LoadHits()) {
1327       AliError("Can not load hits.");
1328     }
1329   }
1330   SetTreeAddress();
1331   
1332   if (fLoader->TreeD() == 0x0 ) {
1333     fLoader->MakeTree("D");
1334     if (fLoader->TreeD() == 0x0 ) {
1335       AliError("Can not get TreeD");
1336       return;
1337     }
1338   }
1339
1340   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1341   GenerNoise(500000); //create teble with noise
1342
1343   //setup TPCDigitsArray 
1344
1345   if(GetDigitsArray()) delete GetDigitsArray();
1346   
1347   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1348   arr->SetClass("AliSimDigits");
1349   arr->Setup(fTPCParam);
1350
1351   arr->MakeTree(fLoader->TreeD());
1352   SetDigitsArray(arr);
1353
1354   fDigitsSwitch=0; // standard digits
1355   // here LHC clock phase
1356   Float_t lhcph = 0.;
1357   switch (fLHCclockPhaseSw){
1358   case 0: 
1359     // no phase
1360     lhcph=0.;
1361     break;
1362   case 1:
1363     // random phase
1364     lhcph = (Int_t)(gRandom->Rndm()/0.25);    
1365     break;
1366   case 2:
1367     lhcph=0.;
1368     // not implemented yet
1369     break;
1370   }
1371   // adding phase to the TreeD user info 
1372   fLoader->TreeD()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
1373   //
1374   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) 
1375     if (IsSectorActive(isec)) {
1376       AliDebug(1,Form("Hits2Digits: Sector %d is active.",isec));
1377       Hits2DigitsSector(isec);
1378     }
1379     else {
1380       AliDebug(1,Form("Hits2Digits: Sector %d is NOT active.",isec));
1381     }
1382   
1383   fLoader->WriteDigits("OVERWRITE"); 
1384   
1385 //this line prevents the crash in the similar one
1386 //on the beginning of this method
1387 //destructor attempts to reset the tree, which is deleted by the loader
1388 //need to be redesign
1389   if(GetDigitsArray()) delete GetDigitsArray();
1390   SetDigitsArray(0x0);
1391   
1392 }
1393
1394 //__________________________________________________________________
1395 void AliTPC::Hits2SDigits2(Int_t eventnumber)  
1396
1397
1398   //-----------------------------------------------------------
1399   //   summable digits - 16 bit "ADC", no noise, no saturation
1400   //-----------------------------------------------------------
1401
1402   //----------------------------------------------------
1403   // Loop over all sectors for a single event
1404   //----------------------------------------------------
1405
1406   AliRunLoader* rl = fLoader->GetRunLoader();
1407
1408   rl->GetEvent(eventnumber);
1409   if (fLoader->TreeH() == 0x0) {
1410     if(fLoader->LoadHits()) {
1411       AliError("Can not load hits.");
1412       return;
1413     }
1414   }
1415   SetTreeAddress();
1416
1417
1418   if (fLoader->TreeS() == 0x0 ) {
1419     fLoader->MakeTree("S");
1420   }
1421   
1422   if(fDefaults == 0) SetDefaults();
1423   
1424   GenerNoise(500000); //create table with noise
1425   //setup TPCDigitsArray 
1426
1427   if(GetDigitsArray()) delete GetDigitsArray();
1428
1429   
1430   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1431   arr->SetClass("AliSimDigits");
1432   arr->Setup(fTPCParam);
1433   arr->MakeTree(fLoader->TreeS());
1434
1435   SetDigitsArray(arr);
1436
1437   fDigitsSwitch=1; // summable digits
1438   
1439     // set zero suppression to "0"
1440   // here LHC clock phase
1441   Float_t lhcph = 0.;
1442   switch (fLHCclockPhaseSw){
1443   case 0: 
1444     // no phase
1445     lhcph=0.;
1446     break;
1447   case 1:
1448     // random phase
1449     lhcph = (Int_t)(gRandom->Rndm()/0.25);    
1450     break;
1451   case 2:
1452     lhcph=0.;
1453     // not implemented yet
1454     break;
1455   }
1456   // adding phase to the TreeS user info 
1457   
1458   fLoader->TreeS()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
1459
1460   fTPCParam->SetZeroSup(0);
1461
1462   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) 
1463     if (IsSectorActive(isec)) {
1464       Hits2DigitsSector(isec);
1465     }
1466
1467   fLoader->WriteSDigits("OVERWRITE");
1468
1469 //this line prevents the crash in the similar one
1470 //on the beginning of this method
1471 //destructor attempts to reset the tree, which is deleted by the loader
1472 //need to be redesign
1473   if(GetDigitsArray()) delete GetDigitsArray();
1474   SetDigitsArray(0x0);
1475 }
1476 //__________________________________________________________________
1477
1478 void AliTPC::Hits2SDigits()  
1479
1480
1481   //-----------------------------------------------------------
1482   //   summable digits - 16 bit "ADC", no noise, no saturation
1483   //-----------------------------------------------------------
1484
1485   if (!fTPCParam->IsGeoRead()){
1486     //
1487     // read transformation matrices for gGeoManager
1488     //
1489     fTPCParam->ReadGeoMatrices();
1490   }
1491   
1492   fLoader->LoadHits("read");
1493   fLoader->LoadSDigits("recreate");
1494   AliRunLoader* runLoader = fLoader->GetRunLoader(); 
1495
1496   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1497     runLoader->GetEvent(iEvent);
1498     SetTreeAddress();
1499     SetActiveSectors();
1500     Hits2SDigits2(iEvent);
1501   }
1502   
1503   fLoader->UnloadHits();
1504   fLoader->UnloadSDigits();
1505   if (fDebugStreamer) {
1506     delete fDebugStreamer;
1507     fDebugStreamer=0;
1508   }    
1509 }
1510 //_____________________________________________________________________________
1511
1512 void AliTPC::Hits2DigitsSector(Int_t isec)
1513 {
1514   //-------------------------------------------------------------------
1515   // TPC conversion from hits to digits.
1516   //------------------------------------------------------------------- 
1517
1518   //-----------------------------------------------------------------
1519   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1520   //-----------------------------------------------------------------
1521
1522   //-------------------------------------------------------
1523   //  Get the access to the track hits
1524   //-------------------------------------------------------
1525
1526   // check if the parameters are set - important if one calls this method
1527   // directly, not from the Hits2Digits
1528
1529   if(fDefaults == 0) SetDefaults();
1530
1531   TTree *tH = fLoader->TreeH(); // pointer to the hits tree
1532   if (tH == 0x0) {
1533     AliFatal("Can not find TreeH in folder");
1534     return;
1535   }
1536
1537   Stat_t ntracks = tH->GetEntries();
1538
1539
1540
1541     TObjArray **row;
1542     
1543     Int_t nrows =fTPCParam->GetNRow(isec);
1544
1545     row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1546     
1547     MakeSector(isec,nrows,tH,ntracks,row);
1548
1549     //--------------------------------------------------------
1550     //   Digitize this sector, row by row
1551     //   row[i] is the pointer to the TObjArray of TVectors,
1552     //   each one containing electrons accepted on this
1553     //   row, assigned into tracks
1554     //--------------------------------------------------------
1555
1556     Int_t i;
1557
1558     if (fDigitsArray->GetTree()==0) {
1559       AliFatal("Tree not set in fDigitsArray");
1560     }
1561
1562     for (i=0;i<nrows;i++){
1563       
1564       AliDigits * dig = fDigitsArray->CreateRow(isec,i); 
1565
1566       DigitizeRow(i,isec,row);
1567
1568       fDigitsArray->StoreRow(isec,i);
1569
1570       Int_t ndig = dig->GetDigitSize(); 
1571         
1572       AliDebug(10,
1573                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1574                     isec,i,ndig));        
1575         
1576       fDigitsArray->ClearRow(isec,i);  
1577
1578    
1579     } // end of the sector digitization
1580
1581     for(i=0;i<nrows+2;i++){
1582       row[i]->Delete();  
1583       delete row[i];   
1584     }
1585       
1586     delete [] row; // delete the array of pointers to TObjArray-s
1587         
1588
1589 } // end of Hits2DigitsSector
1590
1591
1592 //_____________________________________________________________________________
1593 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1594 {
1595   //-----------------------------------------------------------
1596   // Single row digitization, coupling from the neighbouring
1597   // rows taken into account
1598   //-----------------------------------------------------------
1599
1600   //-----------------------------------------------------------------
1601   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1602   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1603   //-----------------------------------------------------------------
1604  
1605   Float_t zerosup = fTPCParam->GetZeroSup();
1606   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor(); 
1607   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); 
1608   AliTPCCalROC * gainROC = gainTPC->GetCalROC(isec);  // pad gains per given sector
1609   AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(isec);  // noise per given sector
1610
1611
1612   fCurrentIndex[1]= isec;
1613   
1614
1615   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1616   Int_t nofTbins = fTPCParam->GetMaxTBin();
1617   Int_t indexRange[4];
1618   //
1619   //  Integrated signal for this row
1620   //  and a single track signal
1621   //    
1622
1623   TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated
1624   TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single
1625   //
1626   TMatrixF &total  = *m1;
1627
1628   //  Array of pointers to the label-signal list
1629
1630   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1631   Float_t  **pList = new Float_t* [nofDigits]; 
1632
1633   Int_t lp;
1634   Int_t i1;   
1635   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1636   //
1637   //calculate signal 
1638   //
1639   Int_t row1=irow;
1640   Int_t row2=irow+2; 
1641   for (Int_t row= row1;row<=row2;row++){
1642     Int_t nTracks= rows[row]->GetEntries();
1643     for (i1=0;i1<nTracks;i1++){
1644       fCurrentIndex[2]= row;
1645       fCurrentIndex[3]=irow+1;
1646       if (row==irow+1){
1647         m2->Zero();  // clear single track signal matrix
1648         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange); 
1649         GetList(trackLabel,nofPads,m2,indexRange,pList);
1650       }
1651       else   GetSignal(rows[row],i1,0,m1,indexRange);
1652     }
1653   }
1654          
1655   Int_t tracks[3];
1656
1657   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1658   Int_t gi=-1;
1659   Float_t fzerosup = zerosup+0.5;
1660   for(Int_t it=0;it<nofTbins;it++){
1661     for(Int_t ip=0;ip<nofPads;ip++){
1662       gi++;
1663       Float_t q=total(ip,it);      
1664       if(fDigitsSwitch == 0){   
1665         Float_t gain = gainROC->GetValue(irow,ip);  // get gain for given - pad-row pad 
1666         Float_t noisePad = noiseROC->GetValue(irow,ip); 
1667         //
1668         q*=gain;
1669         q+=GetNoise()*noisePad;
1670         if(q <=fzerosup) continue; // do not fill zeros
1671         q = TMath::Nint(q);
1672         if(q >= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1;  // saturation
1673
1674       }
1675
1676       else {
1677         if(q <= 0.) continue; // do not fill zeros
1678         if(q>2000.) q=2000.;
1679         q *= 16.;
1680         q = TMath::Nint(q);
1681       }
1682
1683       //
1684       //  "real" signal or electronic noise (list = -1)?
1685       //    
1686
1687       for(Int_t j1=0;j1<3;j1++){
1688         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1689       }
1690
1691 //Begin_Html
1692 /*
1693   <A NAME="AliDigits"></A>
1694   using of AliDigits object
1695 */
1696 //End_Html
1697       dig->SetDigitFast((Short_t)q,it,ip);
1698       if (fDigitsArray->IsSimulated()) {
1699         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1700         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1701         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1702       }
1703     
1704     } // end of loop over time buckets
1705   }  // end of lop over pads 
1706   //
1707   // test
1708   //
1709   //
1710
1711   // glitch filters if normal simulated digits
1712   //
1713   if(!fDigitsSwitch) ((AliSimDigits*)dig)->GlitchFilter();
1714   //
1715   //  This row has been digitized, delete nonused stuff
1716   //
1717
1718   for(lp=0;lp<nofDigits;lp++){
1719     if(pList[lp]) delete [] pList[lp];
1720   }
1721   
1722   delete [] pList;
1723
1724   delete m1;
1725   delete m2;
1726
1727 } // end of DigitizeRow
1728
1729 //_____________________________________________________________________________
1730
1731 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, 
1732              TMatrixF *m1, TMatrixF *m2,Int_t *indexRange)
1733 {
1734
1735   //---------------------------------------------------------------
1736   //  Calculates 2-D signal (pad,time) for a single track,
1737   //  returns a pointer to the signal matrix and the track label 
1738   //  No digitization is performed at this level!!!
1739   //---------------------------------------------------------------
1740
1741   //-----------------------------------------------------------------
1742   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1743   // Modified: Marian Ivanov 
1744   //-----------------------------------------------------------------
1745
1746   TVector *tv;
1747
1748   tv = (TVector*)p1->At(ntr); // pointer to a track
1749   TVector &v = *tv;
1750   
1751   Float_t label = v(0);
1752   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1))/2;
1753
1754   Int_t nElectrons = (tv->GetNrows()-1)/5;
1755   indexRange[0]=9999; // min pad
1756   indexRange[1]=-1; // max pad
1757   indexRange[2]=9999; //min time
1758   indexRange[3]=-1; // max time
1759
1760   TMatrixF &signal = *m1;
1761   TMatrixF &total = *m2;
1762   //
1763   // Get LHC clock phase
1764   //
1765   TParameter<float> *ph;
1766   if(fDigitsSwitch){// s-digits
1767     ph = (TParameter<float>*)fLoader->TreeS()->GetUserInfo()->FindObject("lhcphase0");  
1768   }
1769   else{ // normal digits
1770     ph = (TParameter<float>*)fLoader->TreeD()->GetUserInfo()->FindObject("lhcphase0");
1771   } 
1772   //  Loop over all electrons
1773   //
1774   for(Int_t nel=0; nel<nElectrons; nel++){
1775     Int_t idx=nel*5;
1776     Float_t aval =  v(idx+4);
1777     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac(); 
1778     Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)};
1779     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,
1780                                                             fCurrentIndex[3],ph->GetVal());
1781
1782     Int_t *index = fTPCParam->GetResBin(0);  
1783     Float_t *weight = & (fTPCParam->GetResWeight(0));
1784
1785     if (n>0) for (Int_t i =0; i<n; i++){       
1786       Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1787
1788       if (pad>=0){
1789         Int_t time=index[2];     
1790         Float_t qweight = *(weight)*eltoadcfac;
1791         
1792         if (m1!=0) signal(pad,time)+=qweight;
1793         total(pad,time)+=qweight;
1794         if (indexRange[0]>pad) indexRange[0]=pad;
1795         if (indexRange[1]<pad) indexRange[1]=pad;
1796         if (indexRange[2]>time) indexRange[2]=time;
1797         if (indexRange[3]<time) indexRange[3]=time;
1798         
1799         index+=3;
1800         weight++;       
1801
1802       }  
1803     }
1804   } // end of loop over electrons
1805   
1806   return label; // returns track label when finished
1807 }
1808
1809 //_____________________________________________________________________________
1810 void AliTPC::GetList(Float_t label,Int_t np,TMatrixF *m,
1811                      Int_t *indexRange, Float_t **pList)
1812 {
1813   //----------------------------------------------------------------------
1814   //  Updates the list of tracks contributing to digits for a given row
1815   //----------------------------------------------------------------------
1816
1817   //-----------------------------------------------------------------
1818   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1819   //-----------------------------------------------------------------
1820
1821   TMatrixF &signal = *m;
1822
1823   // lop over nonzero digits
1824
1825   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1826     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1827
1828
1829       // accept only the contribution larger than 500 electrons (1/2 s_noise)
1830
1831       if(signal(ip,it)<0.5) continue; 
1832
1833       Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1834         
1835       if(!pList[globalIndex]){
1836         
1837         // 
1838         // Create new list (6 elements - 3 signals and 3 labels),
1839         //
1840
1841         pList[globalIndex] = new Float_t [6];
1842
1843         // set list to -1 
1844         
1845         *pList[globalIndex] = -1.;
1846         *(pList[globalIndex]+1) = -1.;
1847         *(pList[globalIndex]+2) = -1.;
1848         *(pList[globalIndex]+3) = -1.;
1849         *(pList[globalIndex]+4) = -1.;
1850         *(pList[globalIndex]+5) = -1.;
1851
1852         *pList[globalIndex] = label;
1853         *(pList[globalIndex]+3) = signal(ip,it);
1854       }
1855       else {
1856
1857         // check the signal magnitude
1858
1859         Float_t highest = *(pList[globalIndex]+3);
1860         Float_t middle = *(pList[globalIndex]+4);
1861         Float_t lowest = *(pList[globalIndex]+5);
1862         
1863         //
1864         //  compare the new signal with already existing list
1865         //
1866         
1867         if(signal(ip,it)<lowest) continue; // neglect this track
1868
1869         //
1870
1871         if (signal(ip,it)>highest){
1872           *(pList[globalIndex]+5) = middle;
1873           *(pList[globalIndex]+4) = highest;
1874           *(pList[globalIndex]+3) = signal(ip,it);
1875           
1876           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1877           *(pList[globalIndex]+1) = *pList[globalIndex];
1878           *pList[globalIndex] = label;
1879         }
1880         else if (signal(ip,it)>middle){
1881           *(pList[globalIndex]+5) = middle;
1882           *(pList[globalIndex]+4) = signal(ip,it);
1883           
1884           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1885           *(pList[globalIndex]+1) = label;
1886         }
1887         else{
1888           *(pList[globalIndex]+5) = signal(ip,it);
1889           *(pList[globalIndex]+2) = label;
1890         }
1891       }
1892       
1893     } // end of loop over pads
1894   } // end of loop over time bins
1895
1896 }//end of GetList
1897 //___________________________________________________________________
1898 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1899                         Stat_t ntracks,TObjArray **row)
1900 {
1901
1902   //-----------------------------------------------------------------
1903   // Prepares the sector digitization, creates the vectors of
1904   // tracks for each row of this sector. The track vector
1905   // contains the track label and the position of electrons.
1906   //-----------------------------------------------------------------
1907
1908   // 
1909   // The trasport of the electrons through TPC drift volume
1910   //    Drift (drift velocity + velocity map(not yet implemented)))
1911   //    Application of the random processes (diffusion, gas gain)
1912   //    Systematic effects (ExB effect in drift volume + ROCs)  
1913   //
1914   // Algorithm:
1915   // Loop over primary electrons:
1916   //    Creation of the secondary electrons
1917   //    Loop over electrons (primary+ secondaries)
1918   //        Global coordinate frame:
1919   //          1. Skip electrons if attached  
1920   //          2. ExB effect in drift volume
1921   //             a.) Simulation   calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1922   //             b.) Reconstruction -  calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1923   //          3. Generation of gas gain (Random - Exponential distribution) 
1924   //          4. TransportElectron function (diffusion)
1925   //
1926   //        5. Conversion to the local coordinate frame  pad-row, pad, timebin
1927   //        6. Apply Time0 shift - AliTPCCalPad class 
1928   //            a.) Plus sign in simulation
1929   //            b.) Minus sign in reconstruction 
1930   // end of loop          
1931   //
1932   //-----------------------------------------------------------------
1933   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1934   // Origin: Marian Ivanov,  marian.ivanov@cern.ch
1935   //-----------------------------------------------------------------
1936   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
1937   if (gAlice){ // Set correctly the magnetic field in the ExB calculation
1938     if (!calib->GetExB()){
1939       AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField()); 
1940       if (field) {
1941         calib->SetExBField(field);
1942       }
1943     }
1944   }
1945
1946   Float_t gasgain = fTPCParam->GetGasGain();
1947   gasgain = gasgain/fGainFactor;
1948   Int_t i;
1949   Float_t xyz[5]; 
1950
1951   AliTPChit *tpcHit; // pointer to a sigle TPC hit    
1952   //MI change
1953   TBranch * branch=0;
1954   if (fHitType>1) branch = TH->GetBranch("TPC2");
1955   else branch = TH->GetBranch("TPC");
1956
1957  
1958   //----------------------------------------------
1959   // Create TObjArray-s, one for each row,
1960   // each TObjArray will store the TVectors
1961   // of electrons, one TVectors per each track.
1962   //---------------------------------------------- 
1963     
1964   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
1965   TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
1966
1967   for(i=0; i<nrows+2; i++){
1968     row[i] = new TObjArray;
1969     nofElectrons[i]=0;
1970     tracks[i]=0;
1971   }
1972
1973  
1974
1975   //--------------------------------------------------------------------
1976   //  Loop over tracks, the "track" contains the full history
1977   //--------------------------------------------------------------------
1978   
1979   Int_t previousTrack,currentTrack;
1980   previousTrack = -1; // nothing to store so far!
1981
1982   for(Int_t track=0;track<ntracks;track++){
1983     Bool_t isInSector=kTRUE;
1984     ResetHits();
1985     isInSector = TrackInVolume(isec,track);
1986     if (!isInSector) continue;
1987     //MI change
1988     branch->GetEntry(track); // get next track
1989     
1990     //M.I. changes
1991
1992     tpcHit = (AliTPChit*)FirstHit(-1);
1993
1994     //--------------------------------------------------------------
1995     //  Loop over hits
1996     //--------------------------------------------------------------
1997
1998
1999     while(tpcHit){
2000       
2001       Int_t sector=tpcHit->fSector; // sector number
2002       if(sector != isec){
2003         tpcHit = (AliTPChit*) NextHit();
2004         continue; 
2005       }
2006
2007       // Remove hits which arrive before the TPC opening gate signal
2008       if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
2009           /fTPCParam->GetDriftV()+tpcHit->Time())<fTPCParam->GetGateDelay()) {
2010         tpcHit = (AliTPChit*) NextHit();
2011         continue;
2012       }
2013
2014       currentTrack = tpcHit->Track(); // track number
2015
2016       if(currentTrack != previousTrack){
2017                           
2018         // store already filled fTrack
2019               
2020         for(i=0;i<nrows+2;i++){
2021           if(previousTrack != -1){
2022             if(nofElectrons[i]>0){
2023               TVector &v = *tracks[i];
2024               v(0) = previousTrack;
2025               tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2026               row[i]->Add(tracks[i]);                     
2027             }
2028             else {
2029               delete tracks[i]; // delete empty TVector
2030               tracks[i]=0;
2031             }
2032           }
2033
2034           nofElectrons[i]=0;
2035           tracks[i] = new TVector(601); // TVectors for the next fTrack
2036
2037         } // end of loop over rows
2038                
2039         previousTrack=currentTrack; // update track label 
2040       }
2041            
2042       Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2043
2044       //---------------------------------------------------
2045       //  Calculate the electron attachment probability
2046       //---------------------------------------------------
2047
2048
2049       Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
2050         /fTPCParam->GetDriftV(); 
2051       // in microseconds!       
2052       Float_t attProb = fTPCParam->GetAttCoef()*
2053         fTPCParam->GetOxyCont()*time; //  fraction! 
2054    
2055       //-----------------------------------------------
2056       //  Loop over electrons
2057       //-----------------------------------------------
2058       Int_t index[3];
2059       index[1]=isec;
2060       for(Int_t nel=0;nel<qI;nel++){
2061         // skip if electron lost due to the attachment
2062         if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
2063         
2064         //
2065         // ExB effect
2066         //
2067         Double_t dxyz0[3],dxyz1[3];
2068         dxyz0[0]=tpcHit->X();
2069         dxyz0[1]=tpcHit->Y();
2070         dxyz0[2]=tpcHit->Z();   
2071         if (calib->GetExB()){
2072           calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
2073         }else{
2074           AliError("Not valid ExB calibration");
2075           dxyz1[0]=tpcHit->X();
2076           dxyz1[1]=tpcHit->Y();
2077           dxyz1[2]=tpcHit->Z();         
2078         }
2079         xyz[0]=dxyz1[0];
2080         xyz[1]=dxyz1[1];
2081         xyz[2]=dxyz1[2];        
2082         //
2083         //
2084         //
2085         // protection for the nonphysical avalanche size (10**6 maximum)
2086         //  
2087         Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2088         xyz[3]= (Float_t) (-gasgain*TMath::Log(rn)); 
2089         index[0]=1;
2090           
2091         TransportElectron(xyz,index);    
2092         Int_t rowNumber;
2093         Int_t padrow = fTPCParam->GetPadRow(xyz,index); 
2094         //
2095         // Add Time0 correction due unisochronity
2096         // xyz[0] - pad row coordinate 
2097         // xyz[1] - pad coordinate
2098         // xyz[2] - is in now time bin coordinate system
2099         Float_t correction =0;
2100         if (calib->GetPadTime0()){
2101           if (!calib->GetPadTime0()->GetCalROC(isec)) continue;   
2102           Int_t npads = fTPCParam->GetNPads(isec,padrow);
2103           //      Int_t pad  = TMath::Nint(xyz[1]+fTPCParam->GetNPads(isec,TMath::Nint(xyz[0]))*0.5);
2104           // pad numbering from -npads/2 .. npads/2-1
2105           Int_t pad  = TMath::Nint(xyz[1]+npads/2);
2106           if (pad<0) pad=0;
2107           if (pad>=npads) pad=npads-1;
2108           correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(padrow,pad);
2109           //      printf("%d\t%d\t%d\t%f\n",isec,padrow,pad,correction);
2110           if (fDebugStreamer){
2111             (*fDebugStreamer)<<"Time0"<<
2112               "isec="<<isec<<
2113               "padrow="<<padrow<<
2114               "pad="<<pad<<
2115               "x0="<<xyz[0]<<
2116               "x1="<<xyz[1]<<
2117               "x2="<<xyz[2]<<
2118               "hit.="<<tpcHit<<
2119               "cor="<<correction<<
2120               "\n";
2121           }
2122         }
2123         xyz[2]+=correction;
2124         xyz[2]+=fTPCParam->GetNTBinsL1();    // adding Level 1 time bin offset
2125         //
2126         // Electron track time (for pileup simulation)
2127         xyz[2]+=tpcHit->Time()/fTPCParam->GetTSample(); // adding time of flight
2128         xyz[4] =0;
2129
2130         //
2131         // row 0 - cross talk from the innermost row
2132         // row fNRow+1 cross talk from the outermost row
2133         rowNumber = index[2]+1; 
2134         //transform position to local digit coordinates
2135         //relative to nearest pad row 
2136         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
2137         /*      Float_t x1,y1;
2138         if (isec <fTPCParam->GetNInnerSector()) {
2139           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
2140           y1 = fTPCParam->GetYInner(rowNumber);
2141         }
2142         else{
2143           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
2144           y1 = fTPCParam->GetYOuter(rowNumber);
2145         }
2146         // gain inefficiency at the wires edges - linear
2147         x1=TMath::Abs(x1);
2148         y1-=1.;
2149         if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.)); */
2150         
2151         nofElectrons[rowNumber]++;        
2152         //----------------------------------
2153         // Expand vector if necessary
2154         //----------------------------------
2155         if(nofElectrons[rowNumber]>120){
2156           Int_t range = tracks[rowNumber]->GetNrows();
2157           if((nofElectrons[rowNumber])>(range-1)/5){
2158             
2159             tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons
2160           }
2161         }
2162         
2163         TVector &v = *tracks[rowNumber];
2164         Int_t idx = 5*nofElectrons[rowNumber]-4;
2165         Real_t * position = &(((TVector&)v)(idx)); //make code faster
2166         memcpy(position,xyz,5*sizeof(Float_t));
2167         
2168       } // end of loop over electrons
2169
2170       tpcHit = (AliTPChit*)NextHit();
2171       
2172     } // end of loop over hits
2173   } // end of loop over tracks
2174
2175     //
2176     //   store remaining track (the last one) if not empty
2177     //
2178   
2179   for(i=0;i<nrows+2;i++){
2180     if(nofElectrons[i]>0){
2181       TVector &v = *tracks[i];
2182       v(0) = previousTrack;
2183       tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2184       row[i]->Add(tracks[i]);  
2185     }
2186     else{
2187       delete tracks[i];
2188       tracks[i]=0;
2189     }  
2190   }  
2191   
2192   delete [] tracks;
2193   delete [] nofElectrons;
2194
2195 } // end of MakeSector
2196
2197
2198 //_____________________________________________________________________________
2199 void AliTPC::Init()
2200 {
2201   //
2202   // Initialise TPC detector after definition of geometry
2203   //
2204   AliDebug(1,"*********************************************");
2205 }
2206
2207 //_____________________________________________________________________________
2208 void AliTPC::ResetDigits()
2209 {
2210   //
2211   // Reset number of digits and the digits array for this detector
2212   //
2213   fNdigits   = 0;
2214   if (fDigits)   fDigits->Clear();
2215 }
2216
2217
2218
2219 //_____________________________________________________________________________
2220 void AliTPC::SetSens(Int_t sens)
2221 {
2222
2223   //-------------------------------------------------------------
2224   // Activates/deactivates the sensitive strips at the center of
2225   // the pad row -- this is for the space-point resolution calculations
2226   //-------------------------------------------------------------
2227
2228   //-----------------------------------------------------------------
2229   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2230   //-----------------------------------------------------------------
2231
2232   fSens = sens;
2233 }
2234
2235  
2236 void AliTPC::SetSide(Float_t side=0.)
2237 {
2238   // choice of the TPC side
2239
2240   fSide = side;
2241  
2242 }
2243 //_____________________________________________________________________________
2244
2245 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2246 {
2247   //
2248   // electron transport taking into account:
2249   // 1. diffusion, 
2250   // 2.ExB at the wires
2251   // 3. nonisochronity
2252   //
2253   // xyz and index must be already transformed to system 1
2254   //
2255
2256   fTPCParam->Transform1to2(xyz,index);  // mis-alignment applied in this step
2257   
2258   //add diffusion
2259   Float_t driftl=xyz[2];
2260   if(driftl<0.01) driftl=0.01;
2261   driftl=TMath::Sqrt(driftl);
2262   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2263   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2264   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2265   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2266   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2267
2268   // ExB
2269   
2270   if (fTPCParam->GetMWPCReadout()==kTRUE){
2271     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2272     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2273   }
2274   //add nonisochronity (not implemented yet) 
2275  
2276   
2277 }
2278   
2279 ClassImp(AliTPChit)
2280   //______________________________________________________________________
2281   AliTPChit::AliTPChit()
2282             :AliHit(),
2283              fSector(0),
2284              fPadRow(0),
2285              fQ(0),
2286              fTime(0)
2287 {
2288   //
2289   // default
2290   //
2291
2292 }
2293 //_____________________________________________________________________________
2294 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
2295           :AliHit(shunt,track),
2296              fSector(0),
2297              fPadRow(0),
2298              fQ(0),
2299              fTime(0)
2300 {
2301   //
2302   // Creates a TPC hit object
2303   //
2304   fSector     = vol[0];
2305   fPadRow     = vol[1];
2306   fX          = hits[0];
2307   fY          = hits[1];
2308   fZ          = hits[2];
2309   fQ          = hits[3];
2310   fTime       = hits[4];
2311 }
2312  
2313 //________________________________________________________________________
2314 // Additional code because of the AliTPCTrackHitsV2
2315
2316 void AliTPC::MakeBranch(Option_t *option)
2317 {
2318   //
2319   // Create a new branch in the current Root Tree
2320   // The branch of fHits is automatically split
2321   // MI change 14.09.2000
2322   AliDebug(1,"");
2323   if (fHitType<2) return;
2324   char branchname[10];
2325   sprintf(branchname,"%s2",GetName());  
2326   //
2327   // Get the pointer to the header
2328   const char *cH = strstr(option,"H");
2329   //
2330   if (fTrackHits   && fLoader->TreeH() && cH && fHitType&4) {
2331     AliDebug(1,"Making branch for Type 4 Hits");
2332     fLoader->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2333   }
2334
2335 //   if (fTrackHitsOld   && fLoader->TreeH() && cH && fHitType&2) {    
2336 //     AliDebug(1,"Making branch for Type 2 Hits");
2337 //     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld, 
2338 //                                                    fLoader->TreeH(),fBufferSize,99);
2339 //     fLoader->TreeH()->GetListOfBranches()->Add(branch);
2340 //   }  
2341 }
2342
2343 void AliTPC::SetTreeAddress()
2344 {
2345   //Sets tree address for hits  
2346   if (fHitType<=1) {
2347     if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2348     AliDetector::SetTreeAddress();
2349   }
2350   if (fHitType>1) SetTreeAddress2();
2351 }
2352
2353 void AliTPC::SetTreeAddress2()
2354 {
2355   //
2356   // Set branch address for the TrackHits Tree
2357   // 
2358   AliDebug(1,"");
2359   
2360   TBranch *branch;
2361   char branchname[20];
2362   sprintf(branchname,"%s2",GetName());
2363   //
2364   // Branch address for hit tree
2365   TTree *treeH = fLoader->TreeH();
2366   if ((treeH)&&(fHitType&4)) {
2367     branch = treeH->GetBranch(branchname);
2368     if (branch) {
2369       branch->SetAddress(&fTrackHits);
2370       AliDebug(1,"fHitType&4 Setting");
2371     }
2372     else 
2373       AliDebug(1,"fHitType&4 Failed (can not find branch)");
2374     
2375   }
2376  //  if ((treeH)&&(fHitType&2)) {
2377 //     branch = treeH->GetBranch(branchname);
2378 //     if (branch) {
2379 //       branch->SetAddress(&fTrackHitsOld);
2380 //       AliDebug(1,"fHitType&2 Setting");
2381 //     }
2382 //     else
2383 //       AliDebug(1,"fHitType&2 Failed (can not find branch)");
2384 //   }
2385 }
2386
2387 void AliTPC::FinishPrimary()
2388 {
2389   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();  
2390   //  if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();  
2391 }
2392
2393
2394 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2395
2396   //
2397   // add hit to the list
2398
2399   Int_t rtrack;
2400   if (fIshunt) {
2401     int primary = gAlice->GetMCApp()->GetPrimary(track);
2402     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2403     rtrack=primary;
2404   } else {
2405     rtrack=track;
2406     gAlice->GetMCApp()->FlagTrack(track);
2407   }  
2408   if (fTrackHits && fHitType&4) 
2409     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2410                              hits[1],hits[2],(Int_t)hits[3],hits[4]);
2411  //  if (fTrackHitsOld &&fHitType&2 ) 
2412 //     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2413 //                                 hits[1],hits[2],(Int_t)hits[3]);
2414   
2415 }
2416
2417 void AliTPC::ResetHits()
2418
2419   if (fHitType&1) AliDetector::ResetHits();
2420   if (fHitType>1) ResetHits2();
2421 }
2422
2423 void AliTPC::ResetHits2()
2424 {
2425   //
2426   //reset hits
2427   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2428   // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2429
2430 }   
2431
2432 AliHit* AliTPC::FirstHit(Int_t track)
2433 {
2434   if (fHitType>1) return FirstHit2(track);
2435   return AliDetector::FirstHit(track);
2436 }
2437 AliHit* AliTPC::NextHit()
2438 {
2439   //
2440   // gets next hit
2441   //
2442   if (fHitType>1) return NextHit2();
2443   
2444   return AliDetector::NextHit();
2445 }
2446
2447 AliHit* AliTPC::FirstHit2(Int_t track)
2448 {
2449   //
2450   // Initialise the hit iterator
2451   // Return the address of the first hit for track
2452   // If track>=0 the track is read from disk
2453   // while if track<0 the first hit of the current
2454   // track is returned
2455   // 
2456   if(track>=0) {
2457     gAlice->GetMCApp()->ResetHits();
2458     fLoader->TreeH()->GetEvent(track);
2459   }
2460   //
2461   if (fTrackHits && fHitType&4) {
2462     fTrackHits->First();
2463     return fTrackHits->GetHit();
2464   }
2465  //  if (fTrackHitsOld && fHitType&2) {
2466 //     fTrackHitsOld->First();
2467 //     return fTrackHitsOld->GetHit();
2468 //   }
2469
2470   else return 0;
2471 }
2472
2473 AliHit* AliTPC::NextHit2()
2474 {
2475   //
2476   //Return the next hit for the current track
2477
2478
2479 //   if (fTrackHitsOld && fHitType&2) {
2480 //     fTrackHitsOld->Next();
2481 //     return fTrackHitsOld->GetHit();
2482 //   }
2483   if (fTrackHits) {
2484     fTrackHits->Next();
2485     return fTrackHits->GetHit();
2486   }
2487   else 
2488     return 0;
2489 }
2490
2491 void AliTPC::RemapTrackHitIDs(Int_t *map)
2492 {
2493   //
2494   // remapping
2495   //
2496   if (!fTrackHits) return;
2497   
2498 //   if (fTrackHitsOld && fHitType&2){
2499 //     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2500 //     for (UInt_t i=0;i<arr->GetSize();i++){
2501 //       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2502 //       info->fTrackID = map[info->fTrackID];
2503 //     }
2504 //   }
2505 //  if (fTrackHitsOld && fHitType&4){
2506   if (fTrackHits && fHitType&4){
2507     TClonesArray * arr = fTrackHits->GetArray();;
2508     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2509       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2510       info->SetTrackID(map[info->GetTrackID()]);
2511     }
2512   }
2513 }
2514
2515 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2516 {
2517   //return bool information - is track in given volume
2518   //load only part of the track information 
2519   //return true if current track is in volume
2520   //
2521   //  return kTRUE;
2522  //  if (fTrackHitsOld && fHitType&2) {
2523 //     TBranch * br = fLoader->TreeH()->GetBranch("fTrackHitsInfo");
2524 //     br->GetEvent(track);
2525 //     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2526 //     for (UInt_t j=0;j<ar->GetSize();j++){
2527 //       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2528 //     } 
2529 //   }
2530
2531   if (fTrackHits && fHitType&4) {
2532     TBranch * br1 = fLoader->TreeH()->GetBranch("fVolumes");
2533     TBranch * br2 = fLoader->TreeH()->GetBranch("fNVolumes");    
2534     br2->GetEvent(track);
2535     br1->GetEvent(track);    
2536     Int_t *volumes = fTrackHits->GetVolumes();
2537     Int_t nvolumes = fTrackHits->GetNVolumes();
2538     if (!volumes && nvolumes>0) {
2539       AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes));
2540       return kFALSE;
2541     }
2542     for (Int_t j=0;j<nvolumes; j++)
2543       if (volumes[j]==id) return kTRUE;    
2544   }
2545
2546   if (fHitType&1) {
2547     TBranch * br = fLoader->TreeH()->GetBranch("fSector");
2548     br->GetEvent(track);
2549     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2550       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2551     } 
2552   }
2553   return kFALSE;  
2554
2555 }
2556
2557
2558 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
2559 {
2560   //Makes TPC loader
2561   fLoader = new AliTPCLoader(GetName(),topfoldername);
2562   return fLoader;
2563 }
2564
2565 ////////////////////////////////////////////////////////////////////////
2566 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
2567 //
2568 // load TPC paarmeters from a given file or create new if the object
2569 // is not found there
2570 // 12/05/2003 This method should be moved to the AliTPCLoader
2571 // and one has to decide where to store the TPC parameters
2572 // M.Kowalski
2573   char paramName[50];
2574   sprintf(paramName,"75x40_100x60_150x60");
2575   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2576   if (paramTPC) {
2577     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
2578   } else {
2579     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
2580     //paramTPC = new AliTPCParamSR;
2581     paramTPC = AliTPCcalibDB::Instance()->GetParameters();
2582     if (!paramTPC->IsGeoRead()){
2583       //
2584       // read transformation matrices for gGeoManager
2585       //
2586       paramTPC->ReadGeoMatrices();
2587     }
2588   
2589   }
2590   return paramTPC;
2591
2592 // the older version of parameters can be accessed with this code.
2593 // In some cases, we have old parameters saved in the file but 
2594 // digits were created with new parameters, it can be distinguish 
2595 // by the name of TPC TreeD. The code here is just for the case 
2596 // we would need to compare with old data, uncomment it if needed.
2597 //
2598 //  char paramName[50];
2599 //  sprintf(paramName,"75x40_100x60");
2600 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2601 //  if (paramTPC) {
2602 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
2603 //  } else {
2604 //    sprintf(paramName,"75x40_100x60_150x60");
2605 //    paramTPC=(AliTPCParam*)in->Get(paramName);
2606 //    if (paramTPC) {
2607 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
2608 //    } else {
2609 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2610 //          <<endl;    
2611 //      paramTPC = new AliTPCParamSR;
2612 //    }
2613 //  }
2614 //  return paramTPC;
2615
2616 }
2617
2618