4a60d9e7c5166c095dc012b13270d8423f6dfde4
[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.19.2.4  2000/06/26 07:39:42  kowal2
19 Changes to obey the coding rules
20
21 Revision 1.19.2.3  2000/06/25 08:38:41  kowal2
22 Splitted from AliTPCtracking
23
24 Revision 1.19.2.2  2000/06/16 12:59:28  kowal2
25 Changed parameter settings
26
27 Revision 1.19.2.1  2000/06/09 07:15:07  kowal2
28
29 Defaults loaded automatically (hard-wired)
30 Optional parameters can be set via macro called in the constructor
31
32 Revision 1.19  2000/04/18 19:00:59  fca
33 Small bug fixes to TPC files
34
35 Revision 1.18  2000/04/17 09:37:33  kowal2
36 removed obsolete AliTPCDigitsDisplay.C
37
38 Revision 1.17.2.2  2000/04/10 08:15:12  kowal2
39
40 New, experimental data structure from M. Ivanov
41 New tracking algorithm
42 Different pad geometry for different sectors
43 Digitization rewritten
44
45 Revision 1.17.2.1  2000/04/10 07:56:53  kowal2
46 Not used anymore - removed
47
48 Revision 1.17  2000/01/19 17:17:30  fca
49 Introducing a list of lists of hits -- more hits allowed for detector now
50
51 Revision 1.16  1999/11/05 09:29:23  fca
52 Accept only signals > 0
53
54 Revision 1.15  1999/10/08 06:26:53  fca
55 Removed ClustersIndex - not used anymore
56
57 Revision 1.14  1999/09/29 09:24:33  fca
58 Introduction of the Copyright and cvs Log
59
60 */
61
62 ///////////////////////////////////////////////////////////////////////////////
63 //                                                                           //
64 //  Time Projection Chamber                                                  //
65 //  This class contains the basic functions for the Time Projection Chamber  //
66 //  detector. Functions specific to one particular geometry are              //
67 //  contained in the derived classes                                         //
68 //                                                                           //
69 //Begin_Html
70 /*
71 <img src="picts/AliTPCClass.gif">
72 */
73 //End_Html
74 //                                                                           //
75 //                                                                          //
76 ///////////////////////////////////////////////////////////////////////////////
77
78 //
79
80 #include <TMath.h>
81 #include <TRandom.h>
82 #include <TVector.h>
83 #include <TMatrix.h>
84 #include <TGeometry.h>
85 #include <TNode.h>
86 #include <TTUBS.h>
87 #include <TObjectTable.h>
88 #include "TParticle.h"
89 #include "AliTPC.h"
90 #include <TFile.h>       
91 #include "AliRun.h"
92 #include <iostream.h>
93 #include <fstream.h>
94 #include "AliMC.h"
95
96
97 #include "AliTPCParamSR.h"
98 #include "AliTPCPRF2D.h"
99 #include "AliTPCRF1D.h"
100 #include "AliDigits.h"
101 #include "AliSimDigits.h"
102
103 #include "AliTPCDigitsArray.h"
104 #include "AliCluster.h"
105 #include "AliClusters.h"
106 #include "AliTPCClustersRow.h"
107 #include "AliTPCClustersArray.h"
108
109 #include "AliTPCcluster.h"
110 #include "AliTPCclusterer.h"
111 #include "AliTPCtracker.h"
112
113 #include <TInterpreter.h>
114
115 ClassImp(AliTPC) 
116
117 //_____________________________________________________________________________
118 AliTPC::AliTPC()
119 {
120   //
121   // Default constructor
122   //
123   fIshunt   = 0;
124   fHits     = 0;
125   fDigits   = 0;
126   fNsectors = 0;
127   //MI changes
128   fDigitsArray = 0;
129   fClustersArray = 0;
130   fTPCParam=0;
131 }
132  
133 //_____________________________________________________________________________
134 AliTPC::AliTPC(const char *name, const char *title)
135       : AliDetector(name,title)
136 {
137   //
138   // Standard constructor
139   //
140
141   //
142   // Initialise arrays of hits and digits 
143   fHits     = new TClonesArray("AliTPChit",  176);
144   gAlice->AddHitList(fHits);
145   //MI change  
146   fDigitsArray = 0;
147   fClustersArray= 0;
148   //
149   // Initialise counters
150   fNsectors = 0;
151
152   //
153   fIshunt     =  0;
154   //
155   // Initialise color attributes
156   SetMarkerColor(kYellow);
157
158   //
159   //  Set TPC parameters
160   //
161
162   if (!strcmp(title,"Default")) {  
163      fTPCParam = new AliTPCParamSR;
164   } else {
165     cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
166     fTPCParam=0;
167   }
168
169 }
170
171 //_____________________________________________________________________________
172 AliTPC::~AliTPC()
173 {
174   //
175   // TPC destructor
176   //
177
178   fIshunt   = 0;
179   delete fHits;
180   delete fDigits;
181   delete fTPCParam;
182 }
183
184 //_____________________________________________________________________________
185 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
186 {
187   //
188   // Add a hit to the list
189   //
190   TClonesArray &lhits = *fHits;
191   new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
192 }
193  
194 //_____________________________________________________________________________
195 void AliTPC::BuildGeometry()
196 {
197
198   //
199   // Build TPC ROOT TNode geometry for the event display
200   //
201   TNode *nNode, *nTop;
202   TTUBS *tubs;
203   Int_t i;
204   const int kColorTPC=19;
205   char name[5], title[25];
206   const Double_t kDegrad=TMath::Pi()/180;
207   const Double_t kRaddeg=180./TMath::Pi();
208
209
210   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
211   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
212
213   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
214   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
215
216   Int_t nLo = fTPCParam->GetNInnerSector()/2;
217   Int_t nHi = fTPCParam->GetNOuterSector()/2;  
218
219   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
220   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
221   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
222   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);  
223
224
225   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
226   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
227
228   Double_t rl,ru;
229   
230
231   //
232   // Get ALICE top node
233   //
234
235   nTop=gAlice->GetGeometry()->GetNode("alice");
236
237   //  inner sectors
238
239   rl = fTPCParam->GetInnerRadiusLow();
240   ru = fTPCParam->GetInnerRadiusUp();
241  
242
243   for(i=0;i<nLo;i++) {
244     sprintf(name,"LS%2.2d",i);
245     name[4]='\0';
246     sprintf(title,"TPC low sector %3d",i);
247     title[24]='\0';
248     
249     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
250                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
251     tubs->SetNumberOfDivisions(1);
252     nTop->cd();
253     nNode = new TNode(name,title,name,0,0,0,"");
254     nNode->SetLineColor(kColorTPC);
255     fNodes->Add(nNode);
256   }
257
258   // Outer sectors
259
260   rl = fTPCParam->GetOuterRadiusLow();
261   ru = fTPCParam->GetOuterRadiusUp();
262
263   for(i=0;i<nHi;i++) {
264     sprintf(name,"US%2.2d",i);
265     name[4]='\0';
266     sprintf(title,"TPC upper sector %d",i);
267     title[24]='\0';
268     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
269                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
270     tubs->SetNumberOfDivisions(1);
271     nTop->cd();
272     nNode = new TNode(name,title,name,0,0,0,"");
273     nNode->SetLineColor(kColorTPC);
274     fNodes->Add(nNode);
275   }
276
277 }    
278
279 //_____________________________________________________________________________
280 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
281 {
282   //
283   // Calculate distance from TPC to mouse on the display
284   // Dummy procedure
285   //
286   return 9999;
287 }
288
289 void AliTPC::Clusters2Tracks(TFile *of) {
290   //-----------------------------------------------------------------
291   // This is a track finder.
292   //-----------------------------------------------------------------
293   AliTPCtracker::Clusters2Tracks(fTPCParam,of);
294 }
295
296 //_____________________________________________________________________________
297 void AliTPC::CreateMaterials()
298 {
299   //-----------------------------------------------
300   // Create Materials for for TPC
301   //-----------------------------------------------
302
303   //-----------------------------------------------------------------
304   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
305   //-----------------------------------------------------------------
306
307   Int_t iSXFLD=gAlice->Field()->Integ();
308   Float_t sXMGMX=gAlice->Field()->Max();
309
310   Float_t amat[5]; // atomic numbers
311   Float_t zmat[5]; // z
312   Float_t wmat[5]; // proportions
313
314   Float_t density;
315
316   //  ********************* Gases *******************
317
318   //--------------------------------------------------------------
319   // pure gases
320   //--------------------------------------------------------------
321
322   // Ne
323
324
325   Float_t aNe = 20.18;
326   Float_t zNe = 10.;
327   
328   density = 0.0009;
329
330   AliMaterial(20,"Ne",aNe,zNe,density,999.,999.);
331
332   // Ar
333
334   Float_t aAr = 39.948;
335   Float_t zAr = 18.;
336
337   density = 0.001782;
338  
339   AliMaterial(21,"Ar",aAr,zAr,density,999.,999.);
340
341   Float_t aPure[2];
342   
343   aPure[0] = aNe;
344   aPure[1] = aAr;
345   
346
347   //--------------------------------------------------------------
348   // gases - compounds
349   //--------------------------------------------------------------
350
351   Float_t amol[3];
352
353   //  CO2
354
355   amat[0]=12.011;
356   amat[1]=15.9994;
357
358   zmat[0]=6.;
359   zmat[1]=8.;
360
361   wmat[0]=1.;
362   wmat[1]=2.;
363
364   density=0.001977;
365
366   amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
367
368   AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
369
370   // CF4
371
372   amat[0]=12.011;
373   amat[1]=18.998;
374
375   zmat[0]=6.;
376   zmat[1]=9.;
377  
378   wmat[0]=1.;
379   wmat[1]=4.;
380  
381   density=0.003034;
382
383   amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
384
385   AliMixture(11,"CF4",amat,zmat,density,-2,wmat); 
386
387   // CH4
388
389   amat[0]=12.011;
390   amat[1]=1.;
391
392   zmat[0]=6.;
393   zmat[1]=1.;
394
395   wmat[0]=1.;
396   wmat[1]=4.;
397
398   density=0.000717;
399
400   amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
401
402   AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
403
404   //----------------------------------------------------------------
405   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
406   //----------------------------------------------------------------
407  
408   char namate[21];
409  
410   density = 0.;
411   Float_t am=0;
412   Int_t nc;
413
414   Float_t a,z,rho,absl,x0,buf[1];
415   Int_t nbuf;
416
417   for(nc = 0;nc<fNoComp;nc++)
418     {
419     
420       // retrive material constants
421       
422       gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,x0,absl,buf,nbuf);
423
424       amat[nc] = a;
425       zmat[nc] = z;
426
427       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
428  
429       am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? aPure[nnc] : amol[nnc]); 
430       density += fMixtProp[nc]*rho;  // density of the mixture
431       
432     }
433
434   // mixture proportions by weight!
435
436   for(nc = 0;nc<fNoComp;nc++)
437     {
438
439       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
440
441       wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? aPure[nnc] : amol[nnc])/am;
442
443     }  
444   
445   AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
446   AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
447   AliMixture(33,"Drift gas 3",amat,zmat,density,fNoComp,wmat); 
448
449   AliMedium(2, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
450   AliMedium(3, "Drift gas 2", 32, 0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
451   AliMedium(4, "Drift gas 3", 33, 1, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
452
453   // Air 
454
455   AliMaterial(24, "Air", 14.61, 7.3, .001205, 30420., 67500.);
456
457   AliMedium(24, "Air", 24, 0, iSXFLD, sXMGMX, 10., .1, .1, .1, .1);
458
459   //----------------------------------------------------------------------
460   //               solid materials
461   //----------------------------------------------------------------------
462
463   // Al
464
465   AliMaterial(30, "Al", 26.98, 13., 2.7, 8.9, 37.2);
466
467   AliMedium(0, "Al",30, 0, iSXFLD, sXMGMX, 10., .1, .1, .1,   .1);
468
469   // Si
470
471   AliMaterial(31, "Si", 28.086, 14.,2.33, 9.36, 999.);
472
473   AliMedium(7, "Al",31, 0, iSXFLD, sXMGMX, 10., .1, .1, .1,   .1);
474   
475
476   // Mylar C5H4O2
477
478   amat[0]=12.011;
479   amat[1]=1.;
480   amat[2]=15.9994;
481
482   zmat[0]=6.;
483   zmat[1]=1.;
484   zmat[2]=8.;
485
486   wmat[0]=5.;
487   wmat[1]=4.;
488   wmat[2]=2.; 
489
490   density = 1.39;
491   
492   AliMixture(32, "Mylar",amat,zmat,density,-3,wmat);
493
494   AliMedium(5, "Mylar",32, 0, iSXFLD, sXMGMX, 10., .1, .1, .001, .01);
495
496
497
498
499   // Carbon (normal)
500
501   AliMaterial(33,"C normal",12.011,6.,2.265,18.8,999.);
502
503   AliMedium(6,"C normal",33,0, iSXFLD, sXMGMX, 10., .1, .1, .001, .01);
504
505   // G10 for inner and outr field cage
506   // G10 is 60% SiO2 + 40% epoxy, right now I use A and Z for SiO2
507
508   Float_t rhoFactor;
509
510   amat[0]=28.086;
511   amat[1]=15.9994;
512
513   zmat[0]=14.;
514   zmat[1]=8.;
515
516   wmat[0]=1.;
517   wmat[1]=2.;
518
519   density = 1.7;
520   
521
522   AliMixture(34,"G10 aux.",amat,zmat,density,-2,wmat);
523
524
525   gMC->Gfmate((*fIdmate)[34],namate,a,z,rho,x0,absl,buf,nbuf);
526
527   Float_t thickX0 = 0.0052; // field cage in X0 units
528   
529   Float_t thick = 2.; // in cm
530
531   x0=19.4; // G10 
532
533   rhoFactor = x0*thickX0/thick;
534   density = rho*rhoFactor;
535
536   AliMaterial(35,"G10-fc",a,z,density,999.,999.);
537
538   AliMedium(8,"G10-fc",35,0, iSXFLD, sXMGMX, 10., .1, .1, .001, .01);
539
540   thickX0 = 0.0027; // inner vessel (eta <0.9)
541   thick=0.5;
542   rhoFactor = x0*thickX0/thick;
543   density = rho*rhoFactor;
544
545   AliMaterial(36,"G10-iv",a,z,density,999.,999.);  
546
547   AliMedium(9,"G10-iv",36,0, iSXFLD, sXMGMX, 10., .1, .1, .001, .01);
548
549   //  Carbon fibre  
550   
551   gMC->Gfmate((*fIdmate)[33],namate,a,z,rho,x0,absl,buf,nbuf);
552
553   thickX0 = 0.0133; // outer vessel
554   thick=3.0;
555   rhoFactor = x0*thickX0/thick;
556   density = rho*rhoFactor;
557
558
559   AliMaterial(37,"C-ov",a,z,density,999.,999.);
560
561   AliMedium(10,"C-ov",37,0, iSXFLD, sXMGMX, 10., .1, .1, .001, .01);  
562
563   thickX0=0.015; // inner vessel (cone, eta > 0.9)
564   thick=1.5;
565   rhoFactor = x0*thickX0/thick;
566   density = rho*rhoFactor;
567
568   AliMaterial(38,"C-ivc",a,z,density,999.,999.);
569
570   AliMedium(11,"C-ivc",38,0, iSXFLD, sXMGMX, 10., .1, .1, .001, .01);
571
572   //
573
574   AliMedium(12,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
575     
576 }
577
578
579 void AliTPC::Digits2Clusters(TFile *of)
580 {
581   //-----------------------------------------------------------------
582   // This is a simple cluster finder.
583   //-----------------------------------------------------------------
584   AliTPCclusterer::Digits2Clusters(fTPCParam,of);
585 }
586
587 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
588 extern Double_t SigmaZ2(Double_t, Double_t);
589 //_____________________________________________________________________________
590 void AliTPC::Hits2Clusters(TFile *of)
591 {
592   //--------------------------------------------------------
593   // TPC simple cluster generator from hits
594   // obtained from the TPC Fast Simulator
595   // The point errors are taken from the parametrization
596   //--------------------------------------------------------
597
598   //-----------------------------------------------------------------
599   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
600   //-----------------------------------------------------------------
601   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
602   // Jouri.Belikov@cern.ch
603   //----------------------------------------------------------------
604   
605   /////////////////////////////////////////////////////////////////////////////
606   //
607   //---------------------------------------------------------------------
608   //   ALICE TPC Cluster Parameters
609   //--------------------------------------------------------------------
610        
611   
612
613   // Cluster width in rphi
614   const Float_t kACrphi=0.18322;
615   const Float_t kBCrphi=0.59551e-3;
616   const Float_t kCCrphi=0.60952e-1;
617   // Cluster width in z
618   const Float_t kACz=0.19081;
619   const Float_t kBCz=0.55938e-3;
620   const Float_t kCCz=0.30428;
621
622   TDirectory *savedir=gDirectory; 
623
624   if (!of->IsOpen()) {
625      cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
626      return;
627   }
628
629    if(fTPCParam == 0){
630      printf("AliTPCParam MUST be created firstly\n");
631      return;
632    }
633
634   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
635   //
636   TParticle *particle; // pointer to a given particle
637   AliTPChit *tpcHit; // pointer to a sigle TPC hit
638   TClonesArray *particles; //pointer to the particle list
639   Int_t sector,nhits;
640   Int_t ipart;
641   Float_t xyz[5];
642   Float_t pl,pt,tanth,rpad,ratio;
643   Float_t cph,sph;
644   
645   //---------------------------------------------------------------
646   //  Get the access to the tracks 
647   //---------------------------------------------------------------
648   
649   TTree *tH = gAlice->TreeH();
650   Stat_t ntracks = tH->GetEntries();
651   particles=gAlice->Particles();
652
653   //Switch to the output file
654   of->cd();
655
656   fTPCParam->Write(fTPCParam->GetTitle());
657   AliTPCClustersArray carray;
658   carray.Setup(fTPCParam);
659   carray.SetClusterType("AliTPCcluster");
660   carray.MakeTree();
661
662   Int_t nclusters=0; //cluster counter
663   
664   //------------------------------------------------------------
665   // Loop over all sectors (72 sectors for 20 deg
666   // segmentation for both lower and upper sectors)
667   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
668   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
669   //
670   // First cluster for sector 0 starts at "0"
671   //------------------------------------------------------------
672    
673   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
674     //MI change
675     fTPCParam->AdjustCosSin(isec,cph,sph);
676     
677     //------------------------------------------------------------
678     // Loop over tracks
679     //------------------------------------------------------------
680     
681     for(Int_t track=0;track<ntracks;track++){
682       ResetHits();
683       tH->GetEvent(track);
684       //
685       //  Get number of the TPC hits
686       //
687       nhits=fHits->GetEntriesFast();
688       //
689       // Loop over hits
690       //
691       for(Int_t hit=0;hit<nhits;hit++){
692         tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
693         if (tpcHit->fQ == 0.) continue; //information about track (I.Belikov)
694         sector=tpcHit->fSector; // sector number
695         if(sector != isec) continue; //terminate iteration
696         ipart=tpcHit->fTrack;
697         particle=(TParticle*)particles->UncheckedAt(ipart);
698         pl=particle->Pz();
699         pt=particle->Pt();
700         if(pt < 1.e-9) pt=1.e-9;
701         tanth=pl/pt;
702         tanth = TMath::Abs(tanth);
703         rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
704         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
705
706         //   space-point resolutions
707         
708         sigmaRphi=SigmaY2(rpad,tanth,pt);
709         sigmaZ   =SigmaZ2(rpad,tanth   );
710         
711         //   cluster widths
712         
713         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
714         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
715         
716         // temporary protection
717         
718         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
719         if(sigmaZ < 0.) sigmaZ=0.4e-3;
720         if(clRphi < 0.) clRphi=2.5e-3;
721         if(clZ < 0.) clZ=2.5e-5;
722         
723         //
724         
725         //
726         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
727         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
728         //
729         Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
730         Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
731         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
732           Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
733           fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
734           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
735           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim; 
736         xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigmaZ)); // z
737           if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->fZ; 
738         xyz[2]=tpcHit->fQ;                                     // q
739         xyz[3]=sigmaRphi;                                     // fSigmaY2
740         xyz[4]=sigmaZ;                                        // fSigmaZ2
741
742         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
743         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);     
744
745         Int_t tracks[3]={tpcHit->fTrack, -1, -1};
746         AliTPCcluster cluster(xyz,tracks);
747
748         clrow->InsertCluster(&cluster); nclusters++;
749
750       } // end of loop over hits
751
752     }   // end of loop over tracks
753
754     Int_t nrows=fTPCParam->GetNRow(isec);
755     for (Int_t irow=0; irow<nrows; irow++) {
756         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
757         if (!clrow) continue;
758         carray.StoreRow(isec,irow);
759         carray.ClearRow(isec,irow);
760     }
761
762   } // end of loop over sectors  
763
764   cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
765
766   carray.GetTree()->Write();
767
768   savedir->cd(); //switch back to the input file
769   
770 } // end of function
771
772 //_________________________________________________________________
773 void AliTPC::Hits2ExactClustersSector(Int_t isec)
774 {
775   //--------------------------------------------------------
776   //calculate exact cross point of track and given pad row
777   //resulting values are expressed in "digit" coordinata
778   //--------------------------------------------------------
779
780   //-----------------------------------------------------------------
781   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
782   //-----------------------------------------------------------------
783   //
784   if (fClustersArray==0){    
785     return;
786   }
787   //
788   TParticle *particle; // pointer to a given particle
789   AliTPChit *tpcHit; // pointer to a sigle TPC hit
790   TClonesArray *particles; //pointer to the particle list
791   Int_t sector,nhits;
792   Int_t ipart;
793   const Int_t kcmaxhits=30000;
794   TVector * xxxx = new TVector(kcmaxhits*4);
795   TVector & xxx = *xxxx;
796   Int_t maxhits = kcmaxhits;
797   //construct array for each padrow
798   for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++) 
799     fClustersArray->CreateRow(isec,i);
800   
801   //---------------------------------------------------------------
802   //  Get the access to the tracks 
803   //---------------------------------------------------------------
804   
805   TTree *tH = gAlice->TreeH();
806   Stat_t ntracks = tH->GetEntries();
807   particles=gAlice->Particles();
808   Int_t npart = particles->GetEntriesFast();
809     
810   //------------------------------------------------------------
811   // Loop over tracks
812   //------------------------------------------------------------
813   
814   for(Int_t track=0;track<ntracks;track++){
815     ResetHits();
816     tH->GetEvent(track);
817     //
818     //  Get number of the TPC hits and a pointer
819     //  to the particles
820     //
821     nhits=fHits->GetEntriesFast();
822     //
823     // Loop over hits
824     //
825     Int_t currentIndex=0;
826     Int_t lastrow=-1;  //last writen row
827     for(Int_t hit=0;hit<nhits;hit++){
828       tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
829       if (tpcHit==0) continue;
830       sector=tpcHit->fSector; // sector number
831       if(sector != isec) continue; 
832       ipart=tpcHit->fTrack;
833       if (ipart<npart) particle=(TParticle*)particles->UncheckedAt(ipart);
834       
835       //find row number
836
837       Float_t  x[3]={tpcHit->fX,tpcHit->fY,tpcHit->fZ};
838       Int_t    index[3]={1,isec,0};
839       Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;     
840       if (currentrow<0) continue;
841       if (lastrow<0) lastrow=currentrow;
842       if (currentrow==lastrow){
843         if ( currentIndex>=maxhits){
844           maxhits+=kcmaxhits;
845           xxx.ResizeTo(4*maxhits);
846         }     
847         xxx(currentIndex*4)=x[0];
848         xxx(currentIndex*4+1)=x[1];
849         xxx(currentIndex*4+2)=x[2];     
850         xxx(currentIndex*4+3)=tpcHit->fQ;
851         currentIndex++; 
852       }
853       else 
854         if (currentIndex>2){
855           Float_t sumx=0;
856           Float_t sumx2=0;
857           Float_t sumx3=0;
858           Float_t sumx4=0;
859           Float_t sumy=0;
860           Float_t sumxy=0;
861           Float_t sumx2y=0;
862           Float_t sumz=0;
863           Float_t sumxz=0;
864           Float_t sumx2z=0;
865           Float_t sumq=0;
866           for (Int_t index=0;index<currentIndex;index++){
867             Float_t x,x2,x3,x4;
868             x=x2=x3=x4=xxx(index*4);
869             x2*=x;
870             x3*=x2;
871             x4*=x3;
872             sumx+=x;
873             sumx2+=x2;
874             sumx3+=x3;
875             sumx4+=x4;
876             sumy+=xxx(index*4+1);
877             sumxy+=xxx(index*4+1)*x;
878             sumx2y+=xxx(index*4+1)*x2;
879             sumz+=xxx(index*4+2);
880             sumxz+=xxx(index*4+2)*x;
881             sumx2z+=xxx(index*4+2)*x2;   
882             sumq+=xxx(index*4+3);
883           }
884           Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
885           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
886             sumx2*(sumx*sumx3-sumx2*sumx2);
887           
888           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
889             sumx2*(sumxy*sumx3-sumx2y*sumx2);
890           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
891             sumx2*(sumxz*sumx3-sumx2z*sumx2);
892           
893           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
894             sumx2*(sumx*sumx2y-sumx2*sumxy);
895           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
896             sumx2*(sumx*sumx2z-sumx2*sumxz);
897           
898           Float_t y=detay/det+centralPad;
899           Float_t z=detaz/det;  
900           Float_t by=detby/det; //y angle
901           Float_t bz=detbz/det; //z angle
902           sumy/=Float_t(currentIndex);
903           sumz/=Float_t(currentIndex);
904           AliCluster cl;
905           cl.fX=z;
906           cl.fY=y;
907           cl.fQ=sumq;
908           cl.fSigmaX2=bz;
909           cl.fSigmaY2=by;
910           cl.fTracks[0]=ipart;
911           
912           AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
913           if (row!=0) row->InsertCluster(&cl);
914           currentIndex=0;
915           lastrow=currentrow;
916         } //end of calculating cluster for given row
917         
918         
919         
920     } // end of loop over hits
921   }   // end of loop over tracks 
922   //write padrows to tree 
923   for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
924     fClustersArray->StoreRow(isec,ii);    
925     fClustersArray->ClearRow(isec,ii);        
926   }
927   xxxx->Delete();
928  
929 }
930
931 //__________________________________________________________________  
932 void AliTPC::Hits2Digits()  
933
934  //----------------------------------------------------
935  // Loop over all sectors
936  //----------------------------------------------------
937
938   if(fTPCParam == 0){
939     printf("AliTPCParam MUST be created firstly\n");
940     return;
941   } 
942
943  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
944
945 }
946
947
948 //_____________________________________________________________________________
949 void AliTPC::Hits2DigitsSector(Int_t isec)
950 {
951   //-------------------------------------------------------------------
952   // TPC conversion from hits to digits.
953   //------------------------------------------------------------------- 
954
955   //-----------------------------------------------------------------
956   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
957   //-----------------------------------------------------------------
958
959   //-------------------------------------------------------
960   //  Get the access to the track hits
961   //-------------------------------------------------------
962
963
964   TTree *tH = gAlice->TreeH(); // pointer to the hits tree
965   Stat_t ntracks = tH->GetEntries();
966
967   if( ntracks > 0){
968
969   //------------------------------------------- 
970   //  Only if there are any tracks...
971   //-------------------------------------------
972
973     TObjArray **row;
974     
975       printf("*** Processing sector number %d ***\n",isec);
976
977       Int_t nrows =fTPCParam->GetNRow(isec);
978
979       row= new TObjArray* [nrows];
980     
981       MakeSector(isec,nrows,tH,ntracks,row);
982
983       //--------------------------------------------------------
984       //   Digitize this sector, row by row
985       //   row[i] is the pointer to the TObjArray of TVectors,
986       //   each one containing electrons accepted on this
987       //   row, assigned into tracks
988       //--------------------------------------------------------
989
990       Int_t i;
991
992       if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree();
993
994       for (i=0;i<nrows;i++){
995
996         AliDigits * dig = fDigitsArray->CreateRow(isec,i); 
997
998         DigitizeRow(i,isec,row);
999
1000         fDigitsArray->StoreRow(isec,i);
1001
1002         Int_t ndig = dig->GetDigitSize(); 
1003  
1004         printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1005         
1006         fDigitsArray->ClearRow(isec,i);  
1007
1008    
1009        } // end of the sector digitization
1010
1011       for(i=0;i<nrows;i++){
1012         row[i]->Delete();     
1013       }
1014       
1015        delete [] row; // delete the array of pointers to TObjArray-s
1016         
1017   } // ntracks >0
1018
1019 } // end of Hits2DigitsSector
1020
1021
1022 //_____________________________________________________________________________
1023 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1024 {
1025   //-----------------------------------------------------------
1026   // Single row digitization, coupling from the neighbouring
1027   // rows taken into account
1028   //-----------------------------------------------------------
1029
1030   //-----------------------------------------------------------------
1031   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1032   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1033   //-----------------------------------------------------------------
1034  
1035
1036   Float_t zerosup = fTPCParam->GetZeroSup();
1037   Int_t nrows =fTPCParam->GetNRow(isec);
1038   fCurrentIndex[1]= isec;
1039   
1040
1041   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1042   Int_t nofTbins = fTPCParam->GetMaxTBin();
1043   Int_t indexRange[4];
1044   //
1045   //  Integrated signal for this row
1046   //  and a single track signal
1047   //    
1048   TMatrix *m1   = new TMatrix(0,nofPads,0,nofTbins); // integrated
1049   TMatrix *m2   = new TMatrix(0,nofPads,0,nofTbins); // single
1050   //
1051   TMatrix &total  = *m1;
1052
1053   //  Array of pointers to the label-signal list
1054
1055   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1056   Float_t  **pList = new Float_t* [nofDigits]; 
1057
1058   Int_t lp;
1059   Int_t i1;   
1060   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1061   //
1062   //calculate signal 
1063   //
1064   Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1065   Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1066   for (Int_t row= row1;row<=row2;row++){
1067     Int_t nTracks= rows[row]->GetEntries();
1068     for (i1=0;i1<nTracks;i1++){
1069       fCurrentIndex[2]= row;
1070       fCurrentIndex[3]=irow;
1071       if (row==irow){
1072         m2->Zero();  // clear single track signal matrix
1073         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange); 
1074         GetList(trackLabel,nofPads,m2,indexRange,pList);
1075       }
1076       else   GetSignal(rows[row],i1,0,m1,indexRange);
1077     }
1078   }
1079          
1080   Int_t tracks[3];
1081
1082   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1083   for(Int_t ip=0;ip<nofPads;ip++){
1084     for(Int_t it=0;it<nofTbins;it++){
1085
1086       Float_t q = total(ip,it);
1087
1088       Int_t gi =it*nofPads+ip; // global index
1089
1090       q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); 
1091
1092       q = (Int_t)q;
1093
1094       if(q <=zerosup) continue; // do not fill zeros
1095       if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();  // saturation
1096
1097       //
1098       //  "real" signal or electronic noise (list = -1)?
1099       //    
1100
1101       for(Int_t j1=0;j1<3;j1++){
1102         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1103       }
1104
1105 //Begin_Html
1106 /*
1107   <A NAME="AliDigits"></A>
1108   using of AliDigits object
1109 */
1110 //End_Html
1111       dig->SetDigitFast((Short_t)q,it,ip);
1112       if (fDigitsArray->IsSimulated())
1113         {
1114          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1115          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1116          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1117         }
1118      
1119     
1120     } // end of loop over time buckets
1121   }  // end of lop over pads 
1122
1123   //
1124   //  This row has been digitized, delete nonused stuff
1125   //
1126
1127   for(lp=0;lp<nofDigits;lp++){
1128     if(pList[lp]) delete [] pList[lp];
1129   }
1130   
1131   delete [] pList;
1132
1133   delete m1;
1134   delete m2;
1135   //  delete m3;
1136
1137 } // end of DigitizeRow
1138
1139 //_____________________________________________________________________________
1140
1141 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
1142                           Int_t *indexRange)
1143 {
1144
1145   //---------------------------------------------------------------
1146   //  Calculates 2-D signal (pad,time) for a single track,
1147   //  returns a pointer to the signal matrix and the track label 
1148   //  No digitization is performed at this level!!!
1149   //---------------------------------------------------------------
1150
1151   //-----------------------------------------------------------------
1152   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1153   // Modified: Marian Ivanov 
1154   //-----------------------------------------------------------------
1155
1156   TVector *tv;
1157  
1158   tv = (TVector*)p1->At(ntr); // pointer to a track
1159   TVector &v = *tv;
1160   
1161   Float_t label = v(0);
1162   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
1163
1164   Int_t nElectrons = (tv->GetNrows()-1)/4;
1165   indexRange[0]=9999; // min pad
1166   indexRange[1]=-1; // max pad
1167   indexRange[2]=9999; //min time
1168   indexRange[3]=-1; // max time
1169
1170   //  Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1171
1172   TMatrix &signal = *m1;
1173   TMatrix &total = *m2;
1174   //
1175   //  Loop over all electrons
1176   //
1177   for(Int_t nel=0; nel<nElectrons; nel++){
1178     Int_t idx=nel*4;
1179     Float_t aval =  v(idx+4);
1180     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac(); 
1181     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1182     Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
1183     
1184     if (n>0) for (Int_t i =0; i<n; i++){
1185        Int_t *index = fTPCParam->GetResBin(i);        
1186        Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1187        if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
1188          Int_t time=index[2];    
1189          Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
1190          weight *= eltoadcfac;
1191          
1192          if (m1!=0) signal(pad,time)+=weight; 
1193          total(pad,time)+=weight;
1194          indexRange[0]=TMath::Min(indexRange[0],pad);
1195          indexRange[1]=TMath::Max(indexRange[1],pad);
1196          indexRange[2]=TMath::Min(indexRange[2],time);
1197          indexRange[3]=TMath::Max(indexRange[3],time); 
1198        }         
1199     }
1200   } // end of loop over electrons
1201   
1202   return label; // returns track label when finished
1203 }
1204
1205 //_____________________________________________________________________________
1206 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
1207                      Float_t **pList)
1208 {
1209   //----------------------------------------------------------------------
1210   //  Updates the list of tracks contributing to digits for a given row
1211   //----------------------------------------------------------------------
1212
1213   //-----------------------------------------------------------------
1214   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1215   //-----------------------------------------------------------------
1216
1217   TMatrix &signal = *m;
1218
1219   // lop over nonzero digits
1220
1221   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1222     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1223
1224
1225         // accept only the contribution larger than 500 electrons (1/2 s_noise)
1226
1227         if(signal(ip,it)<0.5) continue; 
1228
1229
1230         Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1231         
1232         if(!pList[globalIndex]){
1233         
1234           // 
1235           // Create new list (6 elements - 3 signals and 3 labels),
1236           //
1237
1238           pList[globalIndex] = new Float_t [6];
1239
1240           // set list to -1 
1241
1242           *pList[globalIndex] = -1.;
1243           *(pList[globalIndex]+1) = -1.;
1244           *(pList[globalIndex]+2) = -1.;
1245           *(pList[globalIndex]+3) = -1.;
1246           *(pList[globalIndex]+4) = -1.;
1247           *(pList[globalIndex]+5) = -1.;
1248
1249
1250           *pList[globalIndex] = label;
1251           *(pList[globalIndex]+3) = signal(ip,it);
1252         }
1253         else{
1254
1255           // check the signal magnitude
1256
1257           Float_t highest = *(pList[globalIndex]+3);
1258           Float_t middle = *(pList[globalIndex]+4);
1259           Float_t lowest = *(pList[globalIndex]+5);
1260
1261           //
1262           //  compare the new signal with already existing list
1263           //
1264
1265           if(signal(ip,it)<lowest) continue; // neglect this track
1266
1267           //
1268
1269           if (signal(ip,it)>highest){
1270             *(pList[globalIndex]+5) = middle;
1271             *(pList[globalIndex]+4) = highest;
1272             *(pList[globalIndex]+3) = signal(ip,it);
1273
1274             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1275             *(pList[globalIndex]+1) = *pList[globalIndex];
1276             *pList[globalIndex] = label;
1277           }
1278           else if (signal(ip,it)>middle){
1279             *(pList[globalIndex]+5) = middle;
1280             *(pList[globalIndex]+4) = signal(ip,it);
1281
1282             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1283             *(pList[globalIndex]+1) = label;
1284           }
1285           else{
1286             *(pList[globalIndex]+5) = signal(ip,it);
1287             *(pList[globalIndex]+2) = label;
1288           }
1289         }
1290
1291     } // end of loop over pads
1292   } // end of loop over time bins
1293
1294
1295
1296 }//end of GetList
1297 //___________________________________________________________________
1298 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1299                         Stat_t ntracks,TObjArray **row)
1300 {
1301
1302   //-----------------------------------------------------------------
1303   // Prepares the sector digitization, creates the vectors of
1304   // tracks for each row of this sector. The track vector
1305   // contains the track label and the position of electrons.
1306   //-----------------------------------------------------------------
1307
1308   //-----------------------------------------------------------------
1309   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1310   //-----------------------------------------------------------------
1311
1312   Float_t gasgain = fTPCParam->GetGasGain();
1313   Int_t i;
1314   Float_t xyz[4]; 
1315
1316   AliTPChit *tpcHit; // pointer to a sigle TPC hit    
1317  
1318   //----------------------------------------------
1319   // Create TObjArray-s, one for each row,
1320   // each TObjArray will store the TVectors
1321   // of electrons, one TVector per each track.
1322   //---------------------------------------------- 
1323     
1324   for(i=0; i<nrows; i++){
1325     row[i] = new TObjArray;
1326   }
1327   Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
1328   TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
1329
1330   //--------------------------------------------------------------------
1331   //  Loop over tracks, the "track" contains the full history
1332   //--------------------------------------------------------------------
1333
1334   Int_t previousTrack,currentTrack;
1335   previousTrack = -1; // nothing to store so far!
1336
1337   for(Int_t track=0;track<ntracks;track++){
1338
1339     ResetHits();
1340
1341     TH->GetEvent(track); // get next track
1342     Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track
1343
1344     if(nhits == 0) continue; // no hits in the TPC for this track
1345
1346     //--------------------------------------------------------------
1347     //  Loop over hits
1348     //--------------------------------------------------------------
1349
1350     for(Int_t hit=0;hit<nhits;hit++){
1351
1352       tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit
1353       
1354       Int_t sector=tpcHit->fSector; // sector number
1355       if(sector != isec) continue; 
1356
1357         currentTrack = tpcHit->fTrack; // track number
1358         if(currentTrack != previousTrack){
1359                           
1360            // store already filled fTrack
1361               
1362            for(i=0;i<nrows;i++){
1363              if(previousTrack != -1){
1364                if(nofElectrons[i]>0){
1365                  TVector &v = *tracks[i];
1366                  v(0) = previousTrack;
1367                  tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1368                  row[i]->Add(tracks[i]);                     
1369                }
1370                else{
1371                  delete tracks[i]; // delete empty TVector
1372                  tracks[i]=0;
1373                }
1374              }
1375
1376              nofElectrons[i]=0;
1377              tracks[i] = new TVector(481); // TVectors for the next fTrack
1378
1379            } // end of loop over rows
1380                
1381            previousTrack=currentTrack; // update track label 
1382         }
1383            
1384         Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1385
1386        //---------------------------------------------------
1387        //  Calculate the electron attachment probability
1388        //---------------------------------------------------
1389
1390
1391         Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->fZ))
1392                                                         /fTPCParam->GetDriftV(); 
1393         // in microseconds!     
1394         Float_t attProb = fTPCParam->GetAttCoef()*
1395           fTPCParam->GetOxyCont()*time; //  fraction! 
1396    
1397         //-----------------------------------------------
1398         //  Loop over electrons
1399         //-----------------------------------------------
1400         Int_t index[3];
1401         index[1]=isec;
1402         for(Int_t nel=0;nel<qI;nel++){
1403           // skip if electron lost due to the attachment
1404           if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1405           xyz[0]=tpcHit->fX;
1406           xyz[1]=tpcHit->fY;
1407           xyz[2]=tpcHit->fZ;      
1408           xyz[3]= (Float_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1409           index[0]=1;
1410           
1411           TransportElectron(xyz,index); //MI change -august       
1412           Int_t rowNumber;
1413           fTPCParam->GetPadRow(xyz,index); //MI change august
1414           rowNumber = index[2];
1415           //transform position to local digit coordinates
1416           //relative to nearest pad row 
1417           if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue;       
1418           nofElectrons[rowNumber]++;      
1419           //----------------------------------
1420           // Expand vector if necessary
1421           //----------------------------------
1422           if(nofElectrons[rowNumber]>120){
1423             Int_t range = tracks[rowNumber]->GetNrows();
1424             if((nofElectrons[rowNumber])>(range-1)/4){
1425         
1426               tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
1427             }
1428           }
1429           
1430           TVector &v = *tracks[rowNumber];
1431           Int_t idx = 4*nofElectrons[rowNumber]-3;
1432
1433           v(idx)=  xyz[0];   // X - pad row coordinate
1434           v(idx+1)=xyz[1];   // Y - pad coordinate (along the pad-row)
1435           v(idx+2)=xyz[2];   // Z - time bin coordinate
1436           v(idx+3)=xyz[3];   // avalanche size  
1437         } // end of loop over electrons
1438         
1439       } // end of loop over hits
1440     } // end of loop over tracks
1441
1442     //
1443     //   store remaining track (the last one) if not empty
1444     //
1445
1446      for(i=0;i<nrows;i++){
1447        if(nofElectrons[i]>0){
1448           TVector &v = *tracks[i];
1449           v(0) = previousTrack;
1450           tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1451           row[i]->Add(tracks[i]);  
1452         }
1453         else{
1454           delete tracks[i];
1455           tracks[i]=0;
1456         }  
1457       }  
1458
1459           delete [] tracks;
1460           delete [] nofElectrons;
1461  
1462
1463 } // end of MakeSector
1464
1465
1466 //_____________________________________________________________________________
1467 void AliTPC::Init()
1468 {
1469   //
1470   // Initialise TPC detector after definition of geometry
1471   //
1472   Int_t i;
1473   //
1474   printf("\n");
1475   for(i=0;i<35;i++) printf("*");
1476   printf(" TPC_INIT ");
1477   for(i=0;i<35;i++) printf("*");
1478   printf("\n");
1479   //
1480   for(i=0;i<80;i++) printf("*");
1481   printf("\n");
1482 }
1483
1484 //_____________________________________________________________________________
1485 void AliTPC::MakeBranch(Option_t* option)
1486 {
1487   //
1488   // Create Tree branches for the TPC.
1489   //
1490   Int_t buffersize = 4000;
1491   char branchname[10];
1492   sprintf(branchname,"%s",GetName());
1493
1494   AliDetector::MakeBranch(option);
1495
1496   char *d = strstr(option,"D");
1497
1498   if (fDigits   && gAlice->TreeD() && d) {
1499     gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1500     printf("Making Branch %s for digits\n",branchname);
1501   }     
1502 }
1503  
1504 //_____________________________________________________________________________
1505 void AliTPC::ResetDigits()
1506 {
1507   //
1508   // Reset number of digits and the digits array for this detector
1509   //
1510   fNdigits   = 0;
1511   if (fDigits)   fDigits->Clear();
1512 }
1513
1514 //_____________________________________________________________________________
1515 void AliTPC::SetSecAL(Int_t sec)
1516 {
1517   //---------------------------------------------------
1518   // Activate/deactivate selection for lower sectors
1519   //---------------------------------------------------
1520
1521   //-----------------------------------------------------------------
1522   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1523   //-----------------------------------------------------------------
1524
1525   fSecAL = sec;
1526 }
1527
1528 //_____________________________________________________________________________
1529 void AliTPC::SetSecAU(Int_t sec)
1530 {
1531   //----------------------------------------------------
1532   // Activate/deactivate selection for upper sectors
1533   //---------------------------------------------------
1534
1535   //-----------------------------------------------------------------
1536   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1537   //-----------------------------------------------------------------
1538
1539   fSecAU = sec;
1540 }
1541
1542 //_____________________________________________________________________________
1543 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
1544 {
1545   //----------------------------------------
1546   // Select active lower sectors
1547   //----------------------------------------
1548
1549   //-----------------------------------------------------------------
1550   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1551   //-----------------------------------------------------------------
1552
1553   fSecLows[0] = s1;
1554   fSecLows[1] = s2;
1555   fSecLows[2] = s3;
1556   fSecLows[3] = s4;
1557   fSecLows[4] = s5;
1558   fSecLows[5] = s6;
1559 }
1560
1561 //_____________________________________________________________________________
1562 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
1563                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10, 
1564                        Int_t s11 , Int_t s12)
1565 {
1566   //--------------------------------
1567   // Select active upper sectors
1568   //--------------------------------
1569
1570   //-----------------------------------------------------------------
1571   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1572   //-----------------------------------------------------------------
1573
1574   fSecUps[0] = s1;
1575   fSecUps[1] = s2;
1576   fSecUps[2] = s3;
1577   fSecUps[3] = s4;
1578   fSecUps[4] = s5;
1579   fSecUps[5] = s6;
1580   fSecUps[6] = s7;
1581   fSecUps[7] = s8;
1582   fSecUps[8] = s9;
1583   fSecUps[9] = s10;
1584   fSecUps[10] = s11;
1585   fSecUps[11] = s12;
1586 }
1587
1588 //_____________________________________________________________________________
1589 void AliTPC::SetSens(Int_t sens)
1590 {
1591
1592   //-------------------------------------------------------------
1593   // Activates/deactivates the sensitive strips at the center of
1594   // the pad row -- this is for the space-point resolution calculations
1595   //-------------------------------------------------------------
1596
1597   //-----------------------------------------------------------------
1598   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1599   //-----------------------------------------------------------------
1600
1601   fSens = sens;
1602 }
1603  
1604 void AliTPC::SetSide(Float_t side=0.)
1605 {
1606   // choice of the TPC side
1607
1608   fSide = side;
1609  
1610 }
1611 //____________________________________________________________________________
1612 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
1613                            Float_t p2,Float_t p3)
1614 {
1615
1616   // gax mixture definition
1617
1618  fNoComp = nc;
1619  
1620  fMixtComp[0]=c1;
1621  fMixtComp[1]=c2;
1622  fMixtComp[2]=c3;
1623
1624  fMixtProp[0]=p1;
1625  fMixtProp[1]=p2;
1626  fMixtProp[2]=p3; 
1627  
1628  
1629 }
1630 //_____________________________________________________________________________
1631
1632 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
1633 {
1634   //
1635   // electron transport taking into account:
1636   // 1. diffusion, 
1637   // 2.ExB at the wires
1638   // 3. nonisochronity
1639   //
1640   // xyz and index must be already transformed to system 1
1641   //
1642
1643   fTPCParam->Transform1to2(xyz,index);
1644   
1645   //add diffusion
1646   Float_t driftl=xyz[2];
1647   if(driftl<0.01) driftl=0.01;
1648   driftl=TMath::Sqrt(driftl);
1649   Float_t sigT = driftl*(fTPCParam->GetDiffT());
1650   Float_t sigL = driftl*(fTPCParam->GetDiffL());
1651   xyz[0]=gRandom->Gaus(xyz[0],sigT);
1652   xyz[1]=gRandom->Gaus(xyz[1],sigT);
1653   xyz[2]=gRandom->Gaus(xyz[2],sigL);
1654
1655   // ExB
1656   
1657   if (fTPCParam->GetMWPCReadout()==kTRUE){
1658     Float_t x1=xyz[0];
1659     fTPCParam->Transform2to2NearestWire(xyz,index);
1660     Float_t dx=xyz[0]-x1;
1661     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
1662   }
1663   //add nonisochronity (not implemented yet)
1664   
1665 }
1666 //_____________________________________________________________________________
1667 void AliTPC::Streamer(TBuffer &R__b)
1668 {
1669   //
1670   // Stream an object of class AliTPC.
1671   //
1672    if (R__b.IsReading()) {
1673       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1674       AliDetector::Streamer(R__b);
1675       if (R__v < 2) return;
1676       R__b >> fNsectors;
1677    } else {
1678       R__b.WriteVersion(AliTPC::IsA());
1679       AliDetector::Streamer(R__b);
1680       R__b << fNsectors;
1681    }
1682 }
1683   
1684 ClassImp(AliTPCdigit)
1685  
1686 //_____________________________________________________________________________
1687 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
1688   AliDigit(tracks)
1689 {
1690   //
1691   // Creates a TPC digit object
1692   //
1693   fSector     = digits[0];
1694   fPadRow     = digits[1];
1695   fPad        = digits[2];
1696   fTime       = digits[3];
1697   fSignal     = digits[4];
1698 }
1699
1700  
1701 ClassImp(AliTPChit)
1702  
1703 //_____________________________________________________________________________
1704 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1705 AliHit(shunt,track)
1706 {
1707   //
1708   // Creates a TPC hit object
1709   //
1710   fSector     = vol[0];
1711   fPadRow     = vol[1];
1712   fX          = hits[0];
1713   fY          = hits[1];
1714   fZ          = hits[2];
1715   fQ          = hits[3];
1716 }
1717  
1718