TPC/TestTPCTrackHits.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /// Macro to compare TClonesArray hits with interpolated hits
  17 /// ConvertHits1 read
  18 ///
  19 /// \author MI
  20
  21 #include "alles.h"
  22 #include "AliObjectArray.h"
  23 #include "AliTPCTrackHits.h"
  24 #include "AliArrayBranch.h"
  25 #include "TestTPCTrackHits.h"
  26
  27 /// \cond CLASSIMP
  28 ClassImp(AliTPChitD)
  29 /// \endcond
  30
  31 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits,  Bool_t debug, TClonesArray *arrd=0);
  32 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
  33
  34 void ConvertHits(const char * benchmark, Bool_t debug)
  35 {
  36
  37   /// convert - not parametrised hits stored in Clones array
  38   /// to
  39
  40   TFile f("galice.root","update");
  41   TClonesArray *arr = new TClonesArray("AliTPChit",100);
  42   TTree * treeCl =(TTree*)f.Get("TreeH0");
  43   TBranch *branch = treeCl->GetBranch("TPC");
  44   AliTPCTrackHits *  myhits = new AliTPCTrackHits;
  45   branch->SetAddress(&arr);
  46   //particle
  47   TTree * treeP = (TTree*)f.Get("TreeK0");
  48   TBranch *branchP = treeP->GetBranch("Particles");
  49   TClonesArray *arrp = new TClonesArray("TParticle",100);
  50   branchP->SetAddress(&arrp);
  51   branchP->GetEvent(0);
  52   //
  53   TFile f2("treeh.root","recreate");
  54   f2.SetCompressionLevel(10);
  55   TTree * treeh2 = new TTree("TreeTPCH","TreeTPCH");
  56   //treeh2->AliBranch("TPC","AliTPCTrackHits", &myhits,4096);
  57   treeh2->GetListOfBranches()->Add(new AliObjectBranch("TPC","AliTPCTrackHits",
  58          &myhits,treeh2,4096,1));
  59
  60   TFile f3("treehcl.root","recreate");
  61   f3.SetCompressionLevel(2);
  62   TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
  63   treeh3->Branch("TPC", &arr,64000,100);
  64
  65   gBenchmark->Start(benchmark);
  66   Int_t trsize = (Int_t)treeCl->GetEntries();
  67   for (Int_t i=0;i<300;i++){
  68     arr->Clear();
  69     Int_t size = branch->GetEvent(i);
  70     //if (size>0){
  71     if ((size>0)&& arr->GetEntriesFast()>0) {
  72       f.cd();
  73       myhits->Clear();
  74       myhits =MakeTrack(arr,arrp,myhits);
  75       CompareHits(arr,myhits,debug);
  76       f2.cd();
  77       treeh2->Fill();
  78       f3.cd();
  79       treeh3->Fill();
  80       if ((i%10)==0) cout<<i<<"\n";
  81     }
  82   }
  83   f3.cd();
  84   treeh3->Write();
  85   f2.cd();
  86   treeh2->Write();
  87   f2.Close();
  88   f.Close();
  89   delete myhits;
  90   gBenchmark->Stop(benchmark);
  91   gBenchmark->Show(benchmark);
  92 }
  93
  94
  95
  96 void CompareHits(const char * benchmark, Bool_t debug)
  97 {
  98   //compare persistent hits
  99   TFile f2("treeh.root");
 100   TTree * treeh2 = (TTree*)f2.Get("TreeTPCH");
 101   AliTPCTrackHits *  myhits = new AliTPCTrackHits ;
 102   AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC");
 103   mybranch->SetAddress(&myhits);
 104
 105
 106   TClonesArray *arr = new TClonesArray("AliTPChit",100);
 107   TFile f3("treehcl.root");
 108   TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
 109   TBranch *branch = treeh3->GetBranch("TPC");
 110   branch->SetAddress(&arr);
 111
 112   cout<<"Lets go!\n";
 113   gBenchmark->Start(benchmark);
 114   Int_t trsize = (Int_t)treeh3->GetEntries();
 115   for (Int_t i=0;i<300;i++){
 116     Int_t size = branch->GetEvent(i);
 117     mybranch->GetEvent(i);
 118     //if (arr){
 119     if ((arr->GetEntriesFast()>0) && size>0) {
 120       CompareHits(arr,myhits,debug);
 121       if ((i%10)==0) cout<<i<<"\n";
 122     }
 123   }
 124   gBenchmark->Stop(benchmark);
 125   gBenchmark->Show(benchmark);
 126 }
 127
 128
 129 void CompareHitsG(const char * benchmark, Bool_t debug)
 130 {
 131   //compare persistent hits
 132   TFile f2("galice.root");
 133   TTree * treeh2 = (TTree*)f2.Get("TreeH0");
 134   AliTPCTrackHits *  myhits = new AliTPCTrackHits ;
 135   AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC2");
 136   mybranch->SetAddress(&myhits);
 137
 138
 139   TClonesArray *arr = new TClonesArray("AliTPChit",100);
 140   //TFile f3("treehcl.root");
 141   //TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
 142   TBranch *branch = treeh2->GetBranch("TPC");
 143   branch->SetAddress(&arr);
 144
 145   TFile f3("treehdelta.root","recreate");
 146   f3.SetCompressionLevel(2);
 147   TTree * treeh3 = new TTree("DelataH","DeltaH");
 148   TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
 149   treeh3->Branch("TPC", &arrd,64000,100);
 150
 151   cout<<"Lets go!\n";
 152   gBenchmark->Start(benchmark);
 153   Int_t trsize = treeh2->GetEntries();
 154   for (Int_t i=0;i<trsize;i++){
 155     Int_t size = branch->GetEvent(i);
 156     mybranch->GetEvent(i);
 157     //if (arr){
 158     if ((arr->GetEntriesFast()>0) && size>0) {
 159       CompareHits(arr,myhits,debug,arrd);
 160     }
 161     if ((i%10)==0) cout<<i<<"\n";
 162     treeh3->Fill();
 163   }
 164   treeh3->Write();
 165   f3.Close();
 166   gBenchmark->Stop(benchmark);
 167   gBenchmark->Show(benchmark);
 168 }
 169
 170
 171
 172
 173 AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)
 174 {
 175   //
 176   //make track wit hits  according
 177   AliTPChit * hit;
 178   //  AliTPCTrackHits * myhits= new AliTPCTrackHits;
 179   myhits->SetHitPrecision(0.002);
 180   myhits->SetStepPrecision(0.003);
 181   myhits->SetMaxDistance(100);
 182   myhits->FlushHitStack(kTRUE);
 183   for (Int_t i=0;i<arr->GetEntriesFast();i++){
 184     hit = (AliTPChit*)arr->At(i);
 185     if (hit){
 186       TParticle *p = (TParticle*)arrp->At(hit->GetTrack());
 187       Float_t momentum = TMath::Sqrt(p->Px()*p->Px()+p->Py()*p->Py());
 188       Float_t ran= 100.*gRandom->Rndm();
 189       if (ran<1.)  myhits->FlushHitStack(kTRUE);
 190       if (momentum<0.01) {//if small momentum - not precise
 191         myhits->SetHitPrecision(0.05);
 192         myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
 193                              hit->Z(), hit->fQ+1000);
 194       }
 195       else {
 196         myhits->SetHitPrecision(0.002);
 197         myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
 198                              hit->Z(), hit->fQ);
 199       }
 200     }
 201   }
 202   myhits->FlushHitStack();
 203   return myhits;
 204 }
 205
 206
 207
 208 void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
 209 {
 210   //
 211   // if debug option  kTRUE
 212   // compare hits and write result  to the stdoutput
 213   AliTPChit * hit, *hit2;
 214   if (arrd) arrd->Clear();
 215
 216   for (Int_t i=0;i<arr->GetEntriesFast();i++){
 217     hit = (AliTPChit*)arr->At(i);
 218     if (hit){
 219       if (i==0) myhits->First();
 220       else myhits->Next();
 221       hit2 = myhits->GetHit();
 222       if (!hit2) {
 223         hit2=0;
 224         if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
 225       }
 226
 227       if (hit&&hit2){
 228         // if (hit){
 229
 230         AliTrackHitsParam *param= myhits->GetParam();
 231         AliHitInfo * info = myhits->GetHitInfo();
 232         //
 233         if (arrd) {
 234           TClonesArray &larrd = *arrd;
 235           AliTPChitD * h = new(larrd[i]) AliTPChitD;
 236           h->SetX(hit->X());
 237           h->SetY(hit->Y());
 238           h->SetZ(hit->Z());
 239           h->SetTrack(hit->GetTrack());
 240           h->fQ = hit->fQ;
 241           h->fSector = hit->fSector;
 242           AliTPChit * hitd = h->GetDelta();
 243           hitd->SetX(hit->X()-hit2->X());
 244           hitd->SetY(hit->Y()-hit2->Y());
 245           hitd->SetZ(hit->Z()-hit2->Z());
 246           hitd->SetTrack(hit->GetTrack()-hit2->GetTrack());
 247           hitd->fQ = hit->fQ-hit2->fQ;
 248           hitd->fSector = hit->fSector-hit2->fSector;
 249         }
 250
 251         if (debug){
 252           Float_t dd =
 253             TMath::Sqrt(
 254                         (hit->X()-hit2->X())*(hit->X()-hit2->X())+
 255                         (hit->Y()-hit2->Y())*(hit->Y()-hit2->Y())+
 256                         (hit->Z()-hit2->Z())*(hit->Z()-hit2->Z()));
 257           printf("C1\t%d\t%d\t%d\t%d\t",
 258                  hit->fSector,hit2->fSector,hit->GetTrack(),hit2->GetTrack());
 259           printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
 260                  hit->X(),hit2->X(), hit->Y(),hit2->Y(),hit->Z(),hit2->Z());
 261           printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
 262                  dd, param->fR,param->fZ,param->fFi,param->fAn,
 263                  param->fAd,param->fTheta);
 264           printf("%d\t%d\t%d\n",
 265                  (Int_t)info->fHitDistance, (Int_t)hit->fQ, (Int_t)hit2->fQ);
 266         }
 267       }
 268     }
 269   }
 270 }
 271
 272
 273
 274
 275 //Filter for paw visualisation of results
 276 //
 277 //sed filter
 278  //sed '/*/d' dd.txt | sed '/C/d'|  cat -n > out.txt
 279  // sed -n '/C2/p' dd.txt | sed 's/C1//' | cat -n >out2.txt
 280  // sed -n '/C3/p' dd.txt | sed 's/C2//' | cat -n >out3.txt
 281
 282 //filter
 283 //myfilter C1 dd.txt >out1.txt
 284 //myfilter C2 dd.txt >out2.txt
 285 //myfilter C3 dd.txt >out3.txt
 286 // sed -n 1,50000p
 287 // sed -n 1,50000p dd.txt | myfilter C1  | sed -n 1,50000p >out1.txt
 288
 289 //myfilter C1 dd.txt | sed -n 1,50000p >out1.txt
 290 //myfilter C2 dd.txt | sed -n 1,50000p >out2.txt
 291 //myfilter C3 dd.txt | sed -n 1,50000p >out3.txt
 292 /*
 293 paw visualisation
 294  Nt/Create 1 'Test of Match' 21 ! ! event v1 v2 t1 t2  x1 x2 y1  y2 z1 z2 dd r z fi an ad theta dist q1 q2
 295   Nt/Read 1 out1.txt
 296  nt
 297  Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r  z fi alfa theta dr1 ddr ddz ddrfi dd
 298   Nt/Read 2 out2.txt
 299
 300  Nt/Create 3 'Test of Match2' 2 ! ! event stacksize
 301   Nt/Read 3 out3.txt
 302 */
 303
 304 void   Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
 305             Double_t fSumX,  Double_t fSumX2, Double_t fSumX3,
 306             Double_t fSumX4, Int_t n,
 307             Double_t &a, Double_t &b, Double_t &c)
 308 {
 309   //
 310   //recalc parameters not fixing origin point
 311   Double_t det =
 312     n* (fSumX2*fSumX4-fSumX3*fSumX3) -
 313     fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
 314     fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
 315
 316   if (TMath::Abs(det)>0){
 317     a =
 318       (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
 319        fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
 320        fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
 321     b=
 322       (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
 323       fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
 324       fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
 325     c=
 326       (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
 327        fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
 328        fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
 329     cout<<a<<"\t"<<b<<"\t"<<c<<"\n";
 330   }
 331 }
 332
 333 void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
 334 {
 335   Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
 336   fSumY = fSumYX = fSumYX2 = fSumX =  fSumX2 = fSumX3 =  fSumX4 =0;
 337   Double_t a2,b2,c2;
 338   for (Int_t i=0;i<n; i++){
 339     Float_t  x = gRandom->Rndm();
 340     Float_t  y = a+b*x+c*x*x;
 341     fSumY+=y;
 342     fSumYX+=y*x;
 343     fSumYX2+=y*x*x;
 344     fSumX += x;
 345     fSumX2 += x*x;
 346     fSumX3 += x*x*x;
 347     fSumX4 += x*x*x*x;
 348   }
 349   Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);
 350 }