Noise table created in Hits2SDigits (M.Ivanov)
[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 /*
17 $Log$
18 Revision 1.48  2001/11/24 16:10:21  kowal2
19 Faster algorithms.
20
21 Revision 1.47  2001/11/19 10:25:34  kowal2
22 Nearest integer instead of integer when converting to ADC counts
23
24 Revision 1.46  2001/11/07 06:47:12  kowal2
25 Removed printouts
26
27 Revision 1.45  2001/11/03 13:33:48  kowal2
28 Updated algorithms in Hits2SDigits, SDigits2Digits,
29 Hits2ExactClusters.
30 Added method Merge
31
32 Revision 1.44  2001/08/30 09:28:48  hristov
33 TTree names are explicitly set via SetName(name) and then Write() is called
34
35 Revision 1.43  2001/07/28 12:02:54  hristov
36 Branch split level set to 99
37
38 Revision 1.42  2001/07/28 11:38:52  hristov
39 Loop variable declared once
40
41 Revision 1.41  2001/07/28 10:53:50  hristov
42 Digitisation done according to the general scheme (M.Ivanov)
43
44 Revision 1.40  2001/07/27 13:03:14  hristov
45 Default Branch split level set to 99
46
47 Revision 1.39  2001/07/26 09:09:34  kowal2
48 Hits2Reco method added
49
50 Revision 1.38  2001/07/20 14:32:43  kowal2
51 Processing of many events possible now
52
53 Revision 1.37  2001/06/12 07:17:18  kowal2
54 Hits2SDigits method implemented (summable digits)
55
56 Revision 1.36  2001/05/16 14:57:25  alibrary
57 New files for folders and Stack
58
59 Revision 1.35  2001/05/08 16:02:22  kowal2
60 Updated material specifications
61
62 Revision 1.34  2001/05/08 15:00:15  hristov
63 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
64
65 Revision 1.33  2001/04/03 12:40:43  kowal2
66 Removed printouts
67
68 Revision 1.32  2001/03/12 17:47:36  hristov
69 Changes needed on Sun with CC 5.0
70
71 Revision 1.31  2001/03/12 08:21:50  kowal2
72 Corrected C++ bug in the material definitions
73
74 Revision 1.30  2001/03/01 17:34:47  kowal2
75 Correction due to the accuracy problem
76
77 Revision 1.29  2001/02/28 16:34:40  kowal2
78 Protection against nonphysical values of the avalanche size,
79 10**6 is the maximum
80
81 Revision 1.28  2001/01/26 19:57:19  hristov
82 Major upgrade of AliRoot code
83
84 Revision 1.27  2001/01/13 17:29:33  kowal2
85 Sun compiler correction
86
87 Revision 1.26  2001/01/10 07:59:43  kowal2
88 Corrections to load points from the noncompressed hits.
89
90 Revision 1.25  2000/11/02 07:25:31  kowal2
91 Changes due to the new hit structure.
92 Memory leak removed.
93
94 Revision 1.24  2000/10/05 16:06:09  kowal2
95 Forward declarations. Changes due to a new class AliComplexCluster.
96
97 Revision 1.23  2000/10/02 21:28:18  fca
98 Removal of useless dependecies via forward declarations
99
100 Revision 1.22  2000/07/10 20:57:39  hristov
101 Update of TPC code and macros by M.Kowalski
102
103 Revision 1.19.2.4  2000/06/26 07:39:42  kowal2
104 Changes to obey the coding rules
105
106 Revision 1.19.2.3  2000/06/25 08:38:41  kowal2
107 Splitted from AliTPCtracking
108
109 Revision 1.19.2.2  2000/06/16 12:59:28  kowal2
110 Changed parameter settings
111
112 Revision 1.19.2.1  2000/06/09 07:15:07  kowal2
113
114 Defaults loaded automatically (hard-wired)
115 Optional parameters can be set via macro called in the constructor
116
117 Revision 1.19  2000/04/18 19:00:59  fca
118 Small bug fixes to TPC files
119
120 Revision 1.18  2000/04/17 09:37:33  kowal2
121 removed obsolete AliTPCDigitsDisplay.C
122
123 Revision 1.17.2.2  2000/04/10 08:15:12  kowal2
124
125 New, experimental data structure from M. Ivanov
126 New tracking algorithm
127 Different pad geometry for different sectors
128 Digitization rewritten
129
130 Revision 1.17.2.1  2000/04/10 07:56:53  kowal2
131 Not used anymore - removed
132
133 Revision 1.17  2000/01/19 17:17:30  fca
134 Introducing a list of lists of hits -- more hits allowed for detector now
135
136 Revision 1.16  1999/11/05 09:29:23  fca
137 Accept only signals > 0
138
139 Revision 1.15  1999/10/08 06:26:53  fca
140 Removed ClustersIndex - not used anymore
141
142 Revision 1.14  1999/09/29 09:24:33  fca
143 Introduction of the Copyright and cvs Log
144
145 */
146
147 ///////////////////////////////////////////////////////////////////////////////
148 //                                                                           //
149 //  Time Projection Chamber                                                  //
150 //  This class contains the basic functions for the Time Projection Chamber  //
151 //  detector. Functions specific to one particular geometry are              //
152 //  contained in the derived classes                                         //
153 //                                                                           //
154 //Begin_Html
155 /*
156 <img src="picts/AliTPCClass.gif">
157 */
158 //End_Html
159 //                                                                           //
160 //                                                                          //
161 ///////////////////////////////////////////////////////////////////////////////
162
163 //
164
165 #include <TMath.h>
166 #include <TRandom.h>
167 #include <TVector.h>
168 #include <TMatrix.h>
169 #include <TGeometry.h>
170 #include <TNode.h>
171 #include <TTUBS.h>
172 #include <TObjectTable.h>
173 #include "TParticle.h"
174 #include "AliTPC.h"
175 #include <TFile.h>  
176 #include <TROOT.h>
177 #include <TSystem.h>     
178 #include "AliRun.h"
179 #include <iostream.h>
180 #include <stdlib.h>
181 #include <fstream.h>
182 #include "AliMC.h"
183 #include "AliMagF.h"
184
185
186 #include "AliTPCParamSR.h"
187 #include "AliTPCPRF2D.h"
188 #include "AliTPCRF1D.h"
189 #include "AliDigits.h"
190 #include "AliSimDigits.h"
191 #include "AliTPCTrackHits.h"
192 #include "AliPoints.h"
193 #include "AliArrayBranch.h"
194
195
196 #include "AliTPCDigitsArray.h"
197 #include "AliComplexCluster.h"
198 #include "AliClusters.h"
199 #include "AliTPCClustersRow.h"
200 #include "AliTPCClustersArray.h"
201
202 #include "AliTPCcluster.h"
203 #include "AliTPCclusterer.h"
204 #include "AliTPCtracker.h"
205
206 #include <TInterpreter.h>
207 #include <TTree.h>
208
209
210
211 ClassImp(AliTPC) 
212
213 //_____________________________________________________________________________
214 AliTPC::AliTPC()
215 {
216   //
217   // Default constructor
218   //
219   fIshunt   = 0;
220   fHits     = 0;
221   fDigits   = 0;
222   fNsectors = 0;
223   //MI changes
224   fDigitsArray = 0;
225   fClustersArray = 0;
226   fDefaults = 0;
227   fTrackHits = 0;  
228   fHitType = 2;  
229   fTPCParam = 0;    
230   fNoiseTable = 0;
231
232 }
233  
234 //_____________________________________________________________________________
235 AliTPC::AliTPC(const char *name, const char *title)
236       : AliDetector(name,title)
237 {
238   //
239   // Standard constructor
240   //
241
242   //
243   // Initialise arrays of hits and digits 
244   fHits     = new TClonesArray("AliTPChit",  176);
245   gAlice->AddHitList(fHits);
246   //MI change  
247   fDigitsArray = 0;
248   fClustersArray= 0;
249   fDefaults = 0;
250   //
251   fTrackHits = new AliTPCTrackHits;  //MI - 13.09.2000
252   fTrackHits->SetHitPrecision(0.002);
253   fTrackHits->SetStepPrecision(0.003);  
254   fTrackHits->SetMaxDistance(100); 
255   fNoiseTable =0;
256
257   fHitType = 2;
258   //
259   // Initialise counters
260   fNsectors = 0;
261
262   //
263   fIshunt     =  0;
264   //
265   // Initialise color attributes
266   SetMarkerColor(kYellow);
267
268   //
269   //  Set TPC parameters
270   //
271
272
273   if (!strcmp(title,"Default")) {       
274     fTPCParam = new AliTPCParamSR;
275   } else {
276     cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
277     fTPCParam=0;
278   }
279
280 }
281
282 //_____________________________________________________________________________
283 AliTPC::~AliTPC()
284 {
285   //
286   // TPC destructor
287   //
288
289   fIshunt   = 0;
290   delete fHits;
291   delete fDigits;
292   delete fTPCParam;
293   delete fTrackHits; //MI 15.09.2000
294   if (fNoiseTable) delete [] fNoiseTable;
295
296 }
297
298 //_____________________________________________________________________________
299 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
300 {
301   //
302   // Add a hit to the list
303   //
304   //  TClonesArray &lhits = *fHits;
305   //  new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
306   if (fHitType&1){
307     TClonesArray &lhits = *fHits;
308     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
309   }
310   if (fHitType&2)
311    AddHit2(track,vol,hits);
312 }
313  
314 //_____________________________________________________________________________
315 void AliTPC::BuildGeometry()
316 {
317
318   //
319   // Build TPC ROOT TNode geometry for the event display
320   //
321   TNode *nNode, *nTop;
322   TTUBS *tubs;
323   Int_t i;
324   const int kColorTPC=19;
325   char name[5], title[25];
326   const Double_t kDegrad=TMath::Pi()/180;
327   const Double_t kRaddeg=180./TMath::Pi();
328
329
330   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
331   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
332
333   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
334   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
335
336   Int_t nLo = fTPCParam->GetNInnerSector()/2;
337   Int_t nHi = fTPCParam->GetNOuterSector()/2;  
338
339   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
340   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
341   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
342   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);  
343
344
345   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
346   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
347
348   Double_t rl,ru;
349   
350
351   //
352   // Get ALICE top node
353   //
354
355   nTop=gAlice->GetGeometry()->GetNode("alice");
356
357   //  inner sectors
358
359   rl = fTPCParam->GetInnerRadiusLow();
360   ru = fTPCParam->GetInnerRadiusUp();
361  
362
363   for(i=0;i<nLo;i++) {
364     sprintf(name,"LS%2.2d",i);
365     name[4]='\0';
366     sprintf(title,"TPC low sector %3d",i);
367     title[24]='\0';
368     
369     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
370                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
371     tubs->SetNumberOfDivisions(1);
372     nTop->cd();
373     nNode = new TNode(name,title,name,0,0,0,"");
374     nNode->SetLineColor(kColorTPC);
375     fNodes->Add(nNode);
376   }
377
378   // Outer sectors
379
380   rl = fTPCParam->GetOuterRadiusLow();
381   ru = fTPCParam->GetOuterRadiusUp();
382
383   for(i=0;i<nHi;i++) {
384     sprintf(name,"US%2.2d",i);
385     name[4]='\0';
386     sprintf(title,"TPC upper sector %d",i);
387     title[24]='\0';
388     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
389                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
390     tubs->SetNumberOfDivisions(1);
391     nTop->cd();
392     nNode = new TNode(name,title,name,0,0,0,"");
393     nNode->SetLineColor(kColorTPC);
394     fNodes->Add(nNode);
395   }
396
397 }    
398
399 //_____________________________________________________________________________
400 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
401 {
402   //
403   // Calculate distance from TPC to mouse on the display
404   // Dummy procedure
405   //
406   return 9999;
407 }
408
409 void AliTPC::Clusters2Tracks(TFile *of) {
410   //-----------------------------------------------------------------
411   // This is a track finder.
412   //-----------------------------------------------------------------
413   AliTPCtracker tracker(fTPCParam);
414   tracker.Clusters2Tracks(gFile,of);
415 }
416
417 //_____________________________________________________________________________
418 void AliTPC::CreateMaterials()
419 {
420   //-----------------------------------------------
421   // Create Materials for for TPC simulations
422   //-----------------------------------------------
423
424   //-----------------------------------------------------------------
425   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
426   //-----------------------------------------------------------------
427
428   Int_t iSXFLD=gAlice->Field()->Integ();
429   Float_t sXMGMX=gAlice->Field()->Max();
430
431   Float_t amat[5]; // atomic numbers
432   Float_t zmat[5]; // z
433   Float_t wmat[5]; // proportions
434
435   Float_t density;
436   Float_t apure[2];
437
438
439   //***************** Gases *************************
440   
441   //-------------------------------------------------
442   // pure gases
443   //-------------------------------------------------
444
445   // Neon
446
447
448   amat[0]= 20.18;
449   zmat[0]= 10.;  
450   density = 0.0009;
451  
452   apure[0]=amat[0];
453
454   AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
455
456   // Argon
457
458   amat[0]= 39.948;
459   zmat[0]= 18.;  
460   density = 0.001782;  
461
462   apure[1]=amat[0];
463
464   AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
465  
466
467   //--------------------------------------------------------------
468   // gases - compounds
469   //--------------------------------------------------------------
470
471   Float_t amol[3];
472
473   // CO2
474
475   amat[0]=12.011;
476   amat[1]=15.9994;
477
478   zmat[0]=6.;
479   zmat[1]=8.;
480
481   wmat[0]=1.;
482   wmat[1]=2.;
483
484   density=0.001977;
485
486   amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
487
488   AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
489   
490   // CF4
491
492   amat[0]=12.011;
493   amat[1]=18.998;
494
495   zmat[0]=6.;
496   zmat[1]=9.;
497  
498   wmat[0]=1.;
499   wmat[1]=4.;
500  
501   density=0.003034;
502
503   amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
504
505   AliMixture(11,"CF4",amat,zmat,density,-2,wmat); 
506
507
508   // CH4
509
510   amat[0]=12.011;
511   amat[1]=1.;
512
513   zmat[0]=6.;
514   zmat[1]=1.;
515
516   wmat[0]=1.;
517   wmat[1]=4.;
518
519   density=0.000717;
520
521   amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
522
523   AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
524
525   //----------------------------------------------------------------
526   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
527   //----------------------------------------------------------------
528
529   char namate[21]; 
530   density = 0.;
531   Float_t am=0;
532   Int_t nc;
533   Float_t rho,absl,X0,buf[1];
534   Int_t nbuf;
535   Float_t a,z;
536
537   for(nc = 0;nc<fNoComp;nc++)
538     {
539     
540       // retrive material constants
541       
542       gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
543
544       amat[nc] = a;
545       zmat[nc] = z;
546
547       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
548  
549       am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]); 
550       density += fMixtProp[nc]*rho;  // density of the mixture
551       
552     }
553
554   // mixture proportions by weight!
555
556   for(nc = 0;nc<fNoComp;nc++)
557     {
558
559       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
560
561       wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? 
562                  apure[nnc] : amol[nnc])/am;
563
564     } 
565
566   // Drift gases 1 - nonsensitive, 2 - sensitive
567
568   AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
569   AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
570
571
572   // Air
573
574   amat[0] = 14.61;
575   zmat[0] = 7.3;
576   density = 0.001205;
577
578   AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.); 
579
580
581   //----------------------------------------------------------------------
582   //               solid materials
583   //----------------------------------------------------------------------
584
585
586   // Kevlar C14H22O2N2
587
588   amat[0] = 12.011;
589   amat[1] = 1.;
590   amat[2] = 15.999;
591   amat[3] = 14.006;
592
593   zmat[0] = 6.;
594   zmat[1] = 1.;
595   zmat[2] = 8.;
596   zmat[3] = 7.;
597
598   wmat[0] = 14.;
599   wmat[1] = 22.;
600   wmat[2] = 2.;
601   wmat[3] = 2.;
602
603   density = 1.45;
604
605   AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);  
606
607   // NOMEX
608
609   amat[0] = 12.011;
610   amat[1] = 1.;
611   amat[2] = 15.999;
612   amat[3] = 14.006;
613
614   zmat[0] = 6.;
615   zmat[1] = 1.;
616   zmat[2] = 8.;
617   zmat[3] = 7.;
618
619   wmat[0] = 14.;
620   wmat[1] = 22.;
621   wmat[2] = 2.;
622   wmat[3] = 2.;
623
624   density = 0.03;
625
626   
627   AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
628
629   // Makrolon C16H18O3
630
631   amat[0] = 12.011;
632   amat[1] = 1.;
633   amat[2] = 15.999;
634
635   zmat[0] = 6.;
636   zmat[1] = 1.;
637   zmat[2] = 8.;
638
639   wmat[0] = 16.;
640   wmat[1] = 18.;
641   wmat[2] = 3.;
642   
643   density = 1.2;
644
645   AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
646   
647   // Mylar C5H4O2
648
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]=5.;
658   wmat[1]=4.;
659   wmat[2]=2.; 
660
661   density = 1.39;
662   
663   AliMixture(37, "Mylar",amat,zmat,density,-3,wmat); 
664
665   // SiO2 - used later for the glass fiber
666
667   amat[0]=28.086;
668   amat[1]=15.9994;
669
670   zmat[0]=14.;
671   zmat[1]=8.;
672
673   wmat[0]=1.;
674   wmat[1]=2.;
675
676
677   AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2)
678
679   // Al
680
681   amat[0] = 26.98;
682   zmat[0] = 13.;
683
684   density = 2.7;
685
686   AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
687
688   // Si
689
690   amat[0] = 28.086;
691   zmat[0] = 14.;
692
693   density = 2.33;
694
695   AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
696
697   // Cu
698
699   amat[0] = 63.546;
700   zmat[0] = 29.;
701
702   density = 8.96;
703
704   AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
705
706   // Tedlar C2H3F
707
708   amat[0] = 12.011;
709   amat[1] = 1.;
710   amat[2] = 18.998;
711
712   zmat[0] = 6.;
713   zmat[1] = 1.;
714   zmat[2] = 9.;
715
716   wmat[0] = 2.;
717   wmat[1] = 3.; 
718   wmat[2] = 1.;
719
720   density = 1.71;
721
722   AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);  
723
724
725   // Plexiglas  C5H8O2
726
727   amat[0]=12.011;
728   amat[1]=1.;
729   amat[2]=15.9994;
730
731   zmat[0]=6.;
732   zmat[1]=1.;
733   zmat[2]=8.;
734
735   wmat[0]=5.;
736   wmat[1]=8.;
737   wmat[2]=2.;
738
739   density=1.18;
740
741   AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
742
743   // Epoxy - C14 H20 O3
744
745   
746   amat[0]=12.011;
747   amat[1]=1.;
748   amat[2]=15.9994;
749
750   zmat[0]=6.;
751   zmat[1]=1.;
752   zmat[2]=8.;
753
754   wmat[0]=14.;
755   wmat[1]=20.;
756   wmat[2]=3.;
757
758   density=1.25;
759
760   AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat);
761
762   // Carbon
763
764   amat[0]=12.011;
765   zmat[0]=6.;
766   density= 2.265;
767
768   AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.);
769
770   // get epoxy
771
772   gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf);
773
774   // Carbon fiber
775
776   wmat[0]=0.644; // by weight!
777   wmat[1]=0.356;
778
779   density=0.5*(1.25+2.265);
780
781   AliMixture(47,"Cfiber",amat,zmat,density,2,wmat);
782
783   // get SiO2
784
785   gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf); 
786
787   wmat[0]=0.725; // by weight!
788   wmat[1]=0.275;
789
790   density=1.7;
791
792   AliMixture(39,"G10",amat,zmat,density,2,wmat);
793
794  
795
796
797   //----------------------------------------------------------
798   // tracking media for gases
799   //----------------------------------------------------------
800
801   AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
802   AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
803   AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
804   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001); 
805
806   //-----------------------------------------------------------  
807   // tracking media for solids
808   //-----------------------------------------------------------
809   
810   AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
811   AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
812   AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
813   AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
814   AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
815   AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
816   AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
817   AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
818   AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
819   AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
820   AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
821   AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
822     
823 }
824
825 void AliTPC::GenerNoise(Int_t tablesize)
826 {
827   //
828   //Generate table with noise
829   //
830   if (fTPCParam==0) {
831     // error message
832     fNoiseDepth=0;
833     return;
834   }
835   if (fNoiseTable)  delete[] fNoiseTable;
836   fNoiseTable = new Float_t[tablesize];
837   fNoiseDepth = tablesize; 
838   fCurrentNoise =0; //!index of the noise in  the noise table 
839   
840   Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
841   for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);      
842 }
843
844 Float_t AliTPC::GetNoise()
845 {
846   // get noise from table
847   //  if ((fCurrentNoise%10)==0) 
848   //  fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
849   if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
850   return fNoiseTable[fCurrentNoise++];
851   //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); 
852 }
853
854
855
856 void AliTPC::Digits2Clusters(TFile *of, Int_t eventnumber)
857 {
858   //-----------------------------------------------------------------
859   // This is a simple cluster finder.
860   //-----------------------------------------------------------------
861   AliTPCclusterer::Digits2Clusters(fTPCParam,of,eventnumber);
862 }
863
864 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
865 extern Double_t SigmaZ2(Double_t, Double_t);
866 //_____________________________________________________________________________
867 void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
868 {
869   //--------------------------------------------------------
870   // TPC simple cluster generator from hits
871   // obtained from the TPC Fast Simulator
872   // The point errors are taken from the parametrization
873   //--------------------------------------------------------
874
875   //-----------------------------------------------------------------
876   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
877   //-----------------------------------------------------------------
878   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
879   // Jouri.Belikov@cern.ch
880   //----------------------------------------------------------------
881   
882   /////////////////////////////////////////////////////////////////////////////
883   //
884   //---------------------------------------------------------------------
885   //   ALICE TPC Cluster Parameters
886   //--------------------------------------------------------------------
887        
888   
889
890   // Cluster width in rphi
891   const Float_t kACrphi=0.18322;
892   const Float_t kBCrphi=0.59551e-3;
893   const Float_t kCCrphi=0.60952e-1;
894   // Cluster width in z
895   const Float_t kACz=0.19081;
896   const Float_t kBCz=0.55938e-3;
897   const Float_t kCCz=0.30428;
898
899   TDirectory *savedir=gDirectory; 
900
901   if (!of->IsOpen()) {
902      cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
903      return;
904   }
905
906   //if(fDefaults == 0) SetDefaults();
907
908   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
909   //
910   TParticle *particle; // pointer to a given particle
911   AliTPChit *tpcHit; // pointer to a sigle TPC hit
912   Int_t sector;
913   Int_t ipart;
914   Float_t xyz[5];
915   Float_t pl,pt,tanth,rpad,ratio;
916   Float_t cph,sph;
917   
918   //---------------------------------------------------------------
919   //  Get the access to the tracks 
920   //---------------------------------------------------------------
921   
922   TTree *tH = gAlice->TreeH();
923   Stat_t ntracks = tH->GetEntries();
924
925   //Switch to the output file
926   of->cd();
927
928   char   cname[100];
929
930   sprintf(cname,"TreeC_TPC_%d",eventn);
931
932   fTPCParam->Write(fTPCParam->GetTitle());
933   AliTPCClustersArray carray;
934   carray.Setup(fTPCParam);
935   carray.SetClusterType("AliTPCcluster");
936   carray.MakeTree();
937
938   Int_t nclusters=0; //cluster counter
939   
940   //------------------------------------------------------------
941   // Loop over all sectors (72 sectors for 20 deg
942   // segmentation for both lower and upper sectors)
943   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
944   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
945   //
946   // First cluster for sector 0 starts at "0"
947   //------------------------------------------------------------
948    
949   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
950     //MI change
951     fTPCParam->AdjustCosSin(isec,cph,sph);
952     
953     //------------------------------------------------------------
954     // Loop over tracks
955     //------------------------------------------------------------
956     
957     for(Int_t track=0;track<ntracks;track++){
958       ResetHits();
959       tH->GetEvent(track);
960       //
961       //  Get number of the TPC hits
962       //
963       // nhits=fHits->GetEntriesFast();
964       //
965      
966        tpcHit = (AliTPChit*)FirstHit(-1);
967
968       // Loop over hits
969       //
970        //   for(Int_t hit=0;hit<nhits;hit++){
971        //tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
972
973        while(tpcHit){
974  
975          if (tpcHit->fQ == 0.) {
976            tpcHit = (AliTPChit*) NextHit();
977            continue; //information about track (I.Belikov)
978          }
979         sector=tpcHit->fSector; // sector number
980
981
982         //      if(sector != isec) continue; //terminate iteration
983
984        if(sector != isec){
985          tpcHit = (AliTPChit*) NextHit();
986          continue; 
987        }
988         ipart=tpcHit->Track();
989         particle=gAlice->Particle(ipart);
990         pl=particle->Pz();
991         pt=particle->Pt();
992         if(pt < 1.e-9) pt=1.e-9;
993         tanth=pl/pt;
994         tanth = TMath::Abs(tanth);
995         rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
996         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
997
998         //   space-point resolutions
999         
1000         sigmaRphi=SigmaY2(rpad,tanth,pt);
1001         sigmaZ   =SigmaZ2(rpad,tanth   );
1002         
1003         //   cluster widths
1004         
1005         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
1006         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
1007         
1008         // temporary protection
1009         
1010         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
1011         if(sigmaZ < 0.) sigmaZ=0.4e-3;
1012         if(clRphi < 0.) clRphi=2.5e-3;
1013         if(clZ < 0.) clZ=2.5e-5;
1014         
1015         //
1016         
1017         //
1018         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
1019         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1020         //
1021         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
1022         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
1023         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
1024           Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
1025           fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
1026           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
1027           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim; 
1028         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
1029           if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z(); 
1030         xyz[2]=sigmaRphi;                                     // fSigmaY2
1031         xyz[3]=sigmaZ;                                        // fSigmaZ2
1032         xyz[4]=tpcHit->fQ;                                    // q
1033
1034         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
1035         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);     
1036
1037         Int_t tracks[3]={tpcHit->Track(), -1, -1};
1038         AliTPCcluster cluster(tracks,xyz);
1039
1040         clrow->InsertCluster(&cluster); nclusters++;
1041
1042         tpcHit = (AliTPChit*)NextHit();
1043         
1044
1045       } // end of loop over hits
1046
1047     }   // end of loop over tracks
1048
1049     Int_t nrows=fTPCParam->GetNRow(isec);
1050     for (Int_t irow=0; irow<nrows; irow++) {
1051         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
1052         if (!clrow) continue;
1053         carray.StoreRow(isec,irow);
1054         carray.ClearRow(isec,irow);
1055     }
1056
1057   } // end of loop over sectors  
1058
1059   cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
1060
1061   carray.GetTree()->SetName(cname);
1062   carray.GetTree()->Write();
1063
1064   savedir->cd(); //switch back to the input file
1065   
1066 } // end of function
1067
1068 //_________________________________________________________________
1069 void AliTPC::Hits2ExactClustersSector(Int_t isec)
1070 {
1071   //--------------------------------------------------------
1072   //calculate exact cross point of track and given pad row
1073   //resulting values are expressed in "digit" coordinata
1074   //--------------------------------------------------------
1075
1076   //-----------------------------------------------------------------
1077   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
1078   //-----------------------------------------------------------------
1079   //
1080   if (fClustersArray==0){    
1081     return;
1082   }
1083   //
1084   TParticle *particle; // pointer to a given particle
1085   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1086   //  Int_t sector,nhits;
1087   Int_t ipart;
1088   const Int_t kcmaxhits=30000;
1089   TVector * xxxx = new TVector(kcmaxhits*4);
1090   TVector & xxx = *xxxx;
1091   Int_t maxhits = kcmaxhits;
1092   //construct array for each padrow
1093   for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++) 
1094     fClustersArray->CreateRow(isec,i);
1095   
1096   //---------------------------------------------------------------
1097   //  Get the access to the tracks 
1098   //---------------------------------------------------------------
1099   
1100   TTree *tH = gAlice->TreeH();
1101   Stat_t ntracks = tH->GetEntries();
1102   Int_t npart = gAlice->GetNtrack();
1103   TBranch * br = tH->GetBranch("fTrackHitsInfo");
1104   //MI change
1105   TBranch * branch=0;
1106   if (fHitType&2) branch = tH->GetBranch("TPC2");
1107   else branch = tH->GetBranch("TPC");
1108
1109   //------------------------------------------------------------
1110   // Loop over tracks
1111   //------------------------------------------------------------
1112   
1113   for(Int_t track=0;track<ntracks;track++){ 
1114     Bool_t isInSector=kTRUE;
1115     ResetHits();
1116         
1117     if (fHitType&2) {
1118       isInSector=kFALSE;    
1119       br->GetEvent(track);
1120       AliObjectArray * ar = fTrackHits->fTrackHitsInfo;
1121       for (UInt_t j=0;j<ar->GetSize();j++){
1122         if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE;
1123       }
1124     }
1125     if (!isInSector) continue;
1126     //MI change
1127     branch->GetEntry(track); // get next track
1128
1129
1130     //
1131     //  Get number of the TPC hits and a pointer
1132     //  to the particles
1133     //
1134     //nhits=fHits->GetEntriesFast();
1135     //
1136     // Loop over hits
1137     //
1138     Int_t currentIndex=0;
1139     Int_t lastrow=-1;  //last writen row
1140
1141     //M.I. changes
1142
1143     tpcHit = (AliTPChit*)FirstHit(-1);
1144     while(tpcHit){
1145       
1146       Int_t sector=tpcHit->fSector; // sector number
1147       if(sector != isec){
1148         tpcHit = (AliTPChit*) NextHit();
1149         continue; 
1150       }
1151
1152       //for(Int_t hit=0;hit<nhits;hit++){
1153       //tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1154       //if (tpcHit==0) continue;
1155       //sector=tpcHit->fSector; // sector number
1156       //if(sector != isec) continue; 
1157       ipart=tpcHit->Track();
1158       if (ipart<npart) particle=gAlice->Particle(ipart);
1159       
1160       //find row number
1161
1162       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
1163       Int_t    index[3]={1,isec,0};
1164       Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;     
1165       if (currentrow<0) {tpcHit = (AliTPChit*)NextHit(); continue;}
1166       if (lastrow<0) lastrow=currentrow;
1167       if (currentrow==lastrow){
1168         if ( currentIndex>=maxhits){
1169           maxhits+=kcmaxhits;
1170           xxx.ResizeTo(4*maxhits);
1171         }     
1172         xxx(currentIndex*4)=x[0];
1173         xxx(currentIndex*4+1)=x[1];
1174         xxx(currentIndex*4+2)=x[2];     
1175         xxx(currentIndex*4+3)=tpcHit->fQ;
1176         currentIndex++; 
1177       }
1178       else 
1179         if (currentIndex>2){
1180           Float_t sumx=0;
1181           Float_t sumx2=0;
1182           Float_t sumx3=0;
1183           Float_t sumx4=0;
1184           Float_t sumy=0;
1185           Float_t sumxy=0;
1186           Float_t sumx2y=0;
1187           Float_t sumz=0;
1188           Float_t sumxz=0;
1189           Float_t sumx2z=0;
1190           Float_t sumq=0;
1191           for (Int_t index=0;index<currentIndex;index++){
1192             Float_t x,x2,x3,x4;
1193             x=x2=x3=x4=xxx(index*4);
1194             x2*=x;
1195             x3*=x2;
1196             x4*=x3;
1197             sumx+=x;
1198             sumx2+=x2;
1199             sumx3+=x3;
1200             sumx4+=x4;
1201             sumy+=xxx(index*4+1);
1202             sumxy+=xxx(index*4+1)*x;
1203             sumx2y+=xxx(index*4+1)*x2;
1204             sumz+=xxx(index*4+2);
1205             sumxz+=xxx(index*4+2)*x;
1206             sumx2z+=xxx(index*4+2)*x2;   
1207             sumq+=xxx(index*4+3);
1208           }
1209           Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
1210           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
1211             sumx2*(sumx*sumx3-sumx2*sumx2);
1212           
1213           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
1214             sumx2*(sumxy*sumx3-sumx2y*sumx2);
1215           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
1216             sumx2*(sumxz*sumx3-sumx2z*sumx2);
1217           
1218           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
1219             sumx2*(sumx*sumx2y-sumx2*sumxy);
1220           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
1221             sumx2*(sumx*sumx2z-sumx2*sumxz);
1222           
1223           Float_t y=detay/det+centralPad;
1224           Float_t z=detaz/det;  
1225           Float_t by=detby/det; //y angle
1226           Float_t bz=detbz/det; //z angle
1227           sumy/=Float_t(currentIndex);
1228           sumz/=Float_t(currentIndex);
1229
1230           AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
1231           if (row!=0) {
1232             AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
1233             //    AliComplexCluster cl;
1234             cl->fX=z;
1235             cl->fY=y;
1236             cl->fQ=sumq;
1237             cl->fSigmaX2=bz;
1238             cl->fSigmaY2=by;
1239             cl->fTracks[0]=ipart;
1240           }
1241           currentIndex=0;
1242           lastrow=currentrow;
1243         } //end of calculating cluster for given row
1244         
1245         
1246       tpcHit = (AliTPChit*)NextHit();
1247     } // end of loop over hits
1248   }   // end of loop over tracks 
1249   //write padrows to tree 
1250   for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1251     fClustersArray->StoreRow(isec,ii);    
1252     fClustersArray->ClearRow(isec,ii);        
1253   }
1254   xxxx->Delete();
1255  
1256 }
1257
1258
1259
1260 //__
1261 void AliTPC::SDigits2Digits2(Int_t eventnumber)  
1262 {
1263   //create digits from summable digits
1264   GenerNoise(500000); //create teble with noise
1265   char  sname[100];
1266   char  dname[100];
1267   sprintf(sname,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
1268   sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1269
1270   //conect tree with sSDigits
1271   TTree *t = (TTree *)gDirectory->Get(sname); 
1272   AliSimDigits digarr, *dummy=&digarr;
1273   t->GetBranch("Segment")->SetAddress(&dummy);
1274   Stat_t nentries = t->GetEntries();
1275
1276   // set zero suppression
1277
1278   fTPCParam->SetZeroSup(2);
1279
1280   // get zero suppression
1281
1282   Int_t zerosup = fTPCParam->GetZeroSup();
1283
1284   //make tree with digits 
1285   
1286   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1287   arr->SetClass("AliSimDigits");
1288   arr->Setup(fTPCParam);
1289   arr->MakeTree(fDigitsFile);
1290   
1291   AliTPCParam * par =fTPCParam;
1292
1293   //Loop over segments of the TPC
1294
1295   for (Int_t n=0; n<nentries; n++) {
1296     t->GetEvent(n);
1297     Int_t sec, row;
1298     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1299       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
1300       continue;
1301     }
1302     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1303     Int_t nrows = digrow->GetNRows();
1304     Int_t ncols = digrow->GetNCols();
1305
1306     digrow->ExpandBuffer();
1307     digarr.ExpandBuffer();
1308     digrow->ExpandTrackBuffer();
1309     digarr.ExpandTrackBuffer();
1310
1311     
1312     Short_t * pamp0 = digarr.GetDigits();
1313     Int_t   * ptracks0 = digarr.GetTracks();
1314     Short_t * pamp1 = digrow->GetDigits();
1315     Int_t   * ptracks1 = digrow->GetTracks();
1316     Int_t  nelems =nrows*ncols;
1317     Int_t saturation = fTPCParam->GetADCSat();
1318     //use internal structure of the AliDigits - for speed reason
1319     //if you cahnge implementation
1320     //of the Alidigits - it must be rewriten -
1321     for (Int_t i= 0; i<nelems; i++){
1322       //      Float_t q = *pamp0;
1323       //q/=16.;  //conversion faktor
1324       //Float_t noise= GetNoise(); 
1325       //q+=noise;      
1326       //q= TMath::Nint(q);
1327       Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
1328       if (q>zerosup){
1329         if (q>saturation) q=saturation;      
1330         *pamp1=(Short_t)q;
1331         //if (ptracks0[0]==0)
1332         //  ptracks1[0]=1;
1333         //else
1334         ptracks1[0]=ptracks0[0];        
1335         ptracks1[nelems]=ptracks0[nelems];
1336         ptracks1[2*nelems]=ptracks0[2*nelems];
1337       }
1338       pamp0++;
1339       pamp1++;
1340       ptracks0++;
1341       ptracks1++;        
1342     }
1343
1344     arr->StoreRow(sec,row);
1345     arr->ClearRow(sec,row);   
1346     // cerr<<sec<<"\t"<<row<<"\n";   
1347   }  
1348
1349     
1350   //write results
1351
1352   
1353   arr->GetTree()->SetName(dname);  
1354   arr->GetTree()->Write();  
1355   delete arr;
1356 }
1357 //_________________________________________
1358 void AliTPC::Merge(TTree * intree, Int_t *mask, Int_t nin, Int_t outid)
1359 {
1360   
1361   //intree - pointer to array of trees with s digits
1362   //mask   - mask for each 
1363   //nin    - number of inputs
1364   //outid  - event  number of the output event
1365
1366   //create digits from summable digits 
1367   //conect tree with sSDigits
1368
1369     
1370   AliSimDigits ** digarr = new AliSimDigits*[nin];
1371   for (Int_t i1=0;i1<nin; i1++){
1372     digarr[i1]=0;
1373     intree[i1].GetBranch("Segment")->SetAddress(&digarr[i1]);
1374   }
1375   Stat_t nentries = intree[0].GetEntries();
1376   
1377   //make tree with digits   
1378   char  dname[100];
1379   sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),outid);
1380   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1381   arr->SetClass("AliSimDigits");
1382   arr->Setup(fTPCParam);
1383   arr->MakeTree(fDigitsFile);  
1384
1385   // set zero suppression
1386
1387   fTPCParam->SetZeroSup(2);
1388
1389   // get zero suppression
1390
1391   Int_t zerosup = fTPCParam->GetZeroSup();
1392
1393
1394   AliTPCParam * par =fTPCParam;
1395
1396   //Loop over segments of the TPC
1397   for (Int_t n=0; n<nentries; n++) {
1398     
1399     for (Int_t i=0;i<nin; i++){ 
1400       //connect proper digits
1401       intree[i].GetEvent(n);      
1402       digarr[i]->ExpandBuffer();
1403       digarr[i]->ExpandTrackBuffer();
1404     }      
1405     Int_t sec, row;
1406     if (!par->AdjustSectorRow(digarr[0]->GetID(),sec,row)) {
1407       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
1408       continue;
1409     }
1410
1411     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1412     Int_t nrows = digrow->GetNRows();
1413     Int_t ncols = digrow->GetNCols();
1414
1415     digrow->ExpandBuffer();
1416     digrow->ExpandTrackBuffer();
1417    
1418     for (Int_t rows=0;rows<nrows; rows++){
1419       for (Int_t col=0;col<ncols; col++){
1420         Float_t q=0;
1421         Int_t label[1000]; // i hope no more than 300 events merged
1422         Int_t labptr = 0;
1423         // looop over digits
1424         for (Int_t i=0;i<nin; i++){ 
1425           q  += digarr[i]->GetDigitFast(rows,col);
1426           for (Int_t tr=0;tr<3;tr++) {
1427             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
1428             if ( lab > 1) {
1429               label[labptr]=lab+mask[i]-2;
1430               labptr++;
1431             }
1432           }
1433         }
1434         //add noise
1435         q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()*16); 
1436         
1437         q/=16;  //conversion faktor
1438         q=(Short_t)q;
1439
1440         if (q> zerosup){
1441
1442           if(q > fTPCParam->GetADCSat()) q = (Short_t)(fTPCParam->GetADCSat());
1443           digrow->SetDigitFast(q,rows,col);  
1444           for (Int_t tr=0;tr<3;tr++){
1445             if (tr<labptr)
1446               ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],rows,col,tr);
1447             else
1448               ((AliSimDigits*)digrow)->SetTrackIDFast(-1,rows,col,tr);
1449           }
1450         }
1451       } 
1452      }         
1453       arr->StoreRow(sec,row);
1454
1455       arr->ClearRow(sec,row);   
1456       // cerr<<sec<<"\t"<<row<<"\n";  
1457   }  
1458   
1459   delete digarr;
1460   arr->GetTree()->SetName(dname); 
1461   arr->GetTree()->Write(); 
1462  
1463   delete arr;  
1464   
1465 }
1466
1467 /*_________________________________________
1468 void AliTPC::SDigits2Digits(Int_t eventnumber)
1469 {
1470
1471
1472   cerr<<"Digitizing TPC...\n";
1473
1474   Hits2Digits(eventnumber);
1475    
1476     
1477   //write results
1478
1479   //  char treeName[100];
1480
1481   //  sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1482   
1483   //  GetDigitsArray()->GetTree()->Write(treeName);  
1484 }
1485 */
1486 //__________________________________________________________________
1487 void AliTPC::SetDefaults(){
1488
1489    
1490    cerr<<"Setting default parameters...\n";
1491
1492   // Set response functions
1493
1494   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1495   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1496   AliTPCPRF2D    * prfouter   = new AliTPCPRF2D;
1497   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1498   rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1499   rf->SetOffset(3*param->GetZSigma());
1500   rf->Update();
1501   
1502   TDirectory *savedir=gDirectory;
1503   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1504   if (!f->IsOpen()) { 
1505     cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ;
1506      exit(3);
1507   }
1508   prfinner->Read("prf_07504_Gati_056068_d02");
1509   prfouter->Read("prf_10006_Gati_047051_d03");
1510   f->Close();
1511   savedir->cd();
1512
1513   param->SetInnerPRF(prfinner);
1514   param->SetOuterPRF(prfouter); 
1515   param->SetTimeRF(rf);
1516
1517   // set fTPCParam
1518
1519   SetParam(param);
1520
1521
1522   fDefaults = 1;
1523
1524 }
1525 //__________________________________________________________________  
1526 void AliTPC::Hits2Digits(Int_t eventnumber)  
1527
1528  //----------------------------------------------------
1529  // Loop over all sectors for a single event
1530  //----------------------------------------------------
1531
1532
1533   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1534   GenerNoise(500000); //create teble with noise
1535
1536   //setup TPCDigitsArray 
1537
1538   if(GetDigitsArray()) delete GetDigitsArray();
1539
1540   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1541   arr->SetClass("AliSimDigits");
1542   arr->Setup(fTPCParam);
1543   arr->MakeTree(fDigitsFile);
1544   SetDigitsArray(arr);
1545
1546   fDigitsSwitch=0; // standard digits
1547
1548   cerr<<"Digitizing TPC...\n";
1549
1550  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1551
1552   // write results
1553
1554   char treeName[100];
1555
1556   sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1557   
1558   GetDigitsArray()->GetTree()->SetName(treeName);  
1559   GetDigitsArray()->GetTree()->Write();  
1560
1561
1562 }
1563
1564
1565
1566 //__________________________________________________________________
1567 void AliTPC::Hits2SDigits2(Int_t eventnumber)  
1568
1569
1570   //-----------------------------------------------------------
1571   //   summable digits - 16 bit "ADC", no noise, no saturation
1572   //-----------------------------------------------------------
1573
1574  //----------------------------------------------------
1575  // Loop over all sectors for a single event
1576  //----------------------------------------------------
1577
1578
1579   if(fDefaults == 0) SetDefaults();
1580   GenerNoise(500000); //create table with noise
1581   //setup TPCDigitsArray 
1582
1583   if(GetDigitsArray()) delete GetDigitsArray();
1584
1585   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1586   arr->SetClass("AliSimDigits");
1587   arr->Setup(fTPCParam);
1588   arr->MakeTree(fDigitsFile);
1589   SetDigitsArray(arr);
1590
1591   cerr<<"Digitizing TPC...\n"; 
1592
1593   fDigitsSwitch=1; // summable digits
1594   
1595     // set zero suppression to "0"
1596
1597   fTPCParam->SetZeroSup(0);
1598
1599  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1600
1601
1602   // write results
1603
1604   char treeName[100];
1605
1606   sprintf(treeName,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
1607   
1608   GetDigitsArray()->GetTree()->SetName(treeName); 
1609   GetDigitsArray()->GetTree()->Write(); 
1610
1611 }
1612
1613
1614
1615 //__________________________________________________________________
1616 void AliTPC::Hits2SDigits()  
1617
1618
1619   //-----------------------------------------------------------
1620   //   summable digits - 16 bit "ADC", no noise, no saturation
1621   //-----------------------------------------------------------
1622
1623  //----------------------------------------------------
1624  // Loop over all sectors for a single event
1625  //----------------------------------------------------
1626   //MI change - for pp run
1627   Int_t eventnumber = gAlice->GetEvNumber();
1628
1629   if(fDefaults == 0) SetDefaults();
1630   GenerNoise(500000); //create table with noise
1631
1632   //setup TPCDigitsArray 
1633
1634   if(GetDigitsArray()) delete GetDigitsArray();
1635
1636   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
1637   arr->SetClass("AliSimDigits");
1638   arr->Setup(fTPCParam);
1639   arr->MakeTree(fDigitsFile);
1640   SetDigitsArray(arr);
1641
1642   cerr<<"Digitizing TPC...\n"; 
1643
1644   //  fDigitsSwitch=1; // summable digits  -for the moment direct
1645
1646  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1647
1648
1649   // write results
1650   char treeName[100];
1651
1652   sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1653   
1654   GetDigitsArray()->GetTree()->SetName(treeName); 
1655   GetDigitsArray()->GetTree()->Write(); 
1656
1657 }
1658
1659
1660 //_____________________________________________________________________________
1661 void AliTPC::Hits2DigitsSector(Int_t isec)
1662 {
1663   //-------------------------------------------------------------------
1664   // TPC conversion from hits to digits.
1665   //------------------------------------------------------------------- 
1666
1667   //-----------------------------------------------------------------
1668   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1669   //-----------------------------------------------------------------
1670
1671   //-------------------------------------------------------
1672   //  Get the access to the track hits
1673   //-------------------------------------------------------
1674
1675   // check if the parameters are set - important if one calls this method
1676   // directly, not from the Hits2Digits
1677
1678   if(fDefaults == 0) SetDefaults();
1679
1680   TTree *tH = gAlice->TreeH(); // pointer to the hits tree
1681   Stat_t ntracks = tH->GetEntries();
1682
1683   if( ntracks > 0){
1684
1685   //------------------------------------------- 
1686   //  Only if there are any tracks...
1687   //-------------------------------------------
1688
1689     TObjArray **row;
1690     
1691     //printf("*** Processing sector number %d ***\n",isec);
1692
1693       Int_t nrows =fTPCParam->GetNRow(isec);
1694
1695       row= new TObjArray* [nrows];
1696     
1697       MakeSector(isec,nrows,tH,ntracks,row);
1698
1699       //--------------------------------------------------------
1700       //   Digitize this sector, row by row
1701       //   row[i] is the pointer to the TObjArray of TVectors,
1702       //   each one containing electrons accepted on this
1703       //   row, assigned into tracks
1704       //--------------------------------------------------------
1705
1706       Int_t i;
1707
1708       if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree(fDigitsFile);
1709
1710       for (i=0;i<nrows;i++){
1711
1712         AliDigits * dig = fDigitsArray->CreateRow(isec,i); 
1713
1714         DigitizeRow(i,isec,row);
1715
1716         fDigitsArray->StoreRow(isec,i);
1717
1718                 Int_t ndig = dig->GetDigitSize(); 
1719         
1720  
1721        //printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1722         
1723         fDigitsArray->ClearRow(isec,i);  
1724
1725    
1726        } // end of the sector digitization
1727
1728       for(i=0;i<nrows;i++){
1729         row[i]->Delete();  
1730         delete row[i];   
1731       }
1732       
1733        delete [] row; // delete the array of pointers to TObjArray-s
1734         
1735   } // ntracks >0
1736
1737 } // end of Hits2DigitsSector
1738
1739
1740 //_____________________________________________________________________________
1741 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1742 {
1743   //-----------------------------------------------------------
1744   // Single row digitization, coupling from the neighbouring
1745   // rows taken into account
1746   //-----------------------------------------------------------
1747
1748   //-----------------------------------------------------------------
1749   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1750   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1751   //-----------------------------------------------------------------
1752  
1753
1754   Float_t zerosup = fTPCParam->GetZeroSup();
1755   Int_t nrows =fTPCParam->GetNRow(isec);
1756   fCurrentIndex[1]= isec;
1757   
1758
1759   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1760   Int_t nofTbins = fTPCParam->GetMaxTBin();
1761   Int_t indexRange[4];
1762   //
1763   //  Integrated signal for this row
1764   //  and a single track signal
1765   //    
1766   TMatrix *m1   = new TMatrix(0,nofPads,0,nofTbins); // integrated
1767   TMatrix *m2   = new TMatrix(0,nofPads,0,nofTbins); // single
1768   //
1769   TMatrix &total  = *m1;
1770
1771   //  Array of pointers to the label-signal list
1772
1773   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1774   Float_t  **pList = new Float_t* [nofDigits]; 
1775
1776   Int_t lp;
1777   Int_t i1;   
1778   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1779   //
1780   //calculate signal 
1781   //
1782   Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1783   Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1784   for (Int_t row= row1;row<=row2;row++){
1785     Int_t nTracks= rows[row]->GetEntries();
1786     for (i1=0;i1<nTracks;i1++){
1787       fCurrentIndex[2]= row;
1788       fCurrentIndex[3]=irow;
1789       if (row==irow){
1790         m2->Zero();  // clear single track signal matrix
1791         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange); 
1792         GetList(trackLabel,nofPads,m2,indexRange,pList);
1793       }
1794       else   GetSignal(rows[row],i1,0,m1,indexRange);
1795     }
1796   }
1797          
1798   Int_t tracks[3];
1799
1800   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1801   for(Int_t ip=0;ip<nofPads;ip++){
1802     for(Int_t it=0;it<nofTbins;it++){
1803
1804       Float_t q = total(ip,it);
1805
1806       Int_t gi =it*nofPads+ip; // global index
1807
1808       if(fDigitsSwitch == 0){
1809
1810         //        q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); 
1811         q+=GetNoise();
1812         q = TMath::Nint(q);
1813
1814         if(q <=zerosup) continue; // do not fill zeros
1815         if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();  // saturation
1816
1817       }
1818
1819       else {
1820        if(q <= 0.) continue; // do not fill zeros
1821        if(q>2000.) q=2000.;
1822        q *= 16.;
1823        q = TMath::Nint(q);
1824       }
1825
1826       //
1827       //  "real" signal or electronic noise (list = -1)?
1828       //    
1829
1830       for(Int_t j1=0;j1<3;j1++){
1831         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1832       }
1833
1834 //Begin_Html
1835 /*
1836   <A NAME="AliDigits"></A>
1837   using of AliDigits object
1838 */
1839 //End_Html
1840       dig->SetDigitFast((Short_t)q,it,ip);
1841       if (fDigitsArray->IsSimulated())
1842         {
1843          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1844          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1845          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1846         }
1847      
1848     
1849     } // end of loop over time buckets
1850   }  // end of lop over pads 
1851
1852   //
1853   //  This row has been digitized, delete nonused stuff
1854   //
1855
1856   for(lp=0;lp<nofDigits;lp++){
1857     if(pList[lp]) delete [] pList[lp];
1858   }
1859   
1860   delete [] pList;
1861
1862   delete m1;
1863   delete m2;
1864   //  delete m3;
1865
1866 } // end of DigitizeRow
1867
1868 //_____________________________________________________________________________
1869
1870 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
1871                           Int_t *indexRange)
1872 {
1873
1874   //---------------------------------------------------------------
1875   //  Calculates 2-D signal (pad,time) for a single track,
1876   //  returns a pointer to the signal matrix and the track label 
1877   //  No digitization is performed at this level!!!
1878   //---------------------------------------------------------------
1879
1880   //-----------------------------------------------------------------
1881   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1882   // Modified: Marian Ivanov 
1883   //-----------------------------------------------------------------
1884
1885   TVector *tv;
1886  
1887   tv = (TVector*)p1->At(ntr); // pointer to a track
1888   TVector &v = *tv;
1889   
1890   Float_t label = v(0);
1891   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
1892
1893   Int_t nElectrons = (tv->GetNrows()-1)/4;
1894   indexRange[0]=9999; // min pad
1895   indexRange[1]=-1; // max pad
1896   indexRange[2]=9999; //min time
1897   indexRange[3]=-1; // max time
1898
1899   //  Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1900
1901   TMatrix &signal = *m1;
1902   TMatrix &total = *m2;
1903   //
1904   //  Loop over all electrons
1905   //
1906   for(Int_t nel=0; nel<nElectrons; nel++){
1907     Int_t idx=nel*4;
1908     Float_t aval =  v(idx+4);
1909     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac(); 
1910     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1911     Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
1912     
1913     if (n>0) for (Int_t i =0; i<n; i++){
1914        Int_t *index = fTPCParam->GetResBin(i);        
1915        Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1916        if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
1917          Int_t time=index[2];    
1918          Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
1919          weight *= eltoadcfac;
1920          
1921          if (m1!=0) signal(pad,time)+=weight; 
1922          total(pad,time)+=weight;
1923          indexRange[0]=TMath::Min(indexRange[0],pad);
1924          indexRange[1]=TMath::Max(indexRange[1],pad);
1925          indexRange[2]=TMath::Min(indexRange[2],time);
1926          indexRange[3]=TMath::Max(indexRange[3],time); 
1927        }         
1928     }
1929   } // end of loop over electrons
1930   
1931   return label; // returns track label when finished
1932 }
1933
1934 //_____________________________________________________________________________
1935 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
1936                      Float_t **pList)
1937 {
1938   //----------------------------------------------------------------------
1939   //  Updates the list of tracks contributing to digits for a given row
1940   //----------------------------------------------------------------------
1941
1942   //-----------------------------------------------------------------
1943   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1944   //-----------------------------------------------------------------
1945
1946   TMatrix &signal = *m;
1947
1948   // lop over nonzero digits
1949
1950   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1951     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1952
1953
1954         // accept only the contribution larger than 500 electrons (1/2 s_noise)
1955
1956         if(signal(ip,it)<0.5) continue; 
1957
1958
1959         Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1960         
1961         if(!pList[globalIndex]){
1962         
1963           // 
1964           // Create new list (6 elements - 3 signals and 3 labels),
1965           //
1966
1967           pList[globalIndex] = new Float_t [6];
1968
1969           // set list to -1 
1970
1971           *pList[globalIndex] = -1.;
1972           *(pList[globalIndex]+1) = -1.;
1973           *(pList[globalIndex]+2) = -1.;
1974           *(pList[globalIndex]+3) = -1.;
1975           *(pList[globalIndex]+4) = -1.;
1976           *(pList[globalIndex]+5) = -1.;
1977
1978
1979           *pList[globalIndex] = label;
1980           *(pList[globalIndex]+3) = signal(ip,it);
1981         }
1982         else{
1983
1984           // check the signal magnitude
1985
1986           Float_t highest = *(pList[globalIndex]+3);
1987           Float_t middle = *(pList[globalIndex]+4);
1988           Float_t lowest = *(pList[globalIndex]+5);
1989
1990           //
1991           //  compare the new signal with already existing list
1992           //
1993
1994           if(signal(ip,it)<lowest) continue; // neglect this track
1995
1996           //
1997
1998           if (signal(ip,it)>highest){
1999             *(pList[globalIndex]+5) = middle;
2000             *(pList[globalIndex]+4) = highest;
2001             *(pList[globalIndex]+3) = signal(ip,it);
2002
2003             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
2004             *(pList[globalIndex]+1) = *pList[globalIndex];
2005             *pList[globalIndex] = label;
2006           }
2007           else if (signal(ip,it)>middle){
2008             *(pList[globalIndex]+5) = middle;
2009             *(pList[globalIndex]+4) = signal(ip,it);
2010
2011             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
2012             *(pList[globalIndex]+1) = label;
2013           }
2014           else{
2015             *(pList[globalIndex]+5) = signal(ip,it);
2016             *(pList[globalIndex]+2) = label;
2017           }
2018         }
2019
2020     } // end of loop over pads
2021   } // end of loop over time bins
2022
2023
2024
2025 }//end of GetList
2026 //___________________________________________________________________
2027 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
2028                         Stat_t ntracks,TObjArray **row)
2029 {
2030
2031   //-----------------------------------------------------------------
2032   // Prepares the sector digitization, creates the vectors of
2033   // tracks for each row of this sector. The track vector
2034   // contains the track label and the position of electrons.
2035   //-----------------------------------------------------------------
2036
2037   //-----------------------------------------------------------------
2038   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2039   //-----------------------------------------------------------------
2040
2041   Float_t gasgain = fTPCParam->GetGasGain();
2042   Int_t i;
2043   Float_t xyz[4]; 
2044
2045   AliTPChit *tpcHit; // pointer to a sigle TPC hit    
2046   //MI change
2047   TBranch * branch=0;
2048   if (fHitType&2) branch = TH->GetBranch("TPC2");
2049   else branch = TH->GetBranch("TPC");
2050
2051  
2052   //----------------------------------------------
2053   // Create TObjArray-s, one for each row,
2054   // each TObjArray will store the TVectors
2055   // of electrons, one TVector per each track.
2056   //---------------------------------------------- 
2057     
2058   Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
2059   TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
2060   for(i=0; i<nrows; i++){
2061     row[i] = new TObjArray;
2062     nofElectrons[i]=0;
2063     tracks[i]=0;
2064   }
2065
2066  
2067
2068   //--------------------------------------------------------------------
2069   //  Loop over tracks, the "track" contains the full history
2070   //--------------------------------------------------------------------
2071
2072   Int_t previousTrack,currentTrack;
2073   previousTrack = -1; // nothing to store so far!
2074
2075   for(Int_t track=0;track<ntracks;track++){
2076     Bool_t isInSector=kTRUE;
2077     ResetHits();
2078
2079     if (fHitType&2) {
2080       isInSector=kFALSE;
2081       TBranch * br = TH->GetBranch("fTrackHitsInfo");
2082       br->GetEvent(track);
2083       AliObjectArray * ar = fTrackHits->fTrackHitsInfo;
2084       for (UInt_t j=0;j<ar->GetSize();j++){
2085         if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE;
2086       }
2087     }
2088     if (!isInSector) continue;
2089     //MI change
2090     branch->GetEntry(track); // get next track
2091
2092     //M.I. changes
2093
2094     tpcHit = (AliTPChit*)FirstHit(-1);
2095
2096     //--------------------------------------------------------------
2097     //  Loop over hits
2098     //--------------------------------------------------------------
2099
2100
2101     while(tpcHit){
2102       
2103       Int_t sector=tpcHit->fSector; // sector number
2104       //      if(sector != isec) continue; 
2105       if(sector != isec){
2106         tpcHit = (AliTPChit*) NextHit();
2107         continue; 
2108       }
2109
2110         currentTrack = tpcHit->Track(); // track number
2111
2112
2113         if(currentTrack != previousTrack){
2114                           
2115            // store already filled fTrack
2116               
2117            for(i=0;i<nrows;i++){
2118              if(previousTrack != -1){
2119                if(nofElectrons[i]>0){
2120                  TVector &v = *tracks[i];
2121                  v(0) = previousTrack;
2122                  tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2123                  row[i]->Add(tracks[i]);                     
2124                }
2125                else{
2126                  delete tracks[i]; // delete empty TVector
2127                  tracks[i]=0;
2128                }
2129              }
2130
2131              nofElectrons[i]=0;
2132              tracks[i] = new TVector(481); // TVectors for the next fTrack
2133
2134            } // end of loop over rows
2135                
2136            previousTrack=currentTrack; // update track label 
2137         }
2138            
2139         Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2140
2141        //---------------------------------------------------
2142        //  Calculate the electron attachment probability
2143        //---------------------------------------------------
2144
2145
2146         Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
2147                                                         /fTPCParam->GetDriftV(); 
2148         // in microseconds!     
2149         Float_t attProb = fTPCParam->GetAttCoef()*
2150           fTPCParam->GetOxyCont()*time; //  fraction! 
2151    
2152         //-----------------------------------------------
2153         //  Loop over electrons
2154         //-----------------------------------------------
2155         Int_t index[3];
2156         index[1]=isec;
2157         for(Int_t nel=0;nel<qI;nel++){
2158           // skip if electron lost due to the attachment
2159           if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
2160           xyz[0]=tpcHit->X();
2161           xyz[1]=tpcHit->Y();
2162           xyz[2]=tpcHit->Z();   
2163           //
2164           // protection for the nonphysical avalanche size (10**6 maximum)
2165           //  
2166           Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2167           xyz[3]= (Float_t) (-gasgain*TMath::Log(rn)); 
2168           index[0]=1;
2169           
2170           TransportElectron(xyz,index); //MI change -august       
2171           Int_t rowNumber;
2172           fTPCParam->GetPadRow(xyz,index); //MI change august
2173           rowNumber = index[2];
2174           //transform position to local digit coordinates
2175           //relative to nearest pad row 
2176           if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue;       
2177           nofElectrons[rowNumber]++;      
2178           //----------------------------------
2179           // Expand vector if necessary
2180           //----------------------------------
2181           if(nofElectrons[rowNumber]>120){
2182             Int_t range = tracks[rowNumber]->GetNrows();
2183             if((nofElectrons[rowNumber])>(range-1)/4){
2184         
2185               tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
2186             }
2187           }
2188           
2189           TVector &v = *tracks[rowNumber];
2190           Int_t idx = 4*nofElectrons[rowNumber]-3;
2191
2192           v(idx)=  xyz[0];   // X - pad row coordinate
2193           v(idx+1)=xyz[1];   // Y - pad coordinate (along the pad-row)
2194           v(idx+2)=xyz[2];   // Z - time bin coordinate
2195           v(idx+3)=xyz[3];   // avalanche size  
2196         } // end of loop over electrons
2197
2198         tpcHit = (AliTPChit*)NextHit();
2199         
2200       } // end of loop over hits
2201     } // end of loop over tracks
2202
2203     //
2204     //   store remaining track (the last one) if not empty
2205     //
2206
2207      for(i=0;i<nrows;i++){
2208        if(nofElectrons[i]>0){
2209           TVector &v = *tracks[i];
2210           v(0) = previousTrack;
2211           tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2212           row[i]->Add(tracks[i]);  
2213         }
2214         else{
2215           delete tracks[i];
2216           tracks[i]=0;
2217         }  
2218       }  
2219
2220           delete [] tracks;
2221           delete [] nofElectrons;
2222  
2223
2224 } // end of MakeSector
2225
2226
2227 //_____________________________________________________________________________
2228 void AliTPC::Init()
2229 {
2230   //
2231   // Initialise TPC detector after definition of geometry
2232   //
2233   Int_t i;
2234   //
2235   if(fDebug) {
2236     printf("\n%s: ",ClassName());
2237     for(i=0;i<35;i++) printf("*");
2238     printf(" TPC_INIT ");
2239     for(i=0;i<35;i++) printf("*");
2240     printf("\n%s: ",ClassName());
2241     //
2242     for(i=0;i<80;i++) printf("*");
2243     printf("\n");
2244   }
2245 }
2246
2247 //_____________________________________________________________________________
2248 void AliTPC::MakeBranch(Option_t* option, const char *file)
2249 {
2250   //
2251   // Create Tree branches for the TPC.
2252   //
2253   Int_t buffersize = 4000;
2254   char branchname[10];
2255   sprintf(branchname,"%s",GetName());
2256
2257   AliDetector::MakeBranch(option,file);
2258
2259   const char *d = strstr(option,"D");
2260
2261   if (fDigits   && gAlice->TreeD() && d) {
2262       MakeBranchInTree(gAlice->TreeD(), 
2263                        branchname, &fDigits, buffersize, file);
2264   }     
2265
2266   if (fHitType&2) MakeBranch2(option,file); // MI change 14.09.2000
2267 }
2268  
2269 //_____________________________________________________________________________
2270 void AliTPC::ResetDigits()
2271 {
2272   //
2273   // Reset number of digits and the digits array for this detector
2274   //
2275   fNdigits   = 0;
2276   if (fDigits)   fDigits->Clear();
2277 }
2278
2279 //_____________________________________________________________________________
2280 void AliTPC::SetSecAL(Int_t sec)
2281 {
2282   //---------------------------------------------------
2283   // Activate/deactivate selection for lower sectors
2284   //---------------------------------------------------
2285
2286   //-----------------------------------------------------------------
2287   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2288   //-----------------------------------------------------------------
2289
2290   fSecAL = sec;
2291 }
2292
2293 //_____________________________________________________________________________
2294 void AliTPC::SetSecAU(Int_t sec)
2295 {
2296   //----------------------------------------------------
2297   // Activate/deactivate selection for upper sectors
2298   //---------------------------------------------------
2299
2300   //-----------------------------------------------------------------
2301   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2302   //-----------------------------------------------------------------
2303
2304   fSecAU = sec;
2305 }
2306
2307 //_____________________________________________________________________________
2308 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2309 {
2310   //----------------------------------------
2311   // Select active lower sectors
2312   //----------------------------------------
2313
2314   //-----------------------------------------------------------------
2315   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2316   //-----------------------------------------------------------------
2317
2318   fSecLows[0] = s1;
2319   fSecLows[1] = s2;
2320   fSecLows[2] = s3;
2321   fSecLows[3] = s4;
2322   fSecLows[4] = s5;
2323   fSecLows[5] = s6;
2324 }
2325
2326 //_____________________________________________________________________________
2327 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2328                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10, 
2329                        Int_t s11 , Int_t s12)
2330 {
2331   //--------------------------------
2332   // Select active upper sectors
2333   //--------------------------------
2334
2335   //-----------------------------------------------------------------
2336   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2337   //-----------------------------------------------------------------
2338
2339   fSecUps[0] = s1;
2340   fSecUps[1] = s2;
2341   fSecUps[2] = s3;
2342   fSecUps[3] = s4;
2343   fSecUps[4] = s5;
2344   fSecUps[5] = s6;
2345   fSecUps[6] = s7;
2346   fSecUps[7] = s8;
2347   fSecUps[8] = s9;
2348   fSecUps[9] = s10;
2349   fSecUps[10] = s11;
2350   fSecUps[11] = s12;
2351 }
2352
2353 //_____________________________________________________________________________
2354 void AliTPC::SetSens(Int_t sens)
2355 {
2356
2357   //-------------------------------------------------------------
2358   // Activates/deactivates the sensitive strips at the center of
2359   // the pad row -- this is for the space-point resolution calculations
2360   //-------------------------------------------------------------
2361
2362   //-----------------------------------------------------------------
2363   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2364   //-----------------------------------------------------------------
2365
2366   fSens = sens;
2367 }
2368
2369  
2370 void AliTPC::SetSide(Float_t side=0.)
2371 {
2372   // choice of the TPC side
2373
2374   fSide = side;
2375  
2376 }
2377 //____________________________________________________________________________
2378 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2379                            Float_t p2,Float_t p3)
2380 {
2381
2382   // gax mixture definition
2383
2384  fNoComp = nc;
2385  
2386  fMixtComp[0]=c1;
2387  fMixtComp[1]=c2;
2388  fMixtComp[2]=c3;
2389
2390  fMixtProp[0]=p1;
2391  fMixtProp[1]=p2;
2392  fMixtProp[2]=p3; 
2393  
2394  
2395 }
2396 //_____________________________________________________________________________
2397
2398 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2399 {
2400   //
2401   // electron transport taking into account:
2402   // 1. diffusion, 
2403   // 2.ExB at the wires
2404   // 3. nonisochronity
2405   //
2406   // xyz and index must be already transformed to system 1
2407   //
2408
2409   fTPCParam->Transform1to2(xyz,index);
2410   
2411   //add diffusion
2412   Float_t driftl=xyz[2];
2413   if(driftl<0.01) driftl=0.01;
2414   driftl=TMath::Sqrt(driftl);
2415   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2416   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2417   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2418   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2419   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2420
2421   // ExB
2422   
2423   if (fTPCParam->GetMWPCReadout()==kTRUE){
2424     Float_t x1=xyz[0];
2425     fTPCParam->Transform2to2NearestWire(xyz,index);
2426     Float_t dx=xyz[0]-x1;
2427     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2428   }
2429   //add nonisochronity (not implemented yet)
2430   
2431 }
2432   
2433 ClassImp(AliTPCdigit)
2434  
2435 //_____________________________________________________________________________
2436 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2437   AliDigit(tracks)
2438 {
2439   //
2440   // Creates a TPC digit object
2441   //
2442   fSector     = digits[0];
2443   fPadRow     = digits[1];
2444   fPad        = digits[2];
2445   fTime       = digits[3];
2446   fSignal     = digits[4];
2447 }
2448
2449  
2450 ClassImp(AliTPChit)
2451  
2452 //_____________________________________________________________________________
2453 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2454 AliHit(shunt,track)
2455 {
2456   //
2457   // Creates a TPC hit object
2458   //
2459   fSector     = vol[0];
2460   fPadRow     = vol[1];
2461   fX          = hits[0];
2462   fY          = hits[1];
2463   fZ          = hits[2];
2464   fQ          = hits[3];
2465 }
2466  
2467
2468 //________________________________________________________________________
2469 // Additional code because of the AliTPCTrackHits
2470
2471 void AliTPC::MakeBranch2(Option_t *option,const char *file)
2472 {
2473   //
2474   // Create a new branch in the current Root Tree
2475   // The branch of fHits is automatically split
2476   // MI change 14.09.2000
2477   if (fHitType&2==0) return;
2478   char branchname[10];
2479   sprintf(branchname,"%s2",GetName());  
2480   //
2481   // Get the pointer to the header
2482   const char *cH = strstr(option,"H");
2483   //
2484   if (fTrackHits   && gAlice->TreeH() && cH) {    
2485     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHits, 
2486                                                    gAlice->TreeH(),fBufferSize,99);
2487     gAlice->TreeH()->GetListOfBranches()->Add(branch);
2488     if (GetDebug()>1) 
2489       printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname);
2490     const char folder [] = "RunMC/Event/Data";
2491     if (GetDebug())
2492       printf("%15s: Publishing %s to %s\n",ClassName(),branchname,folder);
2493     Publish(folder,&fTrackHits,branchname);
2494     if (file) {
2495         TBranch *b = gAlice->TreeH()->GetBranch(branchname);
2496         TDirectory *wd = gDirectory;
2497         b->SetFile(file);
2498         TIter next( b->GetListOfBranches());
2499         while ((b=(TBranch*)next())) {
2500           b->SetFile(file);
2501         }
2502         wd->cd(); 
2503         if (GetDebug()>1) 
2504               cout << "Diverting branch " << branchname << " to file " << file << endl;  
2505     }
2506   }     
2507 }
2508
2509 void AliTPC::SetTreeAddress()
2510 {
2511   if (fHitType&1) AliDetector::SetTreeAddress();
2512   if (fHitType&2) SetTreeAddress2();
2513 }
2514
2515 void AliTPC::SetTreeAddress2()
2516 {
2517   //
2518   // Set branch address for the TrackHits Tree
2519   // 
2520   TBranch *branch;
2521   char branchname[20];
2522   sprintf(branchname,"%s2",GetName());
2523   //
2524   // Branch address for hit tree
2525   TTree *treeH = gAlice->TreeH();
2526   if (treeH) {
2527     branch = treeH->GetBranch(branchname);
2528     if (branch) branch->SetAddress(&fTrackHits);
2529   }
2530 }
2531
2532 void AliTPC::FinishPrimary()
2533 {
2534   if (fTrackHits) fTrackHits->FlushHitStack();  
2535 }
2536
2537
2538 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2539
2540   //
2541   // add hit to the list  
2542   Int_t rtrack;
2543   if (fIshunt) {
2544     int primary = gAlice->GetPrimary(track);
2545     gAlice->Particle(primary)->SetBit(kKeepBit);
2546     rtrack=primary;
2547   } else {
2548     rtrack=track;
2549     gAlice->FlagTrack(track);
2550   }  
2551   //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
2552   //if (hit->fTrack!=rtrack)
2553   //  cout<<"bad track number\n";
2554   if (fTrackHits) 
2555     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2556                              hits[1],hits[2],(Int_t)hits[3]);
2557 }
2558
2559 void AliTPC::ResetHits()
2560 {
2561   if (fHitType&1) AliDetector::ResetHits();
2562   if (fHitType&2) ResetHits2();
2563 }
2564
2565 void AliTPC::ResetHits2()
2566 {
2567   //
2568   //reset hits
2569   if (fTrackHits) fTrackHits->Clear();
2570 }   
2571
2572 AliHit* AliTPC::FirstHit(Int_t track)
2573 {
2574   if (fHitType&2) return FirstHit2(track);
2575   return AliDetector::FirstHit(track);
2576 }
2577 AliHit* AliTPC::NextHit()
2578 {
2579   if (fHitType&2) return NextHit2();
2580   return AliDetector::NextHit();
2581 }
2582
2583 AliHit* AliTPC::FirstHit2(Int_t track)
2584 {
2585   //
2586   // Initialise the hit iterator
2587   // Return the address of the first hit for track
2588   // If track>=0 the track is read from disk
2589   // while if track<0 the first hit of the current
2590   // track is returned
2591   // 
2592   if(track>=0) {
2593     gAlice->ResetHits();
2594     gAlice->TreeH()->GetEvent(track);
2595   }
2596   //
2597   if (fTrackHits) {
2598     fTrackHits->First();
2599     return fTrackHits->GetHit();
2600   }
2601   else return 0;
2602 }
2603
2604 AliHit* AliTPC::NextHit2()
2605 {
2606   //
2607   //Return the next hit for the current track
2608
2609   if (fTrackHits) {
2610     fTrackHits->Next();
2611     return fTrackHits->GetHit();
2612   }
2613   else 
2614     return 0;
2615 }
2616
2617 void AliTPC::LoadPoints(Int_t)
2618 {
2619   //
2620   Int_t a = 0;
2621   /*  if(fHitType==1) return AliDetector::LoadPoints(a);
2622   LoadPoints2(a);
2623   */
2624   if(fHitType==1) AliDetector::LoadPoints(a);
2625   else LoadPoints2(a);
2626    
2627   // LoadPoints3(a);
2628
2629 }
2630
2631
2632 void AliTPC::RemapTrackHitIDs(Int_t *map)
2633 {
2634   if (!fTrackHits) return;
2635   AliObjectArray * arr = fTrackHits->fTrackHitsInfo;
2636   for (UInt_t i=0;i<arr->GetSize();i++){
2637     AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2638     info->fTrackID = map[info->fTrackID];
2639   }
2640   
2641 }
2642
2643
2644 //_____________________________________________________________________________
2645 void AliTPC::LoadPoints2(Int_t)
2646 {
2647   //
2648   // Store x, y, z of all hits in memory
2649   //
2650   if (fTrackHits == 0) return;
2651   //
2652   Int_t nhits = fTrackHits->GetEntriesFast();
2653   if (nhits == 0) return;
2654   Int_t tracks = gAlice->GetNtrack();
2655   if (fPoints == 0) fPoints = new TObjArray(tracks);
2656   AliHit *ahit;
2657   //
2658   Int_t *ntrk=new Int_t[tracks];
2659   Int_t *limi=new Int_t[tracks];
2660   Float_t **coor=new Float_t*[tracks];
2661   for(Int_t i=0;i<tracks;i++) {
2662     ntrk[i]=0;
2663     coor[i]=0;
2664     limi[i]=0;
2665   }
2666   //
2667   AliPoints *points = 0;
2668   Float_t *fp=0;
2669   Int_t trk;
2670   Int_t chunk=nhits/4+1;
2671   //
2672   // Loop over all the hits and store their position
2673   //
2674   ahit = FirstHit2(-1);
2675   //for (Int_t hit=0;hit<nhits;hit++) {
2676   while (ahit){
2677     //    ahit = (AliHit*)fHits->UncheckedAt(hit);
2678     trk=ahit->GetTrack();
2679     if(ntrk[trk]==limi[trk]) {
2680       //
2681       // Initialise a new track
2682       fp=new Float_t[3*(limi[trk]+chunk)];
2683       if(coor[trk]) {
2684         memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2685         delete [] coor[trk];
2686       }
2687       limi[trk]+=chunk;
2688       coor[trk] = fp;
2689     } else {
2690       fp = coor[trk];
2691     }
2692     fp[3*ntrk[trk]  ] = ahit->X();
2693     fp[3*ntrk[trk]+1] = ahit->Y();
2694     fp[3*ntrk[trk]+2] = ahit->Z();
2695     ntrk[trk]++;
2696     ahit = NextHit2();
2697   }
2698   //
2699   for(trk=0; trk<tracks; ++trk) {
2700     if(ntrk[trk]) {
2701       points = new AliPoints();
2702       points->SetMarkerColor(GetMarkerColor());
2703       points->SetMarkerSize(GetMarkerSize());
2704       points->SetDetector(this);
2705       points->SetParticle(trk);
2706       points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
2707       fPoints->AddAt(points,trk);
2708       delete [] coor[trk];
2709       coor[trk]=0;
2710     }
2711   }
2712   delete [] coor;
2713   delete [] ntrk;
2714   delete [] limi;
2715 }
2716
2717
2718 //_____________________________________________________________________________
2719 void AliTPC::LoadPoints3(Int_t)
2720 {
2721   //
2722   // Store x, y, z of all hits in memory
2723   // - only intersection point with pad row
2724   if (fTrackHits == 0) return;
2725   //
2726   Int_t nhits = fTrackHits->GetEntriesFast();
2727   if (nhits == 0) return;
2728   Int_t tracks = gAlice->GetNtrack();
2729   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2730   fPoints->Expand(2*tracks);
2731   AliHit *ahit;
2732   //
2733   Int_t *ntrk=new Int_t[tracks];
2734   Int_t *limi=new Int_t[tracks];
2735   Float_t **coor=new Float_t*[tracks];
2736   for(Int_t i=0;i<tracks;i++) {
2737     ntrk[i]=0;
2738     coor[i]=0;
2739     limi[i]=0;
2740   }
2741   //
2742   AliPoints *points = 0;
2743   Float_t *fp=0;
2744   Int_t trk;
2745   Int_t chunk=nhits/4+1;
2746   //
2747   // Loop over all the hits and store their position
2748   //
2749   ahit = FirstHit2(-1);
2750   //for (Int_t hit=0;hit<nhits;hit++) {
2751
2752   Int_t lastrow = -1;
2753   while (ahit){
2754     //    ahit = (AliHit*)fHits->UncheckedAt(hit);
2755     trk=ahit->GetTrack(); 
2756     Float_t  x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2757     Int_t    index[3]={1,((AliTPChit*)ahit)->fSector,0};
2758     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2759     if (currentrow!=lastrow){
2760       lastrow = currentrow;
2761       //later calculate intersection point           
2762       if(ntrk[trk]==limi[trk]) {
2763         //
2764         // Initialise a new track
2765         fp=new Float_t[3*(limi[trk]+chunk)];
2766         if(coor[trk]) {
2767           memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2768           delete [] coor[trk];
2769         }
2770         limi[trk]+=chunk;
2771         coor[trk] = fp;
2772       } else {
2773         fp = coor[trk];
2774       }
2775       fp[3*ntrk[trk]  ] = ahit->X();
2776       fp[3*ntrk[trk]+1] = ahit->Y();
2777       fp[3*ntrk[trk]+2] = ahit->Z();
2778       ntrk[trk]++;
2779     }
2780     ahit = NextHit2();
2781   }
2782   
2783   //
2784   for(trk=0; trk<tracks; ++trk) {
2785     if(ntrk[trk]) {
2786       points = new AliPoints();
2787       points->SetMarkerColor(GetMarkerColor()+1);
2788       points->SetMarkerStyle(5);
2789       points->SetMarkerSize(0.2);
2790       points->SetDetector(this);
2791       points->SetParticle(trk);
2792       //      points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20);
2793       points->SetPolyMarker(ntrk[trk],coor[trk],30);
2794       fPoints->AddAt(points,tracks+trk);
2795       delete [] coor[trk];
2796       coor[trk]=0;
2797     }
2798   }
2799   delete [] coor;
2800   delete [] ntrk;
2801   delete [] limi;
2802 }
2803
2804
2805
2806 void AliTPC::FindTrackHitsIntersection(TClonesArray * arr)
2807 {
2808
2809   //
2810   //fill clones array with intersection of current point with the
2811   //middle of the row
2812   Int_t sector;
2813   Int_t ipart;
2814   
2815   const Int_t kcmaxhits=30000;
2816   TVector * xxxx = new TVector(kcmaxhits*4);
2817   TVector & xxx = *xxxx;
2818   Int_t maxhits = kcmaxhits;
2819       
2820   //
2821   AliTPChit * tpcHit=0;
2822   tpcHit = (AliTPChit*)FirstHit2(-1);
2823   Int_t currentIndex=0;
2824   Int_t lastrow=-1;  //last writen row
2825
2826   while (tpcHit){
2827     if (tpcHit==0) continue;
2828     sector=tpcHit->fSector; // sector number
2829     ipart=tpcHit->Track();
2830     
2831     //find row number
2832     
2833     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
2834     Int_t    index[3]={1,sector,0};
2835     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;       
2836     if (currentrow<0) continue;
2837     if (lastrow<0) lastrow=currentrow;
2838     if (currentrow==lastrow){
2839       if ( currentIndex>=maxhits){
2840         maxhits+=kcmaxhits;
2841         xxx.ResizeTo(4*maxhits);
2842       }     
2843       xxx(currentIndex*4)=x[0];
2844       xxx(currentIndex*4+1)=x[1];
2845       xxx(currentIndex*4+2)=x[2];       
2846       xxx(currentIndex*4+3)=tpcHit->fQ;
2847       currentIndex++;   
2848     }
2849     else 
2850       if (currentIndex>2){
2851         Float_t sumx=0;
2852         Float_t sumx2=0;
2853         Float_t sumx3=0;
2854         Float_t sumx4=0;
2855         Float_t sumy=0;
2856         Float_t sumxy=0;
2857         Float_t sumx2y=0;
2858         Float_t sumz=0;
2859         Float_t sumxz=0;
2860         Float_t sumx2z=0;
2861         Float_t sumq=0;
2862         for (Int_t index=0;index<currentIndex;index++){
2863           Float_t x,x2,x3,x4;
2864           x=x2=x3=x4=xxx(index*4);
2865           x2*=x;
2866           x3*=x2;
2867           x4*=x3;
2868           sumx+=x;
2869           sumx2+=x2;
2870           sumx3+=x3;
2871           sumx4+=x4;
2872           sumy+=xxx(index*4+1);
2873           sumxy+=xxx(index*4+1)*x;
2874           sumx2y+=xxx(index*4+1)*x2;
2875           sumz+=xxx(index*4+2);
2876           sumxz+=xxx(index*4+2)*x;
2877           sumx2z+=xxx(index*4+2)*x2;     
2878           sumq+=xxx(index*4+3);
2879         }
2880         Float_t centralPad = (fTPCParam->GetNPads(sector,lastrow)-1)/2;
2881         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
2882           sumx2*(sumx*sumx3-sumx2*sumx2);
2883         
2884         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
2885           sumx2*(sumxy*sumx3-sumx2y*sumx2);
2886         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
2887           sumx2*(sumxz*sumx3-sumx2z*sumx2);
2888         
2889         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
2890           sumx2*(sumx*sumx2y-sumx2*sumxy);
2891         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
2892           sumx2*(sumx*sumx2z-sumx2*sumxz);
2893         
2894         Float_t y=detay/det+centralPad;
2895         Float_t z=detaz/det;    
2896         Float_t by=detby/det; //y angle
2897         Float_t bz=detbz/det; //z angle
2898         sumy/=Float_t(currentIndex);
2899         sumz/=Float_t(currentIndex);
2900         
2901         AliComplexCluster cl;
2902         cl.fX=z;
2903         cl.fY=y;
2904         cl.fQ=sumq;
2905         cl.fSigmaX2=bz;
2906         cl.fSigmaY2=by;
2907         cl.fTracks[0]=ipart;
2908         
2909         AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow));
2910         if (row!=0) row->InsertCluster(&cl);
2911         currentIndex=0;
2912         lastrow=currentrow;
2913       } //end of calculating cluster for given row
2914                 
2915   } // end of loop over hits
2916   xxxx->Delete();
2917
2918 }
2919 //_______________________________________________________________________________
2920 void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
2921 {
2922   // produces rec points from digits and writes them on the root file
2923   // AliTPCclusters.root
2924
2925   TDirectory *cwd = gDirectory;
2926
2927
2928   AliTPCParam *dig=(AliTPCParam *)gDirectory->Get("75x40_100x60");
2929   SetParam(dig);
2930   cout<<"AliTPC::Digits2Reco: TPC parameteres have been set"<<endl; 
2931   TFile *out;
2932   if(!gSystem->Getenv("CONFIG_FILE")){
2933     out=TFile::Open("AliTPCclusters.root","recreate");
2934   }
2935   else {
2936     const char *tmp1;
2937     const char *tmp2;
2938     char tmp3[80];
2939     tmp1=gSystem->Getenv("CONFIG_FILE_PREFIX");
2940     tmp2=gSystem->Getenv("CONFIG_OUTDIR");
2941     sprintf(tmp3,"%s%s/AliTPCclusters.root",tmp1,tmp2);
2942     out=TFile::Open(tmp3,"recreate");
2943   }
2944
2945   TStopwatch timer;
2946   cout<<"AliTPC::Digits2Reco - determination of rec points begins"<<endl;
2947   timer.Start();
2948   cwd->cd();
2949   for(Int_t iev=firstevent;iev<lastevent+1;iev++){
2950
2951     printf("Processing event %d\n",iev);
2952     Digits2Clusters(out,iev);
2953   }
2954   cout<<"AliTPC::Digits2Reco - determination of rec points ended"<<endl;
2955   timer.Stop();
2956   timer.Print();
2957   out->Close();
2958   cwd->cd(); 
2959
2960
2961 }