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