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