TPC/AliTPC.cxx

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