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