ITS/TPCtracks.C

   1 struct GoodTrack {
   2   Int_t lab;
   3   Int_t code;
   4   Float_t px,py,pz;
   5   Float_t x,y,z;
   6   Float_t pxg,pyg,pzg,ptg;
   7   Bool_t flag;
   8 };
   9 Int_t good_tracks(GoodTrack *gt, Int_t max);
  10
  11
  12 Int_t TPCtracks() {
  13    Int_t i;
  14    cerr<<"Doing comparison...\n";
  15
  16   // Connect the Root Galice file containing Geometry, Kine and Hits
  17
  18
  19    TFile *cf=TFile::Open("AliTPCclusters.root");
  20    if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
  21    AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
  22    if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }
  23
  24 // Load clusters
  25    AliTPCClustersArray *ca=new AliTPCClustersArray;
  26    ca->Setup(digp);
  27    ca->SetClusterType("AliTPCcluster");
  28    ca->ConnectTree("Segment Tree");
  29    Int_t nentr=Int_t(ca->GetTree()->GetEntries());
  30    for (Int_t i=0; i<nentr; i++) {
  31        ca->LoadEntry(i);
  32    }
  33
  34 // Load tracks
  35    TFile *tf=TFile::Open("AliTPCtracks.root");
  36    if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
  37    TObjArray tarray(2000);
  38    TTree *tracktree=(TTree*)tf->Get("TreeT");
  39    TBranch *tbranch=tracktree->GetBranch("tracks");
  40    nentr=(Int_t)tracktree->GetEntries();
  41    for (i=0; i<nentr; i++) {
  42        AliTPCtrack *iotrack=new AliTPCtrack;
  43        tbranch->SetAddress(&iotrack);
  44        tracktree->GetEvent(i);
  45        iotrack->CookLabel(ca);
  46        tarray.AddLast(iotrack);
  47    }
  48
  49    tf->Close();
  50    delete ca;
  51    cf->Close();
  52    //printf("after cf close\n");
  53
  54    cerr<<"Number of found tracks "<<nentr<<endl;
  55    tarray.Sort();
  56
  57    TFile *pfile = new TFile("tpctracks.root","RECREATE");
  58
  59    TTree tracktree1("TreeT","Tree with TPC tracks");
  60    AliTPCtrack *iotrack=0;
  61    tracktree1.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
  62
  63    for (i=0; i<nentr; i++) {
  64        iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
  65        tracktree1.Fill();
  66    }
  67
  68    tracktree1.Write();
  69    pfile->Close();
  70    tarray.Clear();
  71
  72
  73 /////////////////////////////////////////////////////////////////////////
  74    GoodTrack gt[15000];
  75    Int_t ngood=0;
  76
  77    cerr<<"Marking good tracks (this will take a while)...\n";
  78    ngood=good_tracks(gt,15000);
  79    ofstream out("good_tracks");
  80    if (out) {
  81          for (Int_t ngd=0; ngd<ngood; ngd++)
  82             out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
  83                <<gt[ngd].px <<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
  84                <<gt[ngd].x  <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<' '
  85                <<gt[ngd].pxg  <<' '<<gt[ngd].pyg <<' '<<gt[ngd].pzg <<' '
  86                <<gt[ngd].ptg <<gt[ngd].flag<<endl;
  87    } else cerr<<"Can not open file (good_tracks) !\n";
  88    out.close();
  89
  90    cerr<<"Number of good tracks : "<<ngood<<endl;
  91
  92    printf("before return in tpctracks\n");
  93    return 0;
  94 }
  95
  96 //---------------------------------
  97
  98 Int_t good_tracks(GoodTrack *gt, Int_t max) {
  99    Int_t nt=0;
 100
 101    //TFile *file=TFile::Open("rfio:galice.root");  // via server
 102    TFile *file=TFile::Open("galice.root");
 103    if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
 104
 105    if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
 106      cerr<<"gAlice have not been found on galice.root !\n";
 107      exit(5);
 108    }
 109
 110    gAlice->GetEvent(0);
 111
 112    AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
 113    Int_t ver = TPC->IsVersion();
 114    cerr<<"TPC version "<<ver<<" has been found !\n";
 115
 116    AliTPCParam *digp=(AliTPCParam*)file->Get("75x40_100x60");
 117    if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
 118    TPC->SetParam(digp);
 119
 120    Int_t nrow_up=digp->GetNRowUp();
 121    Int_t nrows=digp->GetNRowLow()+nrow_up;
 122    Int_t zero=digp->GetZeroSup();
 123    Int_t gap=Int_t(0.125*nrows);
 124    Int_t good_number=Int_t(0.4*nrows);
 125
 126    TObjArray *particles=gAlice->Particles();
 127    //Int_t np=particles->GetEntriesFast();
 128    Int_t np=gAlice->GetNtrack(); //FCA correction
 129    Int_t *good=new Int_t[np];
 130    for (Int_t ii=0; ii<np; ii++) good[ii]=0;
 131
 132    //MI change to be possible compile macro
 133    //definition out of the switch statement
 134    Int_t sectors_by_rows=0;
 135    TTree *TD=0;
 136    AliSimDigits da, *digits=&da;
 137    Int_t *count=0;
 138
 139    switch (ver) {
 140    case 1:
 141      {
 142       TFile *cf=TFile::Open("AliTPCclusters.root");
 143       if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
 144       AliTPCClustersArray *ca=new AliTPCClustersArray;
 145       ca->Setup(digp);
 146       ca->SetClusterType("AliTPCcluster");
 147       ca->ConnectTree("Segment Tree");
 148       Int_t nrows=Int_t(ca->GetTree()->GetEntries());
 149       for (Int_t n=0; n<nrows; n++) {
 150           AliSegmentID *s=ca->LoadEntry(n);
 151           Int_t sec,row;
 152           digp->AdjustSectorRow(s->GetID(),sec,row);
 153           AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
 154           Int_t ncl=clrow.GetArray()->GetEntriesFast();
 155           while (ncl--) {
 156               AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
 157               Int_t lab=c->GetLabel(0);
 158               if (lab<0) continue; //noise cluster
 159               lab=TMath::Abs(lab);
 160               if (sec>=digp->GetNInnerSector())
 161               if (row==nrow_up-1    ) good[lab]|=0x1000;
 162               if (sec>=digp->GetNInnerSector())
 163               if (row==nrow_up-1-gap) good[lab]|=0x800;
 164               good[lab]++;
 165           }
 166           ca->ClearRow(sec,row);
 167       }
 168       cf->Close();
 169      }
 170       break;
 171    case 2:
 172       TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
 173       TD->GetBranch("Segment")->SetAddress(&digits);
 174       count = new Int_t[np];
 175       Int_t i;
 176       for (i=0; i<np; i++) count[i]=0;
 177       Int_t sectors_by_rows=(Int_t)TD->GetEntries();
 178       for (i=0; i<sectors_by_rows; i++) {
 179           if (!TD->GetEvent(i)) continue;
 180           Int_t sec,row;
 181           digp->AdjustSectorRow(digits->GetID(),sec,row);
 182           cerr<<sec<<' '<<row<<"                                     \r";
 183           digits->First();
 184           while (digits->Next()) {
 185               Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
 186               Short_t dig = digits->GetDigit(it,ip);
 187               Int_t idx0=digits->GetTrackID(it,ip,0);
 188               Int_t idx1=digits->GetTrackID(it,ip,1);
 189               Int_t idx2=digits->GetTrackID(it,ip,2);
 190               if (idx0>=0 && dig>=zero) count[idx0]+=1;
 191               if (idx1>=0 && dig>=zero) count[idx1]+=1;
 192               if (idx2>=0 && dig>=zero) count[idx2]+=1;
 193           }
 194           for (Int_t j=0; j<np; j++) {
 195               if (count[j]>1) {
 196                  if (sec>=digp->GetNInnerSector())
 197                  if (row==nrow_up-1    ) good[j]|=0x1000;
 198                  if (sec>=digp->GetNInnerSector())
 199                  if (row==nrow_up-1-gap) good[j]|=0x800;
 200                  good[j]++;
 201               }
 202               count[j]=0;
 203           }
 204       }
 205       delete[] count;
 206       break;
 207    default:
 208       cerr<<"Invalid TPC version !\n";
 209       file->Close();
 210       exit(7);
 211    }
 212
 213    TTree *TH=gAlice->TreeH();
 214    //TClonesArray *hits=TPC->Hits();
 215    Int_t npart=TH->GetEntries();
 216
 217    while (npart--) {
 218       AliTPChit *hit0=0;
 219       TPC->ResetHits();
 220       TH->GetEvent(npart);
 221       AliTPChit *hit = (AliTPChit*) TPC->FirstHit(-1);
 222
 223       while(hit) {
 224          if(hit->fQ==0.) break;
 225          hit = (AliTPChit*) TPC->NextHit();
 226       }
 227       if(hit) {
 228          hit0 = new AliTPChit(*hit); // Make copy of hit
 229          hit=hit0;
 230       }
 231       else continue;
 232       AliTPChit *hit1=(AliTPChit*) TPC->NextHit();
 233       if(hit1==0) continue;
 234
 235       if (hit1->fQ != 0.) continue;
 236      // Int_t i=hit->fTrack;
 237            Int_t i=hit->Track();
 238       //TParticle *p = (TParticle*)particles->UncheckedAt(i);
 239           TParticle *p = (TParticle*)gAlice->Particle(i);  //mod. secondonew
 240                                                               //AliTPCcomaprison
 241       if (p->GetFirstMother()>=0) continue;  //secondary particle
 242       if (good[i] < 0x1000+0x800+2+good_number) continue;
 243       if (p->Pt()<0.100) continue;
 244       if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
 245
 246       gt[nt].lab=i;
 247       gt[nt].code=p->GetPdgCode();
 248 //**** px py pz - in global coordinate system, x y z - in local !
 249      // gt[nt].px=hit->fX; gt[nt].py=hit->fY; gt[nt].pz=hit->fZ;  //modificato tenendo conto di AliTPCComparison
 250       gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();
 251       Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
 252       //gt[nt].x = hit1->fX*cs + hit1->fY*sn;
 253      // gt[nt].y =-hit1->fX*sn + hit1->fY*cs;  //modificato tenedo conto di AliTPCComaprison
 254       //gt[nt].z = hit1->fZ;
 255       gt[nt].x = hit1->X()*cs + hit1->Y()*sn;
 256       gt[nt].y =-hit1->X()*sn + hit1->Y()*cs;
 257       gt[nt].z = hit1->Z();
 258       gt[nt].pxg = p->Px();
 259       gt[nt].pyg = p->Py();
 260       gt[nt].pzg = p->Pz();
 261       gt[nt].ptg = p->Pt();
 262       gt[nt].flag = 0;
 263       nt++;
 264
 265         if(hit0) delete hit0;
 266       cerr<<i<<"                \r";
 267       if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
 268    }
 269    delete[] good;
 270
 271    printf("before delete gAlice\n");
 272
 273    delete gAlice; gAlice=0;
 274
 275    printf("after delete gAlice\n");
 276    file->Close();
 277    printf("after file close\n");
 278    return nt;
 279 }
 280
 281 //--------------------------------------
 282