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