ITS/AliITSFindClustersV2.C

   1 #ifndef __CINT__
   2   #include <Riostream.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(Char_t SlowOrFast='f') {
  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 = 0;
  50    if (SlowOrFast=='f') {
  51      branch = pTree->GetBranch("ITSRecPointsF");
  52    }
  53    else {
  54      branch = pTree->GetBranch("ITSRecPoints");
  55    }
  56    if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 6; }
  57    TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
  58    branch->SetAddress(&points);
  59
  60    TClonesArray &cl=*clusters;
  61    Int_t nclusters=0;
  62    Int_t nentr=(Int_t)branch->GetEntries();
  63
  64    cerr<<"Number of entries: "<<nentr<<endl;
  65
  66    //   Float_t lp[5]; Int_t lab[6]; //Why can't it be inside a loop ?
  67    Float_t * lp = new Float_t[5];
  68    Int_t * lab = new Int_t[6];
  69
  70    for (Int_t i=0; i<nentr; i++) {
  71        points->Clear();
  72        branch->GetEvent(i);
  73        Int_t ncl=points->GetEntriesFast(); if (ncl==0){cTree->Fill();continue;}
  74        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
  75        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
  76        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
  77        Double_t yshift = x*rot[0] + y*rot[1];
  78        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
  79        nclusters+=ncl;
  80
  81        Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD
  82        if(lay==4 || lay==3){kmip=280.;};
  83        if(lay==6 || lay==5){kmip=38.;};
  84
  85        for (Int_t j=0; j<ncl; j++) {
  86           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
  87           //Float_t lp[5];
  88           lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
  89           lp[1]=p->GetZ()+zshift;
  90           lp[2]=p->GetSigmaX2();
  91           lp[3]=p->GetSigmaZ2();
  92           lp[4]=p->GetQ(); lp[4]/=kmip;
  93           //Int_t lab[6];
  94           lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
  95           lab[3]=ndet;
  96
  97           Int_t label=lab[0];
  98           if (label>=0) {
  99              TParticle *part=(TParticle*)gAlice->Particle(label);
 100              label=-3;
 101              while (part->P() < 0.005) {
 102                 Int_t m=part->GetFirstMother();
 103                 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
 104                 label=m;
 105                 part=(TParticle*)gAlice->Particle(label);
 106              }
 107              if      (lab[1]<0) lab[1]=label;
 108              else if (lab[2]<0) lab[2]=label;
 109              else cerr<<"No empty labels !\n";
 110           }
 111
 112           new(cl[j]) AliITSclusterV2(lab,lp);
 113        }
 114        cTree->Fill(); clusters->Delete();
 115        points->Delete();
 116    }
 117    cTree->Write();
 118
 119    cerr<<"Number of clusters: "<<nclusters<<endl;
 120
 121    delete [] lp;
 122    delete [] lab;
 123
 124    delete cTree; delete clusters; delete points;
 125
 126    delete gAlice; gAlice=0;
 127
 128    in->Close();
 129    out->Close();
 130
 131    return 0;
 132
 133 }
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147