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