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