TPC/AliBarrelRec_TPCparam.C

   1 /****************************************************************************
   2  * This macro is supposed to do reconstruction in the ITS via Kalman        *
   3  * tracker V2. The ITStracker is feeded with parametrized TPC tracks        *
   4  *                                                                          *
   5  * It does the following steps:                                             *
   6  *             1) TPC tracking parameterization                             *
   7  *             2) ITS cluster finding V2 (via fast points !)                *
   8  *             3) ITS track finding V2                                      *
   9  *             4) Create a reference file with simulation info (p,PDG...)   *
  10  *                                                                          *
  11  * (Origin: A.Dainese, Padova, andrea.dainese@pd,infn.it                    *
  12  *  from AliBarrelReconstruction.C I.Belikov, CERN, Jouri.Belikov@cern.ch)  *
  13  ****************************************************************************/
  14
  15 #ifndef __CINT__
  16 #include <iostream.h>
  17 #include <TFile.h>
  18 #include <TStopwatch.h>
  19 #include <TObject.h>
  20 #include "alles.h"
  21 #include "AliRun.h"
  22 #include "AliHeader.h"
  23 #include "AliGenEventHeader.h"
  24 #include "AliMagF.h"
  25 #include "AliModule.h"
  26 #include "AliArrayI.h"
  27 #include "AliDigits.h"
  28 #include "AliITS.h"
  29 #include "AliTPC.h"
  30 #include "AliITSgeom.h"
  31 #include "AliITSRecPoint.h"
  32 #include "AliITSclusterV2.h"
  33 #include "AliITSsimulationFastPoints.h"
  34 #include "AliITStrackerV2.h"
  35 #include "AliKalmanTrack.h"
  36 #include "AliTPCtrackerParam.h"
  37 #endif
  38
  39 typedef struct {
  40   Int_t lab;
  41   Int_t pdg;
  42   Int_t mumlab;
  43   Int_t mumpdg;
  44   Float_t Vx,Vy,Vz;
  45   Float_t Px,Py,Pz;
  46 } RECTRACK;
  47
  48
  49 Int_t TPCParamTracks(const Char_t *galice,const Char_t *outname,const Int_t coll,const Double_t Bfield,Int_t n);
  50 Int_t ITSFindClusters(const Char_t *inname,const Char_t *outname,Int_t n);
  51 Int_t ITSFindTracks(const Char_t *galice,const Char_t *inname,const Char_t *inname2,const Char_t *outname,Int_t n);
  52 Int_t ITSMakeRefFile(const Char_t *galice, const Char_t *inname, const Char_t *outname, Int_t n);
  53
  54 Int_t AliBarrelRec_TPCparam(Int_t n=1) {
  55
  56   const Char_t *name=" AliBarrelRec_TPCparam";
  57   cerr<<'\n'<<name<<"...\n";
  58   gBenchmark->Start(name);
  59
  60
  61   const Char_t *TPCtrkNameS="AliTPCtracksParam.root";
  62   const Char_t *galiceName="galice.root";
  63   const Char_t *ITSclsName="AliITSclustersV2.root";
  64   const Char_t *ITStrkName="AliITStracksV2.root";
  65   const Char_t *ITSrefName="ITStracksRefFile.root";
  66
  67   // set here the code for the type of collision (needed for TPC tracking
  68   // parameterization). available collisions:
  69   //
  70   // coll = 0 ->   PbPb6000 (HIJING with b<2fm)
  71   const Int_t    collcode = 0;
  72   // set here the value of the magnetic field
  73   const Double_t BfieldValue = 0.4;
  74
  75
  76
  77   AliKalmanTrack::SetConvConst(100/0.299792458/BfieldValue);
  78
  79
  80   // ********** Build TPC tracks with parameterization *********** //
  81   if (TPCParamTracks(galiceName,TPCtrkNameS,collcode,BfieldValue,n)) {
  82     cerr<<"Failed to get TPC hits !\n";
  83     return 1;
  84   }
  85
  86
  87   // ********** Find ITS clusters *********** //
  88   if (ITSFindClusters(galiceName,ITSclsName,n)) {
  89     cerr<<"Failed to get ITS clusters !\n";
  90     return 1;
  91   }
  92
  93
  94   // ********* Find ITS tracks *********** //
  95   if (ITSFindTracks(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,n)) {
  96     cerr<<"Failed to get ITS tracks !\n";
  97     return 1;
  98   }
  99
 100
 101   // ********* Make ITS tracks reference file *********** //
 102   if (ITSMakeRefFile(galiceName,ITStrkName,ITSrefName,n)) {
 103     cerr<<"Failed to get ITS tracks ref file!\n";
 104     return 1;
 105   }
 106
 107   gBenchmark->Stop(name);
 108   gBenchmark->Show(name);
 109
 110   return 0;
 111 }
 112
 113
 114 Int_t TPCParamTracks(const Char_t *galice, const Char_t *outname,
 115                      const Int_t coll,const Double_t Bfield,Int_t n) {
 116
 117   cerr<<"\n*******************************************************************\n";
 118   Int_t rc;
 119
 120   const Char_t *name="TPCParamTracks";
 121   cerr<<'\n'<<name<<"...\n";
 122   gBenchmark->Start(name);
 123
 124   TFile *outfile=TFile::Open(outname,"recreate");
 125   TFile *infile =TFile::Open(galice);
 126
 127   AliTPCtrackerParam tracker(coll,Bfield);
 128   rc = tracker.BuildTPCtracks(infile,outfile,n);
 129
 130   delete gAlice; gAlice=0;
 131
 132   infile->Close();
 133   outfile->Close();
 134
 135   gBenchmark->Stop(name);
 136   gBenchmark->Show(name);
 137
 138   return rc;
 139 }
 140
 141 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
 142
 143
 144   cerr<<"\n*******************************************************************\n";
 145
 146   Int_t rc=0;
 147   const Char_t *name="ITSFindClusters";
 148   cerr<<'\n'<<name<<"...\n";
 149   gBenchmark->Start(name);
 150
 151
 152   // delete reconstruction Tree if it's there
 153   TFile *f =TFile::Open(inname,"update");
 154   f->Delete("TreeR0;*");
 155   f->Close();
 156
 157   TFile *out=TFile::Open(outname,"recreate");
 158   TFile *in =TFile::Open(inname,"update");
 159
 160   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
 161     cerr<<"Can't get gAlice !\n";
 162     return 1;
 163   }
 164
 165
 166   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 167   if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
 168   AliITSgeom *geom=ITS->GetITSgeom();
 169   out->cd();
 170   geom->Write();
 171
 172   Int_t ev=0;
 173   for (ev = 0; ev<n; ev++){
 174     in->cd();   // !!!! MI directory must point to galice. - othervise problem with Tree -connection
 175     gAlice->GetEvent(ev);
 176     //gAlice->TreeR()->Reset();   //reset reconstructed tree
 177
 178
 179     TTree *pTree=gAlice->TreeR();
 180     if (!pTree){
 181       gAlice->MakeTree("R");
 182       pTree = gAlice->TreeR();
 183     }
 184     TBranch *branch=pTree->GetBranch("ITSRecPoints");
 185     if (!branch) {
 186       //if not reconstructed ITS branch do reconstruction
 187       ITS->MakeBranch("R",0);
 188       //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
 189       for (Int_t i=0;i<3;i++) {
 190         ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
 191       }
 192       Int_t nsignal=25;
 193       Int_t size=-1;
 194       Int_t bgr_ev=Int_t(ev/nsignal);
 195       ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
 196       ///////////////////////////////////////////////////////////////////////
 197       gAlice->GetEvent(ev);   //MI comment  - in HitsToFast... they reset treeR to 0
 198       //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
 199       //time
 200     }
 201
 202
 203
 204     out->cd();
 205     TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
 206     char   cname[100];
 207     sprintf(cname,"TreeC_ITS_%d",ev);
 208
 209     TTree *cTree=new TTree(cname,"ITS clusters");
 210     cTree->Branch("Clusters",&clusters);
 211
 212     pTree=gAlice->TreeR();
 213     if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
 214     branch=pTree->GetBranch("ITSRecPoints");
 215     if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
 216     TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
 217     branch->SetAddress(&points);
 218
 219     TClonesArray &cl=*clusters;
 220     Int_t nclusters=0;
 221     Int_t nentr=(Int_t)pTree->GetEntries();
 222     AliITSgeom *geom=ITS->GetITSgeom();
 223
 224     for (Int_t i=0; i<nentr; i++) {
 225       if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
 226       Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
 227       Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
 228       Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
 229       Double_t yshift = x*rot[0] + y*rot[1];
 230       Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
 231       Int_t ncl=points->GetEntriesFast();
 232       nclusters+=ncl;
 233       Float_t lp[5];
 234       Int_t lab[6];
 235       for (Int_t j=0; j<ncl; j++) {
 236         AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
 237         lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
 238         lp[1]=p->GetZ()+zshift;
 239         lp[2]=p->GetSigmaX2();
 240         lp[3]=p->GetSigmaZ2();
 241         lp[4]=p->GetQ();
 242         lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
 243         lab[3]=ndet;
 244
 245         Int_t label=lab[0];
 246         TParticle *part=(TParticle*)gAlice->Particle(label);
 247         label=-3;
 248         while (part->P() < 0.005) {
 249           Int_t m=part->GetFirstMother();
 250           if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
 251           label=m;
 252           part=(TParticle*)gAlice->Particle(label);
 253         }
 254         if      (lab[1]<0) lab[1]=label;
 255         else if (lab[2]<0) lab[2]=label;
 256         else cerr<<"No empty labels !\n";
 257
 258         new(cl[j]) AliITSclusterV2(lab,lp);
 259       }
 260       cTree->Fill(); clusters->Delete();
 261       points->Delete();
 262     }
 263     cTree->Write();
 264     cerr<<"Number of clusters: "<<nclusters<<endl;
 265     delete cTree; delete clusters; delete points;
 266
 267   }
 268
 269
 270   delete gAlice; gAlice=0;
 271   in->Close();
 272   out->Close();
 273   gBenchmark->Stop(name);
 274   gBenchmark->Show(name);
 275
 276   return rc;
 277 }
 278
 279 Int_t ITSFindTracks(const Char_t *galice, const Char_t * inname,
 280                     const Char_t *inname2, const Char_t *outname,
 281                     Int_t n) {
 282
 283
 284   cerr<<"\n*******************************************************************\n";
 285
 286   Int_t rc=0;
 287   const Char_t *name="ITSFindTracks";
 288   cerr<<'\n'<<name<<"...\n";
 289   gBenchmark->Start(name);
 290
 291
 292   TFile *out=TFile::Open(outname,"recreate");
 293   TFile *in =TFile::Open(inname);
 294   TFile *in2 =TFile::Open(inname2);
 295
 296   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
 297   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
 298
 299
 300   for (Int_t ev=0; ev<n; ev++){
 301     AliITStrackerV2 tracker(geom,ev);
 302
 303     TArrayF o;
 304     o.Set(3);
 305     TFile* vtxFile = new TFile(galice);
 306     vtxFile->cd();
 307     AliHeader* header = 0;
 308     TTree* treeE = (TTree*)gDirectory->Get("TE");
 309     treeE->SetBranchAddress("Header",&header);
 310     treeE->GetEntry(ev);
 311     AliGenEventHeader* genHeader = header->GenEventHeader();
 312     if(genHeader) {
 313       // get primary vertex position
 314       genHeader->PrimaryVertex(o);
 315       vtxFile->Close();
 316       delete header;
 317     // set position of primary vertex
 318       Double_t vtx[3];
 319       vtx[0]=o[0]; vtx[1]=o[1]; vtx[2]=o[2];
 320       cerr<<"+++\n+++ Reading primary vertex position from galice.root\n+++\n+++ Setting primary vertex z = "<<vtx[2]<<" cm for ITS tracking\n+++\n";
 321       tracker.SetVertex(vtx);
 322     }else {
 323       cerr<<"+++\n+++ Event header not found in galice.root:\n+++ Primary vertex in (0,0,0) [default]\n+++\n";
 324     }
 325
 326     // setup vertex constraint in the two tracking passes
 327     Int_t flags[2];
 328     flags[0]=1;
 329     tracker.SetupFirstPass(flags);
 330     flags[0]=0;
 331     tracker.SetupSecondPass(flags);
 332
 333     rc=tracker.Clusters2Tracks(in,out);
 334
 335   }
 336
 337   delete gAlice;  gAlice=0;
 338   in->Close();
 339   in2->Close();
 340   out->Close();
 341
 342   gBenchmark->Stop(name);
 343   gBenchmark->Show(name);
 344
 345   return rc;
 346 }
 347
 348
 349 Int_t ITSMakeRefFile(const Char_t *galice, const Char_t *inname,
 350                      const Char_t *outname, Int_t n) {
 351
 352
 353   cerr<<"\n*******************************************************************\n";
 354
 355   Int_t rc=0;
 356   const Char_t *name="ITSMakeRefFile";
 357   cerr<<'\n'<<name<<"...\n";
 358   gBenchmark->Start(name);
 359
 360
 361   TFile *out = TFile::Open(outname,"recreate");
 362   TFile *trk = TFile::Open(inname);
 363   TFile *kin = TFile::Open(galice);
 364
 365
 366   // Get gAlice object from file
 367   if(!(gAlice=(AliRun*)kin->Get("gAlice"))) {
 368     cerr<<"gAlice has not been found on galice.root !\n";
 369     return 1;
 370   }
 371
 372
 373
 374   Int_t label;
 375   TParticle *Part;
 376   TParticle *Mum;
 377   static RECTRACK rectrk;
 378
 379
 380   for(Int_t event=0; event<n; event++){
 381
 382     AliITStrackV2 *itstrack=0;
 383
 384     gAlice->GetEvent(event);
 385
 386     trk->cd();
 387
 388     // Tree with ITS tracks
 389     char tname[100];
 390     sprintf(tname,"TreeT_ITS_%d",event);
 391
 392     TTree *tracktree=(TTree*)trk->Get(tname);
 393     TBranch *tbranch=tracktree->GetBranch("tracks");
 394     Int_t nentr=(Int_t)tracktree->GetEntries();
 395
 396     // Tree for true track parameters
 397     char ttname[100];
 398     sprintf(ttname,"Tree_Ref_%d",event);
 399     TTree *reftree = new TTree(ttname,"Tree with true track params");
 400     reftree->Branch("rectracks",&rectrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
 401
 402     for (Int_t i=0; i<nentr; i++) {
 403       itstrack=new AliITStrackV2;
 404       tbranch->SetAddress(&itstrack);
 405       tracktree->GetEvent(i);
 406       label = TMath::Abs(itstrack->GetLabel());
 407
 408       Part = (TParticle*)gAlice->Particle(label);
 409       rectrk.lab=label;
 410       rectrk.pdg=Part->GetPdgCode();
 411       rectrk.mumlab = Part->GetFirstMother();
 412       if(Part->GetFirstMother()>=0) {
 413         Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother());
 414         rectrk.mumpdg=Mum->GetPdgCode();
 415       } else {
 416         rectrk.mumpdg=-1;
 417       }
 418       rectrk.Vx=Part->Vx();
 419       rectrk.Vy=Part->Vy();
 420       rectrk.Vz=Part->Vz();
 421       rectrk.Px=Part->Px();
 422       rectrk.Py=Part->Py();
 423       rectrk.Pz=Part->Pz();
 424
 425       reftree->Fill();
 426     }
 427
 428     out->cd();
 429     reftree->Write();
 430
 431   }
 432
 433   trk->Close();
 434   kin->Close();
 435   out->Close();
 436
 437   gBenchmark->Stop(name);
 438   gBenchmark->Show(name);
 439
 440
 441   return rc;
 442
 443 }
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455