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[7000];
  75    Int_t ngood=0;
  76
  77    cerr<<"Marking good tracks (this will take a while)...\n";
  78    ngood=good_tracks(gt,7000);
  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    TClonesArray *particles=gAlice->Particles();
 127    Int_t np=particles->GetEntriesFast();
 128    Int_t *good=new Int_t[np];
 129    for (Int_t ii=0; ii<np; ii++) good[ii]=0;
 130
 131    //MI change to be possible compile macro
 132    //definition out of the switch statement
 133    Int_t sectors_by_rows=0;
 134    TTree *TD=0;
 135    AliSimDigits da, *digits=&da;
 136    Int_t *count=0;
 137
 138    switch (ver) {
 139    case 1:
 140      {
 141       TFile *cf=TFile::Open("AliTPCclusters.root");
 142       if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
 143       AliTPCClustersArray *ca=new AliTPCClustersArray;
 144       ca->Setup(digp);
 145       ca->SetClusterType("AliTPCcluster");
 146       ca->ConnectTree("Segment Tree");
 147       Int_t nrows=Int_t(ca->GetTree()->GetEntries());
 148       for (Int_t n=0; n<nrows; n++) {
 149           AliSegmentID *s=ca->LoadEntry(n);
 150           Int_t sec,row;
 151           digp->AdjustSectorRow(s->GetID(),sec,row);
 152           AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
 153           Int_t ncl=clrow.GetArray()->GetEntriesFast();
 154           while (ncl--) {
 155               AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
 156               Int_t lab=c->GetLabel(0);
 157               if (lab<0) continue; //noise cluster
 158               lab=TMath::Abs(lab);
 159               if (sec>=digp->GetNInnerSector())
 160               if (row==nrow_up-1    ) good[lab]|=0x1000;
 161               if (sec>=digp->GetNInnerSector())
 162               if (row==nrow_up-1-gap) good[lab]|=0x800;
 163               good[lab]++;
 164           }
 165           ca->ClearRow(sec,row);
 166       }
 167       cf->Close();
 168      }
 169       break;
 170    case 2:
 171       TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
 172       TD->GetBranch("Segment")->SetAddress(&digits);
 173       count = new Int_t[np];
 174       Int_t i;
 175       for (i=0; i<np; i++) count[i]=0;
 176       Int_t sectors_by_rows=(Int_t)TD->GetEntries();
 177       for (i=0; i<sectors_by_rows; i++) {
 178           if (!TD->GetEvent(i)) continue;
 179           Int_t sec,row;
 180           digp->AdjustSectorRow(digits->GetID(),sec,row);
 181           cerr<<sec<<' '<<row<<"                                     \r";
 182           digits->First();
 183           while (digits->Next()) {
 184               Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
 185               Short_t dig = digits->GetDigit(it,ip);
 186               Int_t idx0=digits->GetTrackID(it,ip,0);
 187               Int_t idx1=digits->GetTrackID(it,ip,1);
 188               Int_t idx2=digits->GetTrackID(it,ip,2);
 189               if (idx0>=0 && dig>=zero) count[idx0]+=1;
 190               if (idx1>=0 && dig>=zero) count[idx1]+=1;
 191               if (idx2>=0 && dig>=zero) count[idx2]+=1;
 192           }
 193           for (Int_t j=0; j<np; j++) {
 194               if (count[j]>1) {
 195                  if (sec>=digp->GetNInnerSector())
 196                  if (row==nrow_up-1    ) good[j]|=0x1000;
 197                  if (sec>=digp->GetNInnerSector())
 198                  if (row==nrow_up-1-gap) good[j]|=0x800;
 199                  good[j]++;
 200               }
 201               count[j]=0;
 202           }
 203       }
 204       delete[] count;
 205       break;
 206    default:
 207       cerr<<"Invalid TPC version !\n";
 208       file->Close();
 209       exit(7);
 210    }
 211
 212    TTree *TH=gAlice->TreeH();
 213    //TClonesArray *hits=TPC->Hits();
 214    Int_t npart=TH->GetEntries();
 215
 216    while (npart--) {
 217       AliTPChit *hit0=0;
 218       TPC->ResetHits();
 219       TH->GetEvent(npart);
 220       AliTPChit *hit = (AliTPChit*) TPC->FirstHit(-1);
 221
 222       while(hit) {
 223          if(hit->fQ==0.) break;
 224          hit = (AliTPChit*) TPC->NextHit();
 225       }
 226       if(hit) {
 227          hit0 = new AliTPChit(*hit); // Make copy of hit
 228          hit=hit0;
 229       }
 230       else continue;
 231       AliTPChit *hit1=(AliTPChit*) TPC->NextHit();
 232       if(hit1==0) continue;
 233
 234       if (hit1->fQ != 0.) continue;
 235      // Int_t i=hit->fTrack;
 236            Int_t i=hit->Track();  //modificata in accordo a nuovo AliTPCComparison
 237       TParticle *p = (TParticle*)particles->UncheckedAt(i);
 238       if (p->GetFirstMother()>=0) continue;  //secondary particle
 239       if (good[i] < 0x1000+0x800+2+good_number) continue;
 240       if (p->Pt()<0.100) continue;
 241       if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
 242
 243       gt[nt].lab=i;
 244       gt[nt].code=p->GetPdgCode();
 245 //**** px py pz - in global coordinate system, x y z - in local !
 246      // gt[nt].px=hit->fX; gt[nt].py=hit->fY; gt[nt].pz=hit->fZ;  //modificato tenendo conto di AliTPCComparison
 247       gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();
 248       Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
 249       //gt[nt].x = hit1->fX*cs + hit1->fY*sn;
 250      // gt[nt].y =-hit1->fX*sn + hit1->fY*cs;  //modificato tenedo conto di AliTPCComaprison
 251       //gt[nt].z = hit1->fZ;
 252       gt[nt].x = hit1->X()*cs + hit1->Y()*sn;
 253       gt[nt].y =-hit1->X()*sn + hit1->Y()*cs;
 254       gt[nt].z = hit1->Z();
 255       gt[nt].pxg = p->Px();
 256       gt[nt].pyg = p->Py();
 257       gt[nt].pzg = p->Pz();
 258       gt[nt].ptg = p->Pt();
 259       gt[nt].flag = 0;
 260       nt++;
 261
 262         if(hit0) delete hit0;
 263       cerr<<i<<"                \r";
 264       if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
 265    }
 266    delete[] good;
 267
 268    printf("before delete gAlice\n");
 269
 270    delete gAlice; gAlice=0;
 271
 272    printf("after delete gAlice\n");
 273    file->Close();
 274    printf("after file close\n");
 275    return nt;
 276 }
 277
 278 //--------------------------------------
 279