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 "AliMagF.h"
  17   #include "AliTPCtracker.h"
  18
  19   #include "AliITS.h"
  20   #include "AliITSgeom.h"
  21   #include "AliITSRecPoint.h"
  22   #include "AliITSclusterV2.h"
  23   #include "AliITSsimulationFastPoints.h"
  24   #include "AliITStrackerV2.h"
  25 #endif
  26
  27 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
  28 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n);
  29 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname,  Int_t n);
  30 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname);
  31
  32 Int_t ITSFindClusters(const Char_t *inname,  const Char_t *outname, Int_t n);
  33 Int_t ITSFindTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
  34 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname);
  35
  36
  37 Int_t AliBarrelReconstruction(Int_t n=1) {
  38    const Char_t *TPCdigName="rfio:galice.root";
  39    const Char_t *TPCclsName="AliTPCclusters.root";
  40    const Char_t *TPCtrkName="AliTPCtracks.root";
  41    const Char_t *TPCtrkNameS="AliTPCtracksSorted.root";
  42
  43
  44    const Char_t *ITSdigName="rfio:galice.root";
  45    const Char_t *ITSclsName="AliITSclustersV2.root";
  46    const Char_t *ITStrkName="AliITStracksV2.root";
  47
  48    AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
  49
  50
  51 // ********** Find TPC clusters *********** //
  52    if (TPCFindClusters(TPCdigName,TPCclsName, n)) {
  53       cerr<<"Failed to get TPC clusters !\n";
  54       return 1;
  55    }
  56
  57 // ********** Find TPC tracks *********** //
  58     if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
  59       cerr<<"Failed to get TPC tracks !\n";
  60       return 1;
  61     }
  62
  63 // ********** Sort and label TPC tracks *********** //
  64    if (TPCSortTracks(TPCclsName,TPCtrkName,TPCtrkNameS,n)) {
  65       cerr<<"Failed to sort TPC tracks !\n";
  66       return 1;
  67     }
  68
  69 // ********** Find ITS clusters *********** //
  70    if (ITSFindClusters(ITSdigName,ITSclsName,n)) {
  71       cerr<<"Failed to get ITS clusters !\n";
  72       return 1;
  73    }
  74
  75 // ********** Find ITS tracks *********** //
  76    {
  77      //TFile *clsFile=TFile::Open(ITSclsName);
  78    if (ITSFindTracks(TPCtrkNameS,ITSclsName,ITStrkName,n)) {
  79       cerr<<"Failed to get ITS tracks !\n";
  80       return 1;
  81    }
  82    //clsFile->Close();
  83    }
  84    return 1;
  85 // ********** Back propagation of the ITS tracks *********** //
  86    {TFile *clsFile=TFile::Open(ITSclsName);
  87    if (ITSPropagateBack(ITStrkName,TPCtrkNameS)) {
  88       cerr<<"Failed to propagate back the ITS tracks !\n";
  89       return 1;
  90    }
  91    clsFile->Close();}
  92
  93
  94 // ********** Back propagation of the TPC tracks *********** //
  95    {TFile *clsFile=TFile::Open(TPCclsName);
  96    if (TPCPropagateBack(TPCtrkName,TPCtrkName)) {
  97       cerr<<"Failed to propagate back the TPC tracks !\n";
  98       return 1;
  99    }
 100    clsFile->Close();}
 101
 102    return 0;
 103 }
 104
 105
 106 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
 107    Int_t rc=0;
 108    const Char_t *name="TPCFindClusters";
 109    cerr<<'\n'<<name<<"...\n";
 110    gBenchmark->Start(name);
 111
 112    TFile *out=TFile::Open(outname,"recreate");
 113    TFile *in =TFile::Open(inname);
 114
 115    AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
 116    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
 117    AliTPCv2 tpc;
 118    tpc.SetParam(param);
 119
 120    //tpc.Digits2Clusters(out); //MI change
 121    for (Int_t i=0;i<n;i++){
 122      printf("Processing event %d\n",i);
 123      tpc.Digits2Clusters(out,i);
 124      //  AliTPCclusterer::Digits2Clusters(dig, out, i);
 125    }
 126    in->Close();
 127    out->Close();
 128    gBenchmark->Stop(name);
 129    gBenchmark->Show(name);
 130
 131    return rc;
 132 }
 133
 134 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n) {
 135    Int_t rc=0;
 136    const Char_t *name="TPCFindTracks";
 137    cerr<<'\n'<<name<<"...\n";
 138    gBenchmark->Start(name);
 139    TFile *out=TFile::Open(outname,"recreate");
 140    TFile *in =TFile::Open(inname);
 141    AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
 142    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
 143
 144    //AliTPCtracker *tracker=new AliTPCtracker(param);
 145    //rc=tracker->Clusters2Tracks(0,out);
 146    //delete tracker;
 147
 148    for (Int_t i=0;i<n;i++){
 149      printf("Processing event %d\n",i);
 150      AliTPCtracker *tracker = new AliTPCtracker(param,i);
 151      //Int_t rc=
 152        tracker->Clusters2Tracks(0,out);
 153      delete tracker;
 154    }
 155
 156    in->Close();
 157    out->Close();
 158    gBenchmark->Stop(name);
 159    gBenchmark->Show(name);
 160
 161    return rc;
 162 }
 163 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t * outname,  Int_t eventn){
 164    Int_t rc=0;
 165    const Char_t *name="TPCSortTracks";
 166    cerr<<'\n'<<name<<"...\n";
 167    gBenchmark->Start(name);
 168
 169    TFile *out =TFile::Open(outname,"recreate");
 170    TFile *in1 =TFile::Open(inname);
 171    TFile *in2 =TFile::Open(inname2);
 172
 173
 174    AliTPCParam *param=(AliTPCParam *)in1->Get("75x40_100x60");
 175    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
 176
 177
 178
 179    // loop over events
 180    for (Int_t event=0;event<eventn; event++){
 181      TObjArray tarray(10000);
 182      AliTPCtrack *iotrack=0;
 183      Int_t i;
 184
 185
 186      in2->cd();
 187      char   tname[100];
 188      sprintf(tname,"TreeT_TPC_%d",event);
 189
 190      TTree *tracktree=(TTree*)in2->Get(tname);
 191      TBranch *tbranch=tracktree->GetBranch("tracks");
 192      Int_t nentr=(Int_t)tracktree->GetEntries();
 193      for (i=0; i<nentr; i++) {
 194        iotrack=new AliTPCtrack;
 195        tbranch->SetAddress(&iotrack);
 196        tracktree->GetEvent(i);
 197        tarray.AddLast(iotrack);
 198      }
 199      tarray.Sort();
 200      // in2->Close();
 201
 202      //assign thacks GEANT labels
 203      in1->cd();
 204      AliTPCtracker *tracker = new AliTPCtracker(param,event);
 205      tracker->LoadInnerSectors();
 206      tracker->LoadOuterSectors();
 207      for (i=0; i<nentr; i++) {
 208        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
 209        tracker->CookLabel(iotrack,0.1);
 210      }
 211      delete tracker;
 212      //in->Close();
 213      //end of GEANT label assignment
 214
 215      tracktree=new TTree(tname,"Tree with TPC tracks");
 216      tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 217      for (i=0; i<nentr; i++) {
 218        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
 219        tracktree->Fill();
 220      }
 221      out->cd();
 222      tracktree->Write();
 223    }
 224
 225    out->Close();
 226    in2->Close();
 227    in1->Close();
 228
 229    gBenchmark->Stop(name);
 230    gBenchmark->Show(name);
 231
 232    return rc;
 233 }
 234
 235 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
 236    Int_t rc=0;
 237    const Char_t *name="ITSFindClusters";
 238    cerr<<'\n'<<name<<"...\n";
 239    gBenchmark->Start(name);
 240    TFile *out=TFile::Open(outname,"recreate");
 241    TFile *in =TFile::Open(inname,"update");
 242
 243    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
 244       cerr<<"Can't get gAlice !\n";
 245       return 1;
 246    }
 247
 248    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 249    if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
 250    AliITSgeom *geom=ITS->GetITSgeom();
 251    out->cd();
 252    geom->Write();
 253
 254    Int_t ev=0;
 255    for (ev = 0; ev<n; ev++){
 256      in->cd();   // !!!! MI directory must point to galice. - othervise problem with Tree -connection
 257      gAlice->GetEvent(ev);
 258      //gAlice->TreeR()->Reset();   //reset reconstructed tree
 259
 260
 261      TTree *pTree=gAlice->TreeR();
 262      if (!pTree){
 263        gAlice->MakeTree("R");
 264        pTree = gAlice->TreeR();
 265      }
 266      TBranch *branch=pTree->GetBranch("ITSRecPoints");
 267      if (!branch) {
 268        //if not reconstructed ITS branch do reconstruction
 269        ITS->MakeBranch("R",0);
 270        //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
 271        for (Int_t i=0;i<3;i++) {
 272          ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
 273        }
 274        Int_t nsignal=25;
 275        Int_t size=-1;
 276        Int_t bgr_ev=Int_t(ev/nsignal);
 277        ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
 278        //////////////////////////////////////////////////////////////////////////
 279        gAlice->GetEvent(ev);   //MI comment  - in HitsToFast... they reset treeR to 0
 280        //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
 281        //time
 282      }
 283
 284
 285
 286      out->cd();
 287      TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
 288      char   cname[100];
 289      sprintf(cname,"TreeC_ITS_%d",ev);
 290
 291      TTree *cTree=new TTree(cname,"ITS clusters");
 292      cTree->Branch("Clusters",&clusters);
 293
 294      pTree=gAlice->TreeR();
 295      if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
 296      branch=pTree->GetBranch("ITSRecPoints");
 297      if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
 298      TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
 299      branch->SetAddress(&points);
 300
 301      TClonesArray &cl=*clusters;
 302      Int_t nclusters=0;
 303      Int_t nentr=(Int_t)pTree->GetEntries();
 304      AliITSgeom *geom=ITS->GetITSgeom();
 305
 306      for (Int_t i=0; i<nentr; i++) {
 307        if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
 308        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
 309        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
 310        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
 311        Double_t yshift = x*rot[0] + y*rot[1];
 312        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
 313        Int_t ncl=points->GetEntriesFast();
 314        nclusters+=ncl;
 315        for (Int_t j=0; j<ncl; j++) {
 316          AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
 317          Float_t lp[5];
 318          lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
 319          lp[1]=p->GetZ()+zshift;
 320          lp[2]=p->GetSigmaX2();
 321          lp[3]=p->GetSigmaZ2();
 322          lp[4]=p->GetQ();
 323          Int_t lab[6];
 324          lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
 325          lab[3]=ndet;
 326
 327          Int_t label=lab[0];
 328          TParticle *part=(TParticle*)gAlice->Particle(label);
 329          label=-3;
 330          while (part->P() < 0.005) {
 331            Int_t m=part->GetFirstMother();
 332            if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
 333            label=m;
 334            part=(TParticle*)gAlice->Particle(label);
 335          }
 336          if      (lab[1]<0) lab[1]=label;
 337          else if (lab[2]<0) lab[2]=label;
 338          else cerr<<"No empty labels !\n";
 339
 340          new(cl[j]) AliITSclusterV2(lab,lp);
 341        }
 342        cTree->Fill(); clusters->Delete();
 343        points->Delete();
 344      }
 345      cTree->Write();
 346      cerr<<"Number of clusters: "<<nclusters<<endl;
 347      delete cTree; delete clusters; delete points;
 348
 349    }
 350
 351    delete gAlice; gAlice=0;
 352    in->Close();
 353    out->Close();
 354    gBenchmark->Stop(name);
 355    gBenchmark->Show(name);
 356
 357    return rc;
 358 }
 359
 360 Int_t ITSFindTracks(const Char_t * inname, const Char_t *inname2, const Char_t *outname, Int_t n) {
 361    Int_t rc=0;
 362    const Char_t *name="ITSFindTracks";
 363    cerr<<'\n'<<name<<"...\n";
 364    gBenchmark->Start(name);
 365
 366
 367    TFile *out=TFile::Open(outname,"recreate");
 368    TFile *in =TFile::Open(inname);
 369    TFile *in2 =TFile::Open(inname2);
 370
 371    AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
 372    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
 373
 374
 375    for (Int_t i=0;i<n;i++){
 376      AliITStrackerV2 tracker(geom,i);
 377      rc=tracker.Clusters2Tracks(in,out);
 378    }
 379
 380    in->Close();
 381    in2->Close();
 382    out->Close();
 383
 384    gBenchmark->Stop(name);
 385    gBenchmark->Show(name);
 386
 387    return rc;
 388 }
 389
 390 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname) {
 391
 392   Int_t rc=0;
 393   /*
 394    const Char_t *name="ITSPropagateBack";
 395    cerr<<'\n'<<name<<"...\n";
 396    gBenchmark->Start(name);
 397
 398    AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
 399    if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
 400    AliITStrackerV2 tracker(geom);
 401
 402    TFile *out=TFile::Open(outname,"update");
 403    TFile *in =TFile::Open(inname);
 404    rc=tracker.PropagateBack(in,out);
 405    in->Close();
 406    out->Close();
 407
 408    gBenchmark->Stop(name);
 409    gBenchmark->Show(name);
 410   */
 411    return rc;
 412 }
 413
 414 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname) {
 415    Int_t rc=0;
 416    const Char_t *name="TPCPropagateBack";
 417    cerr<<'\n'<<name<<"...\n";
 418    gBenchmark->Start(name);
 419
 420    AliTPCParam *param=(AliTPCParam *)gFile->Get("75x40_100x60");
 421    if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
 422    AliTPCtracker *tracker=new AliTPCtracker(param);
 423
 424    TFile *out=TFile::Open(outname,"update");
 425    TFile *in =TFile::Open(inname);
 426    rc=tracker->PropagateBack(in,out);
 427    delete tracker;
 428    in->Close();
 429    out->Close();
 430
 431    gBenchmark->Stop(name);
 432    gBenchmark->Show(name);
 433
 434    return rc;
 435 }
 436