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