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