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