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