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