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