TPC/AliTPC.cxx

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