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