TPC/AliTPCFast.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 //  fast TPC cluster simulation                                              //
  21 //                                                                          //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TParticle.h>
  25 #include <TVector.h>
  26 #include <TRandom.h>
  27
  28 #include "AliRunLoader.h"
  29 #include "AliRun.h"
  30 #include "AliMC.h"
  31 #include "AliTPC.h"
  32 #include "AliTPCParam.h"
  33 #include "AliTPCClustersArray.h"
  34 #include "AliTPCClustersRow.h"
  35 #include "AliTPCcluster.h"
  36 #include "AliComplexCluster.h"
  37 #include "AliTPCFast.h"
  38 #include "AliLog.h"
  39
  40 ClassImp(AliTPCFast)
  41   //____________________________________________________________________
  42 AliTPCFast::AliTPCFast(const AliTPCFast &param)
  43               :TObject(param),fParam(0)
  44
  45 {
  46   //
  47   //  copy constructor - dummy
  48   //
  49   fParam = param.fParam;
  50 }
  51 AliTPCFast & AliTPCFast::operator =(const AliTPCFast & param)
  52 {
  53   //
  54   // assignment operator - dummy
  55   //
  56   fParam=param.fParam;
  57   return (*this);
  58 }
  59
  60 //_____________________________________________________________________________
  61 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
  62 {
  63   //--------------------------------------------------------
  64   // TPC simple cluster generator from hits
  65   // obtained from the TPC Fast Simulator
  66   // The point errors are taken from the parametrization
  67   //--------------------------------------------------------
  68
  69   //-----------------------------------------------------------------
  70   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
  71   //-----------------------------------------------------------------
  72   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
  73   // Jouri.Belikov@cern.ch
  74   //----------------------------------------------------------------
  75
  76   /////////////////////////////////////////////////////////////////////////////
  77   //
  78   //---------------------------------------------------------------------
  79   //   ALICE TPC Cluster Parameters
  80   //--------------------------------------------------------------------
  81
  82
  83
  84   // Cluster width in rphi
  85   const Float_t kACrphi=0.18322;
  86   const Float_t kBCrphi=0.59551e-3;
  87   const Float_t kCCrphi=0.60952e-1;
  88   // Cluster width in z
  89   const Float_t kACz=0.19081;
  90   const Float_t kBCz=0.55938e-3;
  91   const Float_t kCCz=0.30428;
  92
  93
  94   AliLoader* loader = runLoader->GetLoader("TPCLoader");
  95   if (!loader) {
  96     AliError("No TPC loader found");
  97     return;
  98   }
  99   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 100   AliRun* aliRun = runLoader->GetAliRun();
 101   if (!aliRun) {
 102     AliError("Couldn't get AliRun object");
 103     return;
 104   }
 105   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 106   if (!tpc) {
 107     AliError("Couldn't get TPC detector");
 108     return;
 109   }
 110   AliTPCParam* param = tpc->GetParam();
 111   if (!param) {
 112     AliError("No TPC parameters available");
 113     return;
 114   }
 115
 116   //if(fDefaults == 0) SetDefaults();
 117
 118   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
 119   //
 120   TParticle *particle; // pointer to a given particle
 121   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 122   Int_t sector;
 123   Int_t ipart;
 124   Float_t xyz[5];
 125   Float_t pl,pt,tanth,rpad,ratio;
 126   Float_t cph,sph;
 127
 128   //---------------------------------------------------------------
 129   //  Get the access to the tracks
 130   //---------------------------------------------------------------
 131
 132   TTree *tH = loader->TreeH();
 133   if (tH == 0x0)
 134     AliFatal("Can not find TreeH in folder");
 135
 136   tpc->SetTreeAddress();
 137
 138   Stat_t ntracks = tH->GetEntries();
 139
 140   //Switch to the output file
 141
 142   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 143
 144   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
 145
 146   runLoader->CdGAFile();
 147   //param->Write(param->GetTitle());
 148
 149   AliTPCClustersArray carray;
 150   carray.Setup(param);
 151   carray.SetClusterType("AliTPCcluster");
 152   carray.MakeTree(loader->TreeR());
 153
 154   Int_t nclusters=0; //cluster counter
 155
 156   //------------------------------------------------------------
 157   // Loop over all sectors (72 sectors for 20 deg
 158   // segmentation for both lower and upper sectors)
 159   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
 160   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
 161   //
 162   // First cluster for sector 0 starts at "0"
 163   //------------------------------------------------------------
 164
 165   for(Int_t isec=0;isec<param->GetNSector();isec++){
 166     //MI change
 167     param->AdjustCosSin(isec,cph,sph);
 168
 169     //------------------------------------------------------------
 170     // Loop over tracks
 171     //------------------------------------------------------------
 172
 173     for(Int_t track=0;track<ntracks;track++){
 174       tpc->ResetHits();
 175       tH->GetEvent(track);
 176       //
 177       //  Get number of the TPC hits
 178       //
 179        tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 180
 181       // Loop over hits
 182       //
 183        while(tpcHit){
 184
 185          if (tpcHit->fQ == 0.) {
 186            tpcHit = (AliTPChit*) tpc->NextHit();
 187            continue; //information about track (I.Belikov)
 188          }
 189          sector=tpcHit->fSector; // sector number
 190
 191          if(sector != isec){
 192            tpcHit = (AliTPChit*) tpc->NextHit();
 193            continue;
 194          }
 195          ipart=tpcHit->Track();
 196          particle=aliRun->GetMCApp()->Particle(ipart);
 197          pl=particle->Pz();
 198          pt=particle->Pt();
 199          if(pt < 1.e-9) pt=1.e-9;
 200          tanth=pl/pt;
 201          tanth = TMath::Abs(tanth);
 202          rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
 203          ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 204
 205          //   space-point resolutions
 206
 207         sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
 208         sigmaZ   =AliTPCcluster::SigmaZ2(rpad,tanth   );
 209
 210         //   cluster widths
 211
 212         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
 213         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
 214
 215         // temporary protection
 216
 217         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
 218         if(sigmaZ < 0.) sigmaZ=0.4e-3;
 219         if(clRphi < 0.) clRphi=2.5e-3;
 220         if(clZ < 0.) clZ=2.5e-5;
 221
 222         //
 223
 224         //
 225         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
 226         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
 227         //
 228         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
 229         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
 230         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
 231           Float_t alpha=(isec < param->GetNInnerSector()) ?
 232           param->GetInnerAngle() : param->GetOuterAngle();
 233           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
 234           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
 235         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
 236           if (TMath::Abs(xyz[1])>param->GetZLength(isec)) xyz[1]=tpcHit->Z();
 237         xyz[2]=sigmaRphi;                                     // fSigmaY2
 238         xyz[3]=sigmaZ;                                        // fSigmaZ2
 239         xyz[4]=tpcHit->fQ;                                    // q
 240
 241         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
 242         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
 243
 244         Int_t tracks[3]={tpcHit->Track(), -1, -1};
 245         AliTPCcluster cluster(tracks,xyz);
 246
 247         clrow->InsertCluster(&cluster); nclusters++;
 248
 249         tpcHit = (AliTPChit*)tpc->NextHit();
 250
 251
 252       } // end of loop over hits
 253
 254     }   // end of loop over tracks
 255
 256     Int_t nrows=param->GetNRow(isec);
 257     for (Int_t irow=0; irow<nrows; irow++) {
 258         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
 259         if (!clrow) continue;
 260         carray.StoreRow(isec,irow);
 261         carray.ClearRow(isec,irow);
 262     }
 263
 264   } // end of loop over sectors
 265
 266   //  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
 267   loader->WriteRecPoints("OVERWRITE");
 268
 269
 270 } // end of function
 271
 272
 273
 274 //_________________________________________________________________
 275
 276
 277
 278 void  AliTPCFast::Hits2ExactClusters(AliRunLoader* runLoader) const{
 279
 280
 281
 282   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 283   if (!loader) {
 284     AliError("No TPC loader found");
 285     return;
 286   }
 287   AliTPC* tpc = (AliTPC*) gAlice->GetDetector("TPC");
 288   if (!tpc) {
 289     AliError("Couldn't get TPC detector");
 290     return;
 291   }
 292   AliTPCParam* param = tpc->GetParam();
 293   if (!param) {
 294     AliError("No TPC parameters available");
 295     return;
 296   }
 297   //---------------------------------------------------------------
 298   //  Get the access to the tracks
 299   //---------------------------------------------------------------
 300
 301   TTree *tH = loader->TreeH();
 302   if (tH == 0x0) { runLoader->LoadHits("TPC"); loader->TreeH();}
 303
 304   tpc->SetTreeAddress();
 305
 306
 307   //Switch to the output file
 308
 309   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 310
 311   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
 312
 313   runLoader->CdGAFile();
 314   //param->Write(param->GetTitle());
 315
 316   AliTPCClustersArray carray;
 317   carray.Setup(param);
 318   carray.SetClusterType("AliTPCclusterMI");
 319   carray.MakeTree(loader->TreeR());
 320
 321   //------------------------------------------------------------
 322   // Loop over all sectors (72 sectors for 20 deg)
 323   //------------------------------------------------------------
 324
 325   for(Int_t isec=0;isec<param->GetNSector();isec++){
 326     Hits2ExactClustersSector(runLoader, &carray, isec);
 327   } // end of loop over sectors
 328   loader->WriteRecPoints("OVERWRITE");
 329 }
 330
 331
 332
 333
 334 //_________________________________________________________________
 335 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
 336                                           AliTPCClustersArray* clustersArray,
 337                                           Int_t isec) const
 338 {
 339   //--------------------------------------------------------
 340   //calculate exact cross point of track and given pad row
 341   //resulting values are expressed in "local" coordinata
 342   //the sigmaY2 and sigmaZ2 of cluster are the shape of cluster parameters
 343   //                        - thay are used later on for error parameterization
 344   //--------------------------------------------------------
 345
 346   //-----------------------------------------------------------------
 347   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
 348   //-----------------------------------------------------------------
 349   //
 350   if (clustersArray==0){
 351     return;
 352   }
 353   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 354   if (!loader) {
 355     AliError("No TPC loader found");
 356     return;
 357   }
 358   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 359   AliRun* aliRun = runLoader->GetAliRun();
 360   if (!aliRun) {
 361     AliError("Couldn't get AliRun object");
 362     return;
 363   }
 364   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 365   if (!tpc) {
 366     AliError("Couldn't get TPC detector");
 367     return;
 368   }
 369   AliTPCParam* param = tpc->GetParam();
 370   if (!param) {
 371     AliError("No TPC parameters available");
 372     return;
 373   }
 374   //
 375   //
 376   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 377   //  Int_t sector,nhits;
 378   Int_t ipart;
 379   const Int_t kcmaxhits=30000;
 380   TVector * xxxx = new TVector(kcmaxhits*4);
 381   TVector & xxx = *xxxx;
 382   Int_t maxhits = kcmaxhits;
 383   //construct array for each padrow
 384   for (Int_t i=0; i<param->GetNRow(isec);i++)
 385     clustersArray->CreateRow(isec,i);
 386
 387   //---------------------------------------------------------------
 388   //  Get the access to the tracks
 389   //---------------------------------------------------------------
 390
 391   TTree *tH = loader->TreeH();
 392   if (tH == 0x0)
 393     AliFatal("Can not find TreeH in folder");
 394
 395   tpc->SetTreeAddress();
 396
 397   Stat_t ntracks = tH->GetEntries();
 398   //MI change
 399   TBranch * branch=0;
 400   if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
 401   else branch = tH->GetBranch("TPC");
 402
 403   //------------------------------------------------------------
 404   // Loop over tracks
 405   //------------------------------------------------------------
 406
 407   for(Int_t track=0;track<ntracks;track++){
 408     Bool_t isInSector=kTRUE;
 409     tpc->ResetHits();
 410     isInSector = tpc->TrackInVolume(isec,track);
 411     if (!isInSector) continue;
 412     //MI change
 413     branch->GetEntry(track); // get next track
 414     //
 415     // Loop over hits
 416     //
 417     Int_t currentIndex=0;
 418     Int_t lastrow=-1;  //last writen row
 419
 420     //M.I. changes
 421
 422     tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 423     while(tpcHit){
 424
 425       Int_t sector=tpcHit->fSector; // sector number
 426       if(sector != isec){
 427         tpcHit = (AliTPChit*) tpc->NextHit();
 428         continue;
 429       }
 430
 431       ipart=tpcHit->Track();
 432
 433       //find row number
 434
 435       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 436       Int_t    index[3]={1,isec,0};
 437       Int_t    currentrow = param->GetPadRow(x,index);
 438
 439       if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
 440       if (lastrow<0) lastrow=currentrow;
 441       //
 442       //
 443       if (currentrow!=lastrow){
 444         if (currentIndex>2){
 445           Float_t sumx=0;
 446           Float_t sumx2=0;
 447           Float_t sumx3=0;
 448           Float_t sumx4=0;
 449           Float_t sumy=0;
 450           Float_t sumxy=0;
 451           Float_t sumx2y=0;
 452           Float_t sumz=0;
 453           Float_t sumxz=0;
 454           Float_t sumx2z=0;
 455           Float_t sumq=0;
 456           for (Int_t cindex=0;cindex<currentIndex;cindex++){
 457             Float_t x1,x2,x3,x4;
 458             x1=x2=x3=x4=xxx(cindex*4);
 459             x2*=x1;
 460             x3*=x2;
 461             x4*=x3;
 462             sumx+=x1;
 463             sumx2+=x2;
 464             sumx3+=x3;
 465             sumx4+=x4;
 466             sumy+=xxx(cindex*4+1);
 467             sumxy+=xxx(cindex*4+1)*x1;
 468             sumx2y+=xxx(cindex*4+1)*x2;
 469             sumz+=xxx(cindex*4+2);
 470             sumxz+=xxx(cindex*4+2)*x1;
 471             sumx2z+=xxx(cindex*4+2)*x2;
 472             sumq+=xxx(cindex*4+3);
 473           }
 474           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 475             sumx2*(sumx*sumx3-sumx2*sumx2);
 476
 477           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 478             sumx2*(sumxy*sumx3-sumx2y*sumx2);
 479           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 480             sumx2*(sumxz*sumx3-sumx2z*sumx2);
 481
 482           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 483             sumx2*(sumx*sumx2y-sumx2*sumxy);
 484           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 485             sumx2*(sumx*sumx2z-sumx2*sumxz);
 486
 487           if (TMath::Abs(det)<0.00000000000000001){
 488              tpcHit = (AliTPChit*)tpc->NextHit();
 489             continue;
 490           }
 491
 492           Float_t y=detay/det;
 493           Float_t z=detaz/det;
 494           Float_t by=detby/det; //y angle
 495           Float_t bz=detbz/det; //z angle
 496           sumy/=Float_t(currentIndex);
 497           sumz/=Float_t(currentIndex);
 498
 499
 500           //
 501           //
 502           Float_t sign = (tpcHit->Z()>0)? 1.:-1.;
 503           y = (y+0.5)*param->GetPadPitchWidth(isec);
 504           z = z*param->GetZWidth();
 505           //
 506           // y expected shape
 507           Double_t sigmay2 = z*fParam->GetDiffL();
 508           sigmay2 *= sigmay2;
 509           Float_t angulary = by*param->GetPadPitchLength(isec,lastrow);
 510           angulary*=angulary;
 511           angulary/=12;
 512           sigmay2 +=angulary+0.25*param->GetPadPitchWidth(isec)*param->GetPadPitchWidth(isec);
 513           //
 514           // z expected shape
 515           Double_t sigmaz2 = z*fParam->GetDiffT();
 516           sigmaz2 *= sigmaz2;
 517           Float_t angularz = bz*param->GetPadPitchLength(isec,lastrow);
 518           angularz*=angularz;
 519           angularz/=12;
 520           sigmaz2 +=angularz+0.25;
 521           //
 522           sigmaz2 = TMath::Min(sigmaz2,1.);
 523           sigmay2 = TMath::Min(sigmay2,1.);
 524           //
 525           //
 526           z = sign*(param->GetZLength(isec) - z);
 527           if (TMath::Abs(z)< param->GetZLength(isec)-1){
 528             AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
 529             if (row!=0) {
 530               AliTPCclusterMI* cl = new((AliTPCclusterMI*)row->Append()) AliTPCclusterMI ;
 531               //          AliTPCclusterMI cl;
 532               cl->SetZ(z);
 533               cl->SetY(y);
 534               cl->SetQ(sumq);
 535               if (TMath::Abs(sumq)<1000){
 536                 cl->SetMax(Short_t(sumq));
 537               }
 538               else{
 539                 cl->SetMax(0);
 540               }
 541               cl->SetSigmaZ2(sigmaz2);
 542               cl->SetSigmaY2(sigmay2);
 543               cl->SetLabel(ipart,0);
 544               cl->SetLabel(-1,1);
 545               cl->SetLabel(-1,2);
 546               cl->SetType(0);
 547             }
 548           }
 549         } //end of calculating cluster for given row
 550         currentIndex=0;
 551         lastrow=currentrow;
 552       }  //  end of crossing rows
 553       //
 554       if ( currentIndex>=maxhits){
 555         maxhits+=kcmaxhits;
 556         xxx.ResizeTo(4*maxhits);
 557       }
 558       xxx(currentIndex*4)=x[0];
 559       xxx(currentIndex*4+1)=x[1];
 560       xxx(currentIndex*4+2)=x[2];
 561       xxx(currentIndex*4+3)=tpcHit->fQ;
 562       currentIndex++;
 563       tpcHit = (AliTPChit*)tpc->NextHit();
 564     } // end of loop over hits
 565   }   // end of loop over tracks
 566   //write padrows to tree
 567   for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
 568     clustersArray->StoreRow(isec,ii);
 569     clustersArray->ClearRow(isec,ii);
 570   }
 571   xxxx->Delete();
 572
 573 }
 574
 575