ITS/AliITSFindClustersV2.C

   1 #ifndef __CINT__
   2   #include <iostream.h>
   3
   4   #include "AliRun.h"
   5   #include "AliITS.h"
   6   #include "AliITSgeom.h"
   7   #include "AliITSRecPoint.h"
   8   #include "AliITSclusterV2.h"
   9
  10   #include "TFile.h"
  11   #include "TTree.h"
  12   #include "TParticle.h"
  13 #endif
  14
  15 Int_t AliITSFindClustersV2() {
  16 /****************************************************************
  17  *  This macro converts AliITSRecPoint(s) to AliITSclusterV2(s) *
  18  ****************************************************************/
  19    cerr<<"AliITSRecPoint(s) -> AliITSclusterV2(s)...\n";
  20
  21    if (gAlice) {delete gAlice; gAlice=0;}
  22
  23    TFile *in=TFile::Open("galice.root");
  24    if (!in->IsOpen()) { cerr<<"Can't open galice.root !\n"; return 1; }
  25
  26    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
  27       cerr<<"Can't find gAlice !\n";
  28       return 2;
  29    }
  30    gAlice->GetEvent(0);
  31
  32    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
  33    if (!ITS) { cerr<<"Can't find the ITS !\n"; return 3; }
  34    AliITSgeom *geom=ITS->GetITSgeom();
  35
  36    TFile *out=TFile::Open("AliITSclustersV2.root","new");
  37    if (!out->IsOpen()) {
  38       cerr<<"Delete old AliITSclustersV2.root !\n";
  39       return 4;
  40    }
  41    geom->Write();
  42
  43    TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
  44    TTree *cTree=new TTree("TreeC_ITS_0","ITS clusters");
  45    cTree->Branch("Clusters",&clusters);
  46
  47    TTree *pTree=gAlice->TreeR();
  48    if (!pTree) { cerr<<"Can't get TreeR !\n"; return 5; }
  49    TBranch *branch=pTree->GetBranch("ITSRecPoints");
  50    if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 6; }
  51    TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
  52    branch->SetAddress(&points);
  53
  54    TClonesArray &cl=*clusters;
  55    Int_t nclusters=0;
  56    Int_t nentr=(Int_t)pTree->GetEntries();
  57
  58    cerr<<"Number of entries: "<<nentr<<endl;
  59
  60    Float_t lp[5]; Int_t lab[6]; //Why can't it be inside a loop ?
  61
  62    for (Int_t i=0; i<nentr; i++) {
  63        points->Clear();
  64        pTree->GetEvent(i);
  65        Int_t ncl=points->GetEntriesFast(); if (ncl==0){cTree->Fill();continue;}
  66        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
  67        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
  68        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
  69        Double_t yshift = x*rot[0] + y*rot[1];
  70        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
  71        nclusters+=ncl;
  72
  73        Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD
  74        if(lay==4 || lay==3){kmip=280.;};
  75        if(lay==6 || lay==5){kmip=38.;};
  76
  77        for (Int_t j=0; j<ncl; j++) {
  78           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
  79           //Float_t lp[5];
  80           lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
  81           lp[1]=p->GetZ()+zshift;
  82           lp[2]=p->GetSigmaX2();
  83           lp[3]=p->GetSigmaZ2();
  84           lp[4]=p->GetQ(); lp[4]/=kmip;
  85           //Int_t lab[6];
  86           lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
  87           lab[3]=ndet;
  88
  89           Int_t label=lab[0];
  90           if (label>=0) {
  91              TParticle *part=(TParticle*)gAlice->Particle(label);
  92              label=-3;
  93              while (part->P() < 0.005) {
  94                 Int_t m=part->GetFirstMother();
  95                 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
  96                 label=m;
  97                 part=(TParticle*)gAlice->Particle(label);
  98              }
  99              if      (lab[1]<0) lab[1]=label;
 100              else if (lab[2]<0) lab[2]=label;
 101              else cerr<<"No empty labels !\n";
 102           }
 103
 104           new(cl[j]) AliITSclusterV2(lab,lp);
 105        }
 106        cTree->Fill(); clusters->Delete();
 107        points->Delete();
 108    }
 109    cTree->Write();
 110
 111    cerr<<"Number of clusters: "<<nclusters<<endl;
 112
 113    delete cTree; delete clusters; delete points;
 114
 115    delete gAlice; gAlice=0;
 116
 117    in->Close();
 118    out->Close();
 119
 120    return 0;
 121
 122 }
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136