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