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