TPC/AliBarrelReconstruction.C

   1 /****************************************************************************
   2  * This macro is supposed to do reconstruction in the barrel ALICE trackers *
   3  * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
   4  * It does the following steps (April 12, 2001):                            *
   5  *                   1) TPC cluster finding                                 *
   6  *                   2) TPC track finding                                   *
   7  *                   3) ITS cluster finding V2 (via fast points !)          *
   8  *                   4) ITS track finding V2                                *
   9  *                   5) ITS back track propagation V2                       *
  10  *                   6) TPC back track propagation                          *
  11  *                (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch)          *
  12  ****************************************************************************/
  13
  14 #ifndef __CINT__
  15   #include "alles.h"
  16   #include "iostream.h"
  17   #include "AliRun.h"
  18   #include "AliRunLoader.h"
  19   #include "AliLoader.h"
  20   #include "AliTPCLoader.h"
  21   #include "AliITSLoader.h"
  22   #include "AliMagF.h"
  23   #include "AliTPCtracker.h"
  24   #include "AliITS.h"
  25   #include "AliITSgeom.h"
  26   #include "AliITSRecPoint.h"
  27   #include "AliITSclusterV2.h"
  28   #include "AliITSsimulationFastPoints.h"
  29   #include "AliITStrackerV2.h"
  30 #endif
  31
  32
  33
  34 Int_t TPCFindClusters( Int_t n);
  35 Int_t TPCFindTracks(Int_t n);
  36 Int_t TPCSortTracks(const Char_t *outname,  Int_t n);
  37 Int_t TPCPropagateBack();
  38
  39 Int_t ITSFindClusters(Int_t n);
  40 Int_t ITSFindTracks(const Char_t *inname2, Int_t n);
  41 Int_t ITSPropagateBack();
  42
  43 const char* TPCtrkNameS= "TPC.TracksSorted.root";
  44
  45 class AliRunLoader;
  46 class AliTPCLoader;
  47 class AliTPCParam;
  48 AliRunLoader *rl = 0x0;
  49 AliTPCLoader *tpcl = 0x0;
  50 AliTPCParam  *param = 0x0;
  51 Bool_t debug = kFALSE;
  52 Int_t AliBarrelReconstruction(Int_t n=1)
  53  {
  54
  55    if (gAlice)
  56     {
  57       delete gAlice->GetRunLoader();
  58       delete gAlice;//if everything was OK here it is already NULL
  59       gAlice = 0x0;
  60     }
  61
  62    rl = AliRunLoader::Open("galice.root");
  63    if (rl == 0x0)
  64     {
  65       cerr<<"Can not open session"<<endl;
  66       return 1;
  67     }
  68
  69    if (rl->LoadgAlice())
  70     {
  71       cerr<<"Error occured while l"<<endl;
  72       return 1;
  73     }
  74    AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
  75
  76    tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
  77    if (tpcl == 0x0)
  78     {
  79       cerr<<"Can not get TPC Loader"<<endl;
  80       return 1;
  81     }
  82
  83 // ********** Find TPC clusters *********** //
  84    if ( TPCFindClusters(n) )
  85     {
  86       cerr<<"Failed to get TPC clusters !\n";
  87       return 1;
  88     }
  89
  90 // ********** Find TPC tracks *********** //
  91     if (TPCFindTracks(n))
  92      {
  93       cerr<<"Failed to get TPC tracks !\n";
  94       return 1;
  95      }
  96
  97 //   cout<<"Stopping tracking on TPC\n";
  98 // ********** Sort and label TPC tracks *********** //
  99    if (TPCSortTracks(TPCtrkNameS,n)) {
 100       cerr<<"Failed to sort TPC tracks !\n";
 101       return 1;
 102     }
 103
 104 //   cout<<"Stopping tracking on TPC sorted T\n";
 105
 106
 107 // ********** Find ITS clusters *********** //
 108    if (ITSFindClusters(n))
 109     {
 110       cerr<<"Failed to get ITS clusters !\n";
 111       return 1;
 112     }
 113
 114 // ********** Find ITS tracks *********** //
 115
 116    if (ITSFindTracks(TPCtrkNameS,n))
 117     {
 118       cerr<<"Failed to get ITS tracks !\n";
 119       return 1;
 120     }
 121    //clsFile->Close();
 122    cout<<"Stopping on ITSFindTracks\n";
 123    delete rl;
 124    rl = 0x0;
 125    return 0;
 126    return 1;
 127
 128 // ********** Back propagation of the ITS tracks *********** //
 129    if ( ITSPropagateBack() ) {
 130       cerr<<"Failed to propagate back the ITS tracks !\n";
 131       return 1;
 132    }
 133
 134
 135 // ********** Back propagation of the TPC tracks *********** //
 136
 137    if (TPCPropagateBack()) {
 138       cerr<<"Failed to propagate back the TPC tracks !\n";
 139       return 1;
 140    }
 141
 142    return 0;
 143 }
 144
 145
 146 Int_t TPCFindClusters(Int_t n)
 147  {
 148    Int_t rc=0;
 149    const Char_t *name="TPCFindClusters";
 150    cerr<<'\n'<<name<<"...\n";
 151    rl->CdGAFile();
 152    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
 153    if (!param)
 154     {
 155      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
 156      if (!param)
 157       {
 158         cerr<<"TPC parameters have not been found !\n";
 159         return 1;
 160       }
 161     }
 162 //   param->Dump();
 163    param->Update();
 164
 165    gBenchmark->Start(name);
 166
 167    AliTPCv2 tpc;
 168    tpc.SetParam(param);
 169    tpc.SetLoader(tpcl);
 170
 171    tpcl->LoadDigits("read");
 172    tpcl->LoadRecPoints("recreate");
 173
 174
 175    //tpc.Digits2Clusters(out); //MI change
 176    for (Int_t i=0;i<n;i++)
 177     {
 178       printf("Processing event %d\n",i);
 179       tpc.Digits2Clusters(i);
 180      //     AliTPCclusterer::Digits2Clusters(dig, out, i);
 181     }
 182
 183    tpcl->UnloadDigits();
 184    tpcl->UnloadRecPoints();
 185    gBenchmark->Stop(name);
 186    gBenchmark->Show(name);
 187
 188    return rc;
 189 }
 190
 191 Int_t TPCFindTracks(Int_t n) {
 192
 193    Int_t rc=0;
 194    AliTPCtracker *tracker = 0x0;
 195    const Char_t *name="TPCFindTracks";
 196    rl->GetEvent(0);
 197    rl->CdGAFile();
 198    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
 199
 200    if (!param)
 201     {
 202      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
 203      if (!param)
 204       {
 205         cerr<<"TPC parameters have not been found !\n";
 206         return 1;
 207       }
 208     }
 209 //   param->Dump();
 210    param->Update();
 211
 212    gBenchmark->Start(name);
 213
 214    tpcl->LoadRecPoints("read");
 215    tpcl->LoadTracks("recreate");
 216
 217    //AliTPCtracker *tracker=new AliTPCtracker(param);
 218    //rc=tracker->Clusters2Tracks(0,out);
 219    //delete tracker;
 220
 221    for (Int_t i=0;i<n;i++){
 222      rl->GetEvent(i);
 223      printf("Processing event %d\n",i);
 224
 225      tracker = new AliTPCtracker(param, i, AliConfig::fgkDefaultEventFolderName.Data());
 226      //Int_t rc=
 227      tracker->Clusters2Tracks();
 228      delete tracker;
 229    }
 230    tpcl->UnloadRecPoints();
 231    tpcl->UnloadTracks();
 232
 233    gBenchmark->Stop(name);
 234    gBenchmark->Show(name);
 235
 236    return rc;
 237 }
 238
 239
 240 Int_t TPCSortTracks(const Char_t * outname,  Int_t eventn){
 241    Int_t rc=0;
 242    const Char_t *name="TPCSortTracks";
 243
 244    cerr<<'\n'<<name<<"...\n";
 245    rl->CdGAFile();
 246    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
 247    if (!param)
 248     {
 249      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
 250      if (!param)
 251       {
 252         cerr<<"TPC parameters have not been found !\n";
 253         return 1;
 254       }
 255     }
 256    param->Update();
 257
 258    gBenchmark->Start(name);
 259
 260
 261    AliRunLoader* rl2 = AliRunLoader::Open("galice.root","tmp");
 262
 263    AliLoader* tpcl2 = (AliTPCLoader*)rl2->GetLoader("TPCLoader");
 264    cout<<"tpcl2->SetTracksFileName("<<outname<<");\n";
 265    tpcl2->SetTracksFileName(TString(outname));
 266    tpcl2->LoadTracks("recreate");
 267
 268    // loop over events
 269    for (Int_t event=0;event<eventn; event++)
 270     {
 271
 272      rl->GetEvent(event);
 273      rl2->GetEvent(event);
 274
 275      TObjArray tarray(10000);
 276      AliTPCtrack *iotrack=0;
 277      Int_t i;
 278
 279      if (tpcl->TreeT() == 0x0)   tpcl->LoadTracks("read");
 280      TTree *tracktree=tpcl->TreeT();
 281      if (tracktree == 0x0)
 282       {
 283         cerr<<"Can not get TreeT for event "<<event<<endl;
 284         continue;
 285       }
 286      TBranch *tbranch=tracktree->GetBranch("tracks");
 287      Int_t nentr=(Int_t)tracktree->GetEntries();
 288      cout<<"IN Tracks nentr = "<<nentr<<endl;
 289      for (i=0; i<nentr; i++)
 290       {
 291        iotrack=new AliTPCtrack;
 292        tbranch->SetAddress(&iotrack);
 293        tracktree->GetEvent(i);
 294        tarray.AddLast(iotrack);
 295       }
 296      tarray.Sort();
 297
 298      //assign thacks GEANT labels
 299      cout<<"Running Tracker\n";
 300      AliTPCtracker *tracker = new AliTPCtracker(param, event);
 301
 302      cout<<"Load Sectors\n";
 303
 304      tracker->LoadInnerSectors();
 305      tracker->LoadOuterSectors();
 306
 307      cout<<"Cooking Labels\n";
 308
 309      for (i=0; i<nentr; i++)
 310       {
 311        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
 312        tracker->CookLabel(iotrack,0.1);
 313       }
 314      cout<<"deleting tracker\n";
 315      delete tracker;
 316      //in->Close();
 317      //end of GEANT label assignment
 318
 319
 320      if (tpcl2->TreeT() == 0x0)   tpcl2->MakeTree("T");
 321      tracktree= tpcl2->TreeT();
 322      if (tracktree == 0x0)
 323       {
 324         cerr<<"Can not get TreeT for Sorted tracks\n";
 325         return 1;
 326       }
 327      tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 328      for (i=0; i<nentr; i++) {
 329        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
 330        tracktree->Fill();
 331      }
 332
 333      cout<<" Tracks File Name is "<<tpcl2->GetTracksFileName()<<endl;
 334      tpcl2->WriteTracks("OVERWRITE");
 335    }
 336
 337
 338    delete rl2;
 339
 340    tpcl->UnloadTracks();
 341    tpcl->UnloadRecPoints();
 342
 343    gBenchmark->Stop(name);
 344    gBenchmark->Show(name);
 345
 346    return rc;
 347 }
 348
 349 Int_t ITSFindClusters(Int_t n)
 350  {
 351    Int_t rc=0;
 352    const Char_t *name="ITSFindClusters";
 353    cerr<<'\n'<<name<<"...\n";
 354    gBenchmark->Start(name);
 355
 356    if (rl->GetEvent(0))
 357     {
 358       cerr<<"Problems\n";
 359       return 1;
 360     }
 361
 362    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
 363    if (itsl == 0x0)
 364     {
 365       cerr<<"Can not get ITS Loader from Run Loader\n";
 366       return 1;
 367     }
 368
 369    gAlice=rl->GetAliRun();
 370    if (!gAlice) {
 371       cerr<<"Can't get gAlice !\n";
 372       return 1;
 373    }
 374    if (rl->TreeK() == 0x0) rl->LoadKinematics();
 375    if (rl->TreeE() == 0x0) rl->LoadHeader();
 376
 377    itsl->LoadRawClusters("recreate");
 378
 379    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 380    if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
 381    AliITSgeom *geom=ITS->GetITSgeom();
 382
 383    Int_t ev=0;
 384    for (ev = 0; ev<n; ev++)
 385    {
 386
 387      rl->GetEvent(ev);
 388      TBranch *branch = 0x0;
 389      TTree *pTree = 0x0;
 390
 391      pTree=itsl->TreeR();
 392      if (pTree == 0x0)
 393       {
 394         itsl->LoadRecPoints("read");
 395         pTree=itsl->TreeR();
 396         if (pTree == 0x0)
 397          {
 398            ::Error("AliBarrelReonstruction.C::ITSFindClusters",
 399                    "Can not get TreeR for event %d",ev);
 400            return 1;
 401          }
 402       }
 403
 404      pTree=itsl->TreeR();
 405      branch=(pTree)?pTree->GetBranch("ITSRecPoints"):0x0;
 406
 407      if (branch== 0x0) {
 408        //if not reconstructed ITS branch do reconstruction
 409        if (debug) ::Info("AliBarrelReconstruction.C","Did not get ITSRecPoints from TreeR.");
 410        if (debug) ::Info("AliBarrelReconstruction.C","Making branch and Producing RecPoints");
 411        itsl->SetRecPointsFileOption("recreate");
 412        if (itsl->TreeR()==0x0) itsl->MakeTree("R");
 413
 414        ITS->MakeBranch("RF");
 415        ITS->SetTreeAddress();
 416        itsl->LoadHits();
 417
 418        //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
 419        for (Int_t i=0;i<3;i++) {
 420          ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
 421        }
 422        Int_t nsignal=25;
 423        Int_t size=-1;
 424        Int_t bgr_ev=Int_t(ev/nsignal);
 425
 426        ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
 427        //////////////////////////////////////////////////////////////////////////
 428        //MI comment  - in HitsToFast... they reset treeR to 0
 429        //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
 430        //time
 431        itsl->UnloadHits();
 432        ITS->SetTreeAddress();
 433      }
 434
 435      TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
 436
 437      if (itsl->TreeC() == 0x0) itsl->MakeTree("C");
 438      TTree *cTree = itsl->TreeC();
 439      TBranch* b = cTree->GetBranch("Clusters");
 440      if (b == 0x0)
 441        cTree->Branch("Clusters",&clusters);
 442      else b->SetAddress(&clusters);
 443
 444
 445      if (itsl->TreeR() == 0x0) itsl->LoadRecPoints();
 446
 447      if (branch == 0x0)
 448       {//it means that we produced FastRecPoints above
 449        pTree=itsl->TreeR();
 450        if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
 451        branch=pTree->GetBranch("ITSRecPointsF");
 452        if (!branch) { cerr<<"Can't get Fast RecPoints branch named ITSRecPointsF  !\n"; return 1;}
 453       }
 454      TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
 455      branch->SetAddress(&points);
 456
 457      TClonesArray &cl=*clusters;
 458      Int_t nclusters=0;
 459      Int_t nentr=(Int_t)branch->GetEntries();
 460      AliITSgeom *geom=ITS->GetITSgeom();
 461
 462      cout<<"\n\n Number of entries in TreeR = "<<nentr<<endl;
 463      for (Int_t i=0; i<nentr; i++)
 464       {
 465        if (!branch->GetEvent(i))
 466         {
 467           cTree->Fill();
 468           continue;
 469         }
 470
 471        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
 472        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
 473        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
 474        Double_t yshift = x*rot[0] + y*rot[1];
 475        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
 476        Int_t ncl=points->GetEntriesFast();
 477        nclusters+=ncl;
 478
 479        if (debug) ::Info("AliBarrelReconstruction.C",
 480               "i=%d lay=%d lad=%d det=%d NRP=%d",
 481                i,   lay,   lad,   det,   ncl);
 482
 483        for (Int_t j=0; j<ncl; j++)
 484         {
 485           Float_t lp[5];
 486           Int_t lab[4];
 487
 488           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
 489           lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
 490           lp[1]=p->GetZ()+zshift;
 491           lp[2]=p->GetSigmaX2();
 492           lp[3]=p->GetSigmaZ2();
 493           lp[4]=p->GetQ();
 494           lab[0]=p->GetLabel(0);
 495           lab[1]=p->GetLabel(1);
 496           lab[2]=p->GetLabel(2);
 497           lab[3]=ndet;
 498
 499           Int_t label=lab[0];
 500           if (label<0) continue;
 501           TParticle *part=(TParticle*)rl->Stack()->Particle(label);
 502           label=-3;
 503           while (part->P() < 0.005) {
 504             Int_t m=part->GetFirstMother();
 505             if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
 506             label=m;
 507             part=(TParticle*)gAlice->Particle(label);
 508           }
 509           if (lab[1]<0) lab[1]=label;
 510           else if (lab[2]<0) lab[2]=label;
 511                else cerr<<"No empty labels !\n";
 512
 513           new(cl[j]) AliITSclusterV2(lab,lp);
 514          }
 515        cTree->Fill(); clusters->Delete();
 516        points->Delete();
 517      }
 518      itsl->WriteRawClusters("OVERWRITE");
 519
 520      cerr<<"Number of clusters: "<<nclusters<<endl;
 521      delete clusters; delete points;
 522
 523    }
 524    itsl->UnloadRecPoints();
 525    itsl->UnloadRawClusters();
 526    rl->UnloadKinematics();
 527
 528    gBenchmark->Stop(name);
 529    gBenchmark->Show(name);
 530
 531    return rc;
 532 }
 533
 534 Int_t ITSFindTracks(const Char_t *inname2, Int_t n)
 535 {
 536
 537    Int_t rc=0;
 538    const Char_t *name="ITSFindTracks";
 539    cerr<<'\n'<<name<<"...\n";
 540    gBenchmark->Start(name);
 541
 542    rl->GetEvent(0);
 543
 544    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
 545    if ( (itsl == 0x0) || (tpcl == 0x0))
 546     {
 547       cerr<<"Can not get ITS Loader from Run Loader\n";
 548       return 1;
 549     }
 550
 551    tpcl->SetTracksFileName(TString(inname2));
 552    if(rl->GetAliRun() == 0x0) rl->LoadgAlice();
 553
 554    AliITS* ITS = rl->GetAliRun()->GetModule("ITS");
 555
 556    rl->CdGAFile();
 557    AliITSgeom *geom=ITS->GetITSgeom();
 558    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
 559
 560    itsl->LoadTracks("recreate");
 561    itsl->LoadRawClusters("read");
 562    for (Int_t i=0;i<n;i++)
 563     {
 564       rl->GetEvent(i);
 565       AliITStrackerV2 tracker(geom,i,AliConfig::fgkDefaultEventFolderName);
 566       rc=tracker.Clusters2Tracks();
 567     }
 568
 569
 570    gBenchmark->Stop(name);
 571    gBenchmark->Show(name);
 572
 573    return rc;
 574 }
 575
 576 Int_t ITSPropagateBack()
 577 {
 578
 579   Int_t rc=0;
 580   /*
 581    const Char_t *name="ITSPropagateBack";
 582    cerr<<'\n'<<name<<"...\n";
 583    gBenchmark->Start(name);
 584
 585    AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
 586    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
 587    AliITStrackerV2 tracker(geom);
 588
 589    TFile *out=TFile::Open(outname,"update");
 590    TFile *in =TFile::Open(inname);
 591    rc=tracker.PropagateBack(in,out);
 592    in->Close();
 593    out->Close();
 594
 595    gBenchmark->Stop(name);
 596    gBenchmark->Show(name);
 597   */
 598    return rc;
 599 }
 600
 601 Int_t TPCPropagateBack()
 602 {
 603    Int_t rc=0;
 604    const Char_t *name="TPCPropagateBack";
 605    cerr<<'\n'<<name<<"...\n";
 606    gBenchmark->Start(name);
 607
 608    param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
 609    if (!param)
 610     {
 611      param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
 612      if (!param)
 613       {
 614         cerr<<"TPC parameters have not been found !\n";
 615         return 1;
 616       }
 617     }
 618 //   param->Dump();
 619    param->Update();
 620
 621    AliTPCtracker *tracker = new AliTPCtracker(param, AliConfig::fgkDefaultEventFolderName);
 622
 623 //   TFile *out=TFile::Open(outname,"update");
 624 //   TFile *in =TFile::Open(inname);
 625
 626    rc=tracker->PropagateBack();
 627    delete tracker;
 628
 629    gBenchmark->Stop(name);
 630    gBenchmark->Show(name);
 631
 632    return rc;
 633 }