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     AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField());
1805     if (field) {
1806       calib->SetExBField(field->SolenoidField());
1807     }
1808   }
1809
1810   Float_t gasgain = fTPCParam->GetGasGain();
1811   gasgain = gasgain/fGainFactor;
1812   Int_t i;
1813   Float_t xyz[5];
1814
1815   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1816   //MI change
1817   TBranch * branch=0;
1818   if (fHitType>1) branch = TH->GetBranch("TPC2");
1819   else branch = TH->GetBranch("TPC");
1820
1821
1822   //----------------------------------------------
1823   // Create TObjArray-s, one for each row,
1824   // each TObjArray will store the TVectors
1825   // of electrons, one TVectors per each track.
1826   //----------------------------------------------
1827
1828   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
1829   TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
1830
1831   for(i=0; i<nrows+2; i++){
1832     row[i] = new TObjArray;
1833     nofElectrons[i]=0;
1834     tracks[i]=0;
1835   }
1836
1837
1838
1839   //--------------------------------------------------------------------
1840   //  Loop over tracks, the "track" contains the full history
1841   //--------------------------------------------------------------------
1842
1843   Int_t previousTrack,currentTrack;
1844   previousTrack = -1; // nothing to store so far!
1845
1846   for(Int_t track=0;track<ntracks;track++){
1847     Bool_t isInSector=kTRUE;
1848     ResetHits();
1849     isInSector = TrackInVolume(isec,track);
1850     if (!isInSector) continue;
1851     //MI change
1852     branch->GetEntry(track); // get next track
1853
1854     //M.I. changes
1855
1856     tpcHit = (AliTPChit*)FirstHit(-1);
1857
1858     //--------------------------------------------------------------
1859     //  Loop over hits
1860     //--------------------------------------------------------------
1861
1862
1863     while(tpcHit){
1864
1865       Int_t sector=tpcHit->fSector; // sector number
1866       if(sector != isec){
1867         tpcHit = (AliTPChit*) NextHit();
1868         continue;
1869       }
1870
1871       // Remove hits which arrive before the TPC opening gate signal
1872       if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
1873           /fTPCParam->GetDriftV()+tpcHit->Time())<fTPCParam->GetGateDelay()) {
1874         tpcHit = (AliTPChit*) NextHit();
1875         continue;
1876       }
1877
1878       currentTrack = tpcHit->Track(); // track number
1879
1880       if(currentTrack != previousTrack){
1881
1882         // store already filled fTrack
1883
1884         for(i=0;i<nrows+2;i++){
1885           if(previousTrack != -1){
1886             if(nofElectrons[i]>0){
1887               TVector &v = *tracks[i];
1888               v(0) = previousTrack;
1889               tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
1890               row[i]->Add(tracks[i]);
1891             }
1892             else {
1893               delete tracks[i]; // delete empty TVector
1894               tracks[i]=0;
1895             }
1896           }
1897
1898           nofElectrons[i]=0;
1899           tracks[i] = new TVector(601); // TVectors for the next fTrack
1900
1901         } // end of loop over rows
1902
1903         previousTrack=currentTrack; // update track label
1904       }
1905
1906       Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1907
1908       //---------------------------------------------------
1909       //  Calculate the electron attachment probability
1910       //---------------------------------------------------
1911
1912
1913       Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
1914         /fTPCParam->GetDriftV();
1915       // in microseconds!
1916       Float_t attProb = fTPCParam->GetAttCoef()*
1917         fTPCParam->GetOxyCont()*time; //  fraction!
1918
1919       //-----------------------------------------------
1920       //  Loop over electrons
1921       //-----------------------------------------------
1922       Int_t index[3];
1923       index[1]=isec;
1924       for(Int_t nel=0;nel<qI;nel++){
1925         // skip if electron lost due to the attachment
1926         if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1927
1928         //
1929         // ExB effect
1930         //
1931         Double_t dxyz0[3],dxyz1[3];
1932         dxyz0[0]=tpcHit->X();
1933         dxyz0[1]=tpcHit->Y();
1934         dxyz0[2]=tpcHit->Z();
1935         if (calib->GetExB()){
1936           calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1937         }else{
1938           AliError("Not valid ExB calibration");
1939           dxyz1[0]=tpcHit->X();
1940           dxyz1[1]=tpcHit->Y();
1941           dxyz1[2]=tpcHit->Z();
1942         }
1943         xyz[0]=dxyz1[0];
1944         xyz[1]=dxyz1[1];
1945         xyz[2]=dxyz1[2];
1946         //
1947         //
1948         //
1949         // protection for the nonphysical avalanche size (10**6 maximum)
1950         //
1951         Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
1952         xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
1953         index[0]=1;
1954
1955         TransportElectron(xyz,index);
1956         Int_t rowNumber;
1957         Int_t padrow = fTPCParam->GetPadRow(xyz,index);
1958         //
1959         // Add Time0 correction due unisochronity
1960         // xyz[0] - pad row coordinate
1961         // xyz[1] - pad coordinate
1962         // xyz[2] - is in now time bin coordinate system
1963         Float_t correction =0;
1964         if (calib->GetPadTime0()){
1965           if (!calib->GetPadTime0()->GetCalROC(isec)) continue;
1966           Int_t npads = fTPCParam->GetNPads(isec,padrow);
1967           //      Int_t pad  = TMath::Nint(xyz[1]+fTPCParam->GetNPads(isec,TMath::Nint(xyz[0]))*0.5);
1968           // pad numbering from -npads/2 .. npads/2-1
1969           Int_t pad  = TMath::Nint(xyz[1]+npads/2);
1970           if (pad<0) pad=0;
1971           if (pad>=npads) pad=npads-1;
1972           correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(padrow,pad);
1973           //      printf("%d\t%d\t%d\t%f\n",isec,padrow,pad,correction);
1974           if (fDebugStreamer){
1975             (*fDebugStreamer)<<"Time0"<<
1976               "isec="<<isec<<
1977               "padrow="<<padrow<<
1978               "pad="<<pad<<
1979               "x0="<<xyz[0]<<
1980               "x1="<<xyz[1]<<
1981               "x2="<<xyz[2]<<
1982               "hit.="<<tpcHit<<
1983               "cor="<<correction<<
1984               "\n";
1985           }
1986         }
1987         xyz[2]+=correction;
1988         xyz[2]+=fTPCParam->GetNTBinsL1();    // adding Level 1 time bin offset
1989         //
1990         // Electron track time (for pileup simulation)
1991         xyz[2]+=tpcHit->Time()/fTPCParam->GetTSample(); // adding time of flight
1992         xyz[4] =0;
1993
1994         //
1995         // row 0 - cross talk from the innermost row
1996         // row fNRow+1 cross talk from the outermost row
1997         rowNumber = index[2]+1;
1998         //transform position to local digit coordinates
1999         //relative to nearest pad row
2000         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
2001         /*      Float_t x1,y1;
2002         if (isec <fTPCParam->GetNInnerSector()) {
2003           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
2004           y1 = fTPCParam->GetYInner(rowNumber);
2005         }
2006         else{
2007           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
2008           y1 = fTPCParam->GetYOuter(rowNumber);
2009         }
2010         // gain inefficiency at the wires edges - linear
2011         x1=TMath::Abs(x1);
2012         y1-=1.;
2013         if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.)); */
2014
2015         nofElectrons[rowNumber]++;
2016         //----------------------------------
2017         // Expand vector if necessary
2018         //----------------------------------
2019         if(nofElectrons[rowNumber]>120){
2020           Int_t range = tracks[rowNumber]->GetNrows();
2021           if((nofElectrons[rowNumber])>(range-1)/5){
2022
2023             tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons
2024           }
2025         }
2026
2027         TVector &v = *tracks[rowNumber];
2028         Int_t idx = 5*nofElectrons[rowNumber]-4;
2029         Real_t * position = &(((TVector&)v)(idx)); //make code faster
2030         memcpy(position,xyz,5*sizeof(Float_t));
2031
2032       } // end of loop over electrons
2033
2034       tpcHit = (AliTPChit*)NextHit();
2035
2036     } // end of loop over hits
2037   } // end of loop over tracks
2038
2039     //
2040     //   store remaining track (the last one) if not empty
2041     //
2042
2043   for(i=0;i<nrows+2;i++){
2044     if(nofElectrons[i]>0){
2045       TVector &v = *tracks[i];
2046       v(0) = previousTrack;
2047       tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2048       row[i]->Add(tracks[i]);
2049     }
2050     else{
2051       delete tracks[i];
2052       tracks[i]=0;
2053     }
2054   }
2055
2056   delete [] tracks;
2057   delete [] nofElectrons;
2058
2059 } // end of MakeSector
2060
2061
2062 //_____________________________________________________________________________
2063 void AliTPC::Init()
2064 {
2065   //
2066   // Initialise TPC detector after definition of geometry
2067   //
2068   AliDebug(1,"*********************************************");
2069 }
2070
2071 //_____________________________________________________________________________
2072 void AliTPC::ResetDigits()
2073 {
2074   //
2075   // Reset number of digits and the digits array for this detector
2076   //
2077   fNdigits   = 0;
2078   if (fDigits)   fDigits->Clear();
2079 }
2080
2081
2082
2083 //_____________________________________________________________________________
2084 void AliTPC::SetSens(Int_t sens)
2085 {
2086
2087   //-------------------------------------------------------------
2088   // Activates/deactivates the sensitive strips at the center of
2089   // the pad row -- this is for the space-point resolution calculations
2090   //-------------------------------------------------------------
2091
2092   //-----------------------------------------------------------------
2093   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2094   //-----------------------------------------------------------------
2095
2096   fSens = sens;
2097 }
2098
2099
2100 void AliTPC::SetSide(Float_t side=0.)
2101 {
2102   // choice of the TPC side
2103
2104   fSide = side;
2105
2106 }
2107 //_____________________________________________________________________________
2108
2109 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2110 {
2111   //
2112   // electron transport taking into account:
2113   // 1. diffusion,
2114   // 2.ExB at the wires
2115   // 3. nonisochronity
2116   //
2117   // xyz and index must be already transformed to system 1
2118   //
2119
2120   fTPCParam->Transform1to2(xyz,index);  // mis-alignment applied in this step
2121
2122   //add diffusion
2123   Float_t driftl=xyz[2];
2124   if(driftl<0.01) driftl=0.01;
2125   driftl=TMath::Sqrt(driftl);
2126   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2127   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2128   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2129   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2130   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2131
2132   // ExB
2133
2134   if (fTPCParam->GetMWPCReadout()==kTRUE){
2135     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2136     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2137   }
2138   //add nonisochronity (not implemented yet)
2139
2140
2141 }
2142
2143 ClassImp(AliTPChit)
2144   //______________________________________________________________________
2145   AliTPChit::AliTPChit()
2146             :AliHit(),
2147              fSector(0),
2148              fPadRow(0),
2149              fQ(0),
2150              fTime(0)
2151 {
2152   //
2153   // default
2154   //
2155
2156 }
2157 //_____________________________________________________________________________
2158 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
2159           :AliHit(shunt,track),
2160              fSector(0),
2161              fPadRow(0),
2162              fQ(0),
2163              fTime(0)
2164 {
2165   //
2166   // Creates a TPC hit object
2167   //
2168   fSector     = vol[0];
2169   fPadRow     = vol[1];
2170   fX          = hits[0];
2171   fY          = hits[1];
2172   fZ          = hits[2];
2173   fQ          = hits[3];
2174   fTime       = hits[4];
2175 }
2176
2177 //________________________________________________________________________
2178 // Additional code because of the AliTPCTrackHitsV2
2179
2180 void AliTPC::MakeBranch(Option_t *option)
2181 {
2182   //
2183   // Create a new branch in the current Root Tree
2184   // The branch of fHits is automatically split
2185   // MI change 14.09.2000
2186   AliDebug(1,"");
2187   if (fHitType<2) return;
2188   char branchname[10];
2189   sprintf(branchname,"%s2",GetName());
2190   //
2191   // Get the pointer to the header
2192   const char *cH = strstr(option,"H");
2193   //
2194   if (fTrackHits   && fLoader->TreeH() && cH && fHitType&4) {
2195     AliDebug(1,"Making branch for Type 4 Hits");
2196     fLoader->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2197   }
2198
2199 //   if (fTrackHitsOld   && fLoader->TreeH() && cH && fHitType&2) {
2200 //     AliDebug(1,"Making branch for Type 2 Hits");
2201 //     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2202 //                                                    fLoader->TreeH(),fBufferSize,99);
2203 //     fLoader->TreeH()->GetListOfBranches()->Add(branch);
2204 //   }
2205 }
2206
2207 void AliTPC::SetTreeAddress()
2208 {
2209   //Sets tree address for hits
2210   if (fHitType<=1) {
2211     if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2212     AliDetector::SetTreeAddress();
2213   }
2214   if (fHitType>1) SetTreeAddress2();
2215 }
2216
2217 void AliTPC::SetTreeAddress2()
2218 {
2219   //
2220   // Set branch address for the TrackHits Tree
2221   //
2222   AliDebug(1,"");
2223
2224   TBranch *branch;
2225   char branchname[20];
2226   sprintf(branchname,"%s2",GetName());
2227   //
2228   // Branch address for hit tree
2229   TTree *treeH = fLoader->TreeH();
2230   if ((treeH)&&(fHitType&4)) {
2231     branch = treeH->GetBranch(branchname);
2232     if (branch) {
2233       branch->SetAddress(&fTrackHits);
2234       AliDebug(1,"fHitType&4 Setting");
2235     }
2236     else
2237       AliDebug(1,"fHitType&4 Failed (can not find branch)");
2238
2239   }
2240  //  if ((treeH)&&(fHitType&2)) {
2241 //     branch = treeH->GetBranch(branchname);
2242 //     if (branch) {
2243 //       branch->SetAddress(&fTrackHitsOld);
2244 //       AliDebug(1,"fHitType&2 Setting");
2245 //     }
2246 //     else
2247 //       AliDebug(1,"fHitType&2 Failed (can not find branch)");
2248 //   }
2249 }
2250
2251 void AliTPC::FinishPrimary()
2252 {
2253   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();
2254   //  if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();
2255 }
2256
2257
2258 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2259 {
2260   //
2261   // add hit to the list
2262
2263   Int_t rtrack;
2264   if (fIshunt) {
2265     int primary = gAlice->GetMCApp()->GetPrimary(track);
2266     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2267     rtrack=primary;
2268   } else {
2269     rtrack=track;
2270     gAlice->GetMCApp()->FlagTrack(track);
2271   }
2272   if (fTrackHits && fHitType&4)
2273     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2274                              hits[1],hits[2],(Int_t)hits[3],hits[4]);
2275  //  if (fTrackHitsOld &&fHitType&2 )
2276 //     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2277 //                                 hits[1],hits[2],(Int_t)hits[3]);
2278
2279 }
2280
2281 void AliTPC::ResetHits()
2282 {
2283   if (fHitType&1) AliDetector::ResetHits();
2284   if (fHitType>1) ResetHits2();
2285 }
2286
2287 void AliTPC::ResetHits2()
2288 {
2289   //
2290   //reset hits
2291   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2292   // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2293
2294 }
2295
2296 AliHit* AliTPC::FirstHit(Int_t track)
2297 {
2298   if (fHitType>1) return FirstHit2(track);
2299   return AliDetector::FirstHit(track);
2300 }
2301 AliHit* AliTPC::NextHit()
2302 {
2303   //
2304   // gets next hit
2305   //
2306   if (fHitType>1) return NextHit2();
2307
2308   return AliDetector::NextHit();
2309 }
2310
2311 AliHit* AliTPC::FirstHit2(Int_t track)
2312 {
2313   //
2314   // Initialise the hit iterator
2315   // Return the address of the first hit for track
2316   // If track>=0 the track is read from disk
2317   // while if track<0 the first hit of the current
2318   // track is returned
2319   //
2320   if(track>=0) {
2321     gAlice->GetMCApp()->ResetHits();
2322     fLoader->TreeH()->GetEvent(track);
2323   }
2324   //
2325   if (fTrackHits && fHitType&4) {
2326     fTrackHits->First();
2327     return fTrackHits->GetHit();
2328   }
2329  //  if (fTrackHitsOld && fHitType&2) {
2330 //     fTrackHitsOld->First();
2331 //     return fTrackHitsOld->GetHit();
2332 //   }
2333
2334   else return 0;
2335 }
2336
2337 AliHit* AliTPC::NextHit2()
2338 {
2339   //
2340   //Return the next hit for the current track
2341
2342
2343 //   if (fTrackHitsOld && fHitType&2) {
2344 //     fTrackHitsOld->Next();
2345 //     return fTrackHitsOld->GetHit();
2346 //   }
2347   if (fTrackHits) {
2348     fTrackHits->Next();
2349     return fTrackHits->GetHit();
2350   }
2351   else
2352     return 0;
2353 }
2354
2355 void AliTPC::RemapTrackHitIDs(Int_t *map)
2356 {
2357   //
2358   // remapping
2359   //
2360   if (!fTrackHits) return;
2361
2362 //   if (fTrackHitsOld && fHitType&2){
2363 //     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2364 //     for (UInt_t i=0;i<arr->GetSize();i++){
2365 //       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2366 //       info->fTrackID = map[info->fTrackID];
2367 //     }
2368 //   }
2369 //  if (fTrackHitsOld && fHitType&4){
2370   if (fTrackHits && fHitType&4){
2371     TClonesArray * arr = fTrackHits->GetArray();;
2372     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2373       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2374       info->SetTrackID(map[info->GetTrackID()]);
2375     }
2376   }
2377 }
2378
2379 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2380 {
2381   //return bool information - is track in given volume
2382   //load only part of the track information
2383   //return true if current track is in volume
2384   //
2385   //  return kTRUE;
2386  //  if (fTrackHitsOld && fHitType&2) {
2387 //     TBranch * br = fLoader->TreeH()->GetBranch("fTrackHitsInfo");
2388 //     br->GetEvent(track);
2389 //     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2390 //     for (UInt_t j=0;j<ar->GetSize();j++){
2391 //       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2392 //     }
2393 //   }
2394
2395   if (fTrackHits && fHitType&4) {
2396     TBranch * br1 = fLoader->TreeH()->GetBranch("fVolumes");
2397     TBranch * br2 = fLoader->TreeH()->GetBranch("fNVolumes");
2398     br2->GetEvent(track);
2399     br1->GetEvent(track);
2400     Int_t *volumes = fTrackHits->GetVolumes();
2401     Int_t nvolumes = fTrackHits->GetNVolumes();
2402     if (!volumes && nvolumes>0) {
2403       AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes));
2404       return kFALSE;
2405     }
2406     for (Int_t j=0;j<nvolumes; j++)
2407       if (volumes[j]==id) return kTRUE;
2408   }
2409
2410   if (fHitType&1) {
2411     TBranch * br = fLoader->TreeH()->GetBranch("fSector");
2412     br->GetEvent(track);
2413     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2414       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2415     }
2416   }
2417   return kFALSE;
2418
2419 }
2420
2421
2422 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
2423 {
2424   //Makes TPC loader
2425   fLoader = new AliTPCLoader(GetName(),topfoldername);
2426   return fLoader;
2427 }
2428
2429 ////////////////////////////////////////////////////////////////////////
2430 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
2431 //
2432 // load TPC paarmeters from a given file or create new if the object
2433 // is not found there
2434 // 12/05/2003 This method should be moved to the AliTPCLoader
2435 // and one has to decide where to store the TPC parameters
2436 // M.Kowalski
2437   char paramName[50];
2438   sprintf(paramName,"75x40_100x60_150x60");
2439   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2440   if (paramTPC) {
2441     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
2442   } else {
2443     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
2444     //paramTPC = new AliTPCParamSR;
2445     paramTPC = AliTPCcalibDB::Instance()->GetParameters();
2446     if (!paramTPC->IsGeoRead()){
2447       //
2448       // read transformation matrices for gGeoManager
2449       //
2450       paramTPC->ReadGeoMatrices();
2451     }
2452
2453   }
2454   return paramTPC;
2455
2456 // the older version of parameters can be accessed with this code.
2457 // In some cases, we have old parameters saved in the file but
2458 // digits were created with new parameters, it can be distinguish
2459 // by the name of TPC TreeD. The code here is just for the case
2460 // we would need to compare with old data, uncomment it if needed.
2461 //
2462 //  char paramName[50];
2463 //  sprintf(paramName,"75x40_100x60");
2464 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2465 //  if (paramTPC) {
2466 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
2467 //  } else {
2468 //    sprintf(paramName,"75x40_100x60_150x60");
2469 //    paramTPC=(AliTPCParam*)in->Get(paramName);
2470 //    if (paramTPC) {
2471 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
2472 //    } else {
2473 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2474 //          <<endl;
2475 //      paramTPC = new AliTPCParamSR;
2476 //    }
2477 //  }
2478 //  return paramTPC;
2479
2480 }
2481
2482