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