STEER/AliITSclusterComparison.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3
   4 #include "AliITStrackerMI.h"
   5 #include "AliITSclusterV2.h"
   6
   7 #include "AliITShit.h"
   8 #include "AliITS.h"
   9 #include "AliITSgeom.h"
  10 //#include "AliITSgeometryDetail.h"
  11 //#include "AliITSparameter.h"
  12 #include "alles.h"
  13 #include "TFile.h"
  14 #include "TStopwatch.h"
  15 #include "Rtypes.h"
  16 #include "TTree.h"
  17
  18
  19 #include "AliRunLoader.h"
  20 #include "AliStack.h"
  21 #include "TF1.h"
  22 #include "AliTrackReference.h"
  23 #endif
  24 #include "AliITSclusterComparison.h"
  25
  26
  27 ClassImp(AliITSCI)
  28 ClassImp(AliITSClusterErrAnal)
  29
  30
  31
  32 AliITSClusterErrAnal::AliITSClusterErrAnal(Char_t *chloader )
  33 {
  34   //
  35   //SET Input loaders
  36   if (gAlice){
  37     delete gAlice->GetRunLoader();
  38     delete gAlice;
  39     gAlice = 0x0;
  40   }
  41   fRunLoader = AliRunLoader::Open(chloader);
  42   if (fRunLoader == 0x0){
  43     cerr<<"Can not open session"<<endl;
  44     return;
  45   }
  46   fITSLoader = fRunLoader->GetLoader("ITSLoader");
  47   if (fITSLoader == 0x0){
  48     cerr<<"Can not get ITS Loader"<<endl;
  49     return ;
  50   }
  51   if (fRunLoader->LoadgAlice()){
  52     cerr<<"Error occured while l"<<endl;
  53     return;
  54   }
  55   fRunLoader->CdGAFile();
  56   fITS = (AliITS*)gAlice->GetDetector("ITS");
  57   if (!fITS) {
  58     cerr<<"AliITSclusterComparison.C : Can not find the ITS detector !"<<endl;
  59   }
  60   AliITSgeom *geom = fITS->GetITSgeom();
  61   //An instance of the ITS tracker
  62   fTracker = new AliITStrackerMI(geom);
  63   //
  64   //
  65   AliITSCI * clinfo = new AliITSCI();
  66   fFileA  = new TFile("itsclusteranal.root","recreate");
  67   fFileA->cd();
  68   fTreeA  = new TTree("itscl","itscl");
  69   fTreeA->Branch("itscl","AliITSCI",&clinfo);
  70
  71   AliITSclusterV2 * cl = new AliITSclusterV2;
  72   fTreeB  = new TTree("Clusters","Clusters");
  73   fTreeB->Branch("cl","AliITSclusterV2",&cl);
  74
  75   fClusters = 0;
  76   delete clinfo;
  77 }
  78
  79 void AliITSClusterErrAnal::SetIO(Int_t event)
  80 {
  81   //
  82   //set input output for given event
  83   fRunLoader->SetEventNumber(event);
  84   fRunLoader->LoadHeader();
  85   fRunLoader->LoadKinematics();
  86   fRunLoader->LoadTrackRefs();
  87   fITSLoader->LoadHits();
  88   fITSLoader->LoadRecPoints("read");
  89   //
  90   fStack = fRunLoader->Stack();
  91   fHitTree = fITSLoader->TreeH();
  92   fClusterTree = fITSLoader->TreeR();
  93   fReferenceTree = fRunLoader->TreeTR();
  94   //
  95 }
  96
  97 void AliITSClusterErrAnal::LoadClusters()
  98 {
  99   //
 100   //
 101   // Load clusters
 102   if (fClusters) {
 103     fClusters->Delete();
 104     delete fClusters;
 105   }
 106   Int_t nparticles = fStack->GetNtrack();
 107   fClusters = new TObjArray(nparticles);
 108   //
 109   TObjArray *ClusterArray = new TClonesArray("AliITSclusterV2");
 110   TObjArray carray(2000);
 111   TBranch *branch=fClusterTree->GetBranch("Clusters");
 112   if (!branch) {
 113     Error("ReadClusters","Can't get the branch !");
 114     return;
 115   }
 116   branch->SetAddress(&ClusterArray);
 117   TBranch *brcl = fTreeB->GetBranch("cl");
 118   AliITSclusterV2* cluster= new AliITSclusterV2;
 119   brcl->SetAddress(&cluster);
 120   Int_t nentries = (Int_t)fClusterTree->GetEntries();
 121   for (Int_t i=0;i<nentries;i++){
 122     fClusterTree->GetEvent(i);
 123     Int_t nCluster = ClusterArray->GetEntriesFast();
 124     SignDeltas(ClusterArray,-3.8);
 125     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
 126       cluster = (AliITSclusterV2*)ClusterArray->UncheckedAt(iCluster);
 127       carray.AddLast(new AliITSclusterV2(*cluster));
 128       fTreeB->Fill();
 129     }
 130   }
 131   //
 132   Int_t nClusters = carray.GetEntriesFast();
 133   printf("Total number of clusters %d \n", nClusters);
 134   fTracker->LoadClusters(fClusterTree);
 135   //
 136   //
 137   //SORT clusters
 138   //
 139   Int_t all=0;
 140   Int_t noisecl =0;
 141   for (Int_t i=0;i<nClusters;i++){
 142     AliITSclusterV2 *cl = (AliITSclusterV2 *) carray.UncheckedAt(i);
 143     if (cl->GetLabel(0)<0) noisecl++;
 144     //
 145     for (Int_t itrack=0; itrack<3;itrack++){
 146       Int_t lab = cl->GetLabel(itrack);
 147       if (lab>=0){
 148         TObjArray * array = (TObjArray*)fClusters->At(lab);
 149         if (array==0){
 150           array = new TObjArray(20);
 151           fClusters->AddAt(array,lab);
 152         }
 153         array->AddLast(cl);
 154         all++;
 155       }
 156     }
 157   }
 158   printf("Total number of track clusters %d \n", all);
 159   printf("Total number of noise clusters %d \n", noisecl);
 160 }
 161
 162 void AliITSClusterErrAnal::SignDeltas( TObjArray *ClusterArray, Float_t vz)
 163 {
 164   //
 165   //
 166   Int_t entries = ClusterArray->GetEntriesFast();
 167   if (entries<4) return;
 168   AliITSclusterV2* cluster = (AliITSclusterV2*)ClusterArray->At(0);
 169   Int_t layer = cluster->GetLayer();
 170   if (layer>1) return;
 171   Int_t index[10000];
 172   Int_t ncandidates=0;
 173   Float_t r = (layer>0)? 7:4;
 174   //
 175   for (Int_t i=0;i<entries;i++){
 176     AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(i);
 177     Float_t nz = 1+(cl0->GetZ()-vz)/r;
 178     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
 179     index[ncandidates] = i;  //candidate to belong to delta electron track
 180     ncandidates++;
 181     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
 182       cl0->SetDeltaProbability(1);
 183     }
 184   }
 185   //
 186   //
 187   //
 188   for (Int_t i=0;i<ncandidates;i++){
 189     AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(index[i]);
 190     if (cl0->GetDeltaProbability()>0.8) continue;
 191     //
 192     Int_t ncl = 0;
 193     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
 194     sumy=sumz=sumy2=sumyz=sumw=0.0;
 195     for (Int_t j=0;j<ncandidates;j++){
 196       if (i==j) continue;
 197       AliITSclusterV2* cl1 = (AliITSclusterV2*)ClusterArray->At(index[j]);
 198       //
 199       Float_t dz = cl0->GetZ()-cl1->GetZ();
 200       Float_t dy = cl0->GetY()-cl1->GetY();
 201       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
 202         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
 203         y[ncl] = cl1->GetY();
 204         z[ncl] = cl1->GetZ();
 205         sumy+= y[ncl]*weight;
 206         sumz+= z[ncl]*weight;
 207         sumy2+=y[ncl]*y[ncl]*weight;
 208         sumyz+=y[ncl]*z[ncl]*weight;
 209         sumw+=weight;
 210         ncl++;
 211       }
 212     }
 213     if (ncl<4) continue;
 214     Float_t det = sumw*sumy2  - sumy*sumy;
 215     Float_t delta=1000;
 216     if (TMath::Abs(det)>0.01){
 217       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
 218       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
 219       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
 220     }
 221     else{
 222       Float_t z0  = sumyz/sumy;
 223       delta = TMath::Abs(cl0->GetZ()-z0);
 224     }
 225     if ( delta<0.05) {
 226       cl0->SetDeltaProbability(1-20.*delta);
 227     }
 228   }
 229 }
 230
 231 void AliITSClusterErrAnal::LoadParticles()
 232 {
 233   //
 234   //
 235   //  fReferences = new TObjArray;
 236 }
 237
 238 void AliITSClusterErrAnal::SortReferences()
 239 {
 240   //
 241   //
 242   //
 243   printf("Sorting references\n");
 244   fReferences = new TObjArray;
 245   Int_t ntracks = fStack->GetNtrack();
 246   fReferences->Expand(ntracks);
 247   Int_t nentries = (Int_t)fReferenceTree->GetEntries();
 248   TClonesArray * arr = new TClonesArray("AliTrackReference");
 249   TBranch * br = fReferenceTree->GetBranch("ITS");
 250   br->SetAddress(&arr);
 251   //
 252   TClonesArray *labarr=0;
 253   Int_t nreferences=0;
 254   Int_t nreftracks=0;
 255   for (Int_t iprim=0;iprim<nentries;iprim++){
 256     if (br->GetEntry(iprim)){
 257       for (Int_t iref=0;iref<arr->GetEntriesFast();iref++){
 258         AliTrackReference *ref =(AliTrackReference*)arr->At(iref);
 259         if (!ref) continue;
 260         Int_t lab = ref->GetTrack();
 261         //if ( (lab<0) || (lab>ntracks)) continue;
 262         //
 263         if (fReferences->At(lab)==0) {
 264           labarr = new TClonesArray("AliTrackReference");
 265           fReferences->AddAt(labarr,lab);
 266           nreftracks++;
 267         }
 268         TClonesArray &larr = *labarr;
 269         new(larr[larr.GetEntriesFast()]) AliTrackReference(*ref);
 270         nreferences++;
 271       }
 272     }
 273   }
 274   printf("Total number of references = \t%d\n", nreferences);
 275   printf("Total number of tracks with references = \t%d\n", nreftracks);
 276   printf("End - Sorting references\n");
 277
 278 }
 279
 280
 281 void AliITSClusterErrAnal::GetNTeor(Int_t layer, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
 282 {
 283   //
 284   //get "mean shape"
 285   if (layer==0){
 286     ny = 1.+TMath::Abs(phi)*3.2;
 287     nz = 1.+TMath::Abs(theta)*0.34;
 288     return;
 289   }
 290   if (layer==1){
 291     ny = 1.+TMath::Abs(phi)*3.2;
 292     nz = 1.+TMath::Abs(theta)*0.28;
 293     return;
 294   }
 295
 296   if (layer>3){
 297     ny = 2.02+TMath::Abs(phi)*1.95;
 298     nz = 2.02+TMath::Abs(phi)*2.35;
 299     return;
 300   }
 301   ny  = 6.6-2.7*abs(phi);
 302   nz  = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
 303 }
 304
 305
 306 Int_t AliITSClusterErrAnal::GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi,Float_t expQ, Float_t &erry, Float_t &errz)
 307 {
 308   //calculate cluster position error
 309   //
 310   Float_t nz,ny;
 311   GetNTeor(layer, theta,phi,ny,nz);
 312   erry   = TMath::Sqrt(cl->GetSigmaY2());
 313   errz   = TMath::Sqrt(cl->GetSigmaZ2());
 314   //
 315   // PIXELS
 316   if (layer<2){
 317
 318     if (TMath::Abs(ny-cl->GetNy())>0.6)  {
 319       if (ny<cl->GetNy()){
 320         erry*=0.4+TMath::Abs(ny-cl->GetNy());
 321         errz*=0.4+TMath::Abs(ny-cl->GetNy());
 322       }else{
 323         erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
 324         errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
 325       }
 326     }
 327     if (TMath::Abs(nz-cl->GetNz())>1.)  {
 328       erry*=TMath::Abs(nz-cl->GetNz());
 329       errz*=TMath::Abs(nz-cl->GetNz());
 330     }
 331     erry*=0.85;
 332     errz*=0.85;
 333     erry= TMath::Min(erry,float(0.005));
 334     errz= TMath::Min(errz,float(0.03));
 335     return 10;
 336   }
 337   //STRIPS
 338   if (layer>3){
 339     if (cl->GetNy()==100||cl->GetNz()==100){
 340       erry = 0.004;
 341       errz = 0.2;
 342       return 100;
 343     }
 344     if (cl->GetNy()+cl->GetNz()>12){
 345       erry = 0.06;
 346       errz = 0.57;
 347       return 100;
 348     }
 349     Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
 350     Float_t chargematch = TMath::Max(double(normq/expQ),2.);
 351     //
 352     if (cl->GetType()==1 || cl->GetType()==10 ){
 353       if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
 354         errz = 0.043;
 355         erry = 0.00094;
 356         return 101;
 357       }
 358       if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
 359         errz = 0.06;
 360         erry =0.0013;
 361         return 102;
 362       }
 363       erry = 0.0027;
 364       errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15);
 365       return 103;
 366     }
 367     if (cl->GetType()==2 || cl->GetType()==11 ){
 368       erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05);
 369       errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5);
 370       return 104;
 371     }
 372
 373     if (cl->GetType()>100 ){
 374       if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
 375         errz = 0.05;
 376         erry = 0.00096;
 377         return 105;
 378       }
 379       if (cl->GetNy()+cl->GetNz()-nz-ny<1){
 380         errz = 0.10;
 381         erry = 0.0025;
 382         return 106;
 383       }
 384
 385       errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4);
 386       erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05);
 387       return 107;
 388     }
 389     Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
 390     if (diff<1) diff=1;
 391     if (diff>4) diff=4;
 392
 393     if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
 394       errz = 0.14*diff;
 395       erry = 0.003*diff;
 396       return 108;
 397     }
 398     erry = 0.04*diff;
 399     errz = 0.06*diff;
 400     return 109;
 401   }
 402
 403   //DRIFTS
 404   Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
 405   Float_t chargematch = normq/expQ;
 406   Float_t factorz=1;
 407   Int_t   cnz = cl->GetNz()%10;
 408   //charge match
 409   if (cl->GetType()==1){
 410     if (chargematch<1.25){
 411       erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
 412     }
 413     else{
 414       erry = 0.003*chargematch;
 415       if (cl->GetNz()==3) erry*=1.5;
 416     }
 417     if (chargematch<1.0){
 418       errz =  0.0011*(1.+6./cl->GetQ());
 419     }
 420     else{
 421       errz = 0.002*(1+2*(chargematch-1.));
 422     }
 423     if (cnz>nz+0.6) {
 424       erry*=(cnz-nz+0.5);
 425       errz*=1.4*(cnz-nz+0.5);
 426     }
 427   }
 428   if (cl->GetType()>1){
 429     if (chargematch<1){
 430       erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
 431       errz =  0.0016*(1.+6./cl->GetQ());
 432     }
 433     else{
 434       errz = 0.0014*(1+3*(chargematch-1.));
 435       erry = 0.003*(1+3*(chargematch-1.));
 436     }
 437     if (cnz>nz+0.6) {
 438       erry*=(cnz-nz+0.5);
 439       errz*=1.4*(cnz-nz+0.5);
 440     }
 441   }
 442
 443   if (TMath::Abs(cl->GetY())>2.5){
 444     factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
 445   }
 446   if (TMath::Abs(cl->GetY())<1){
 447     factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
 448   }
 449   factorz= TMath::Min(factorz,float(4.));
 450   errz*=factorz;
 451
 452   erry= TMath::Min(erry,float(0.05));
 453   errz= TMath::Min(errz,float(0.05));
 454   return 200;
 455 }
 456
 457
 458 AliITSclusterV2 * AliITSClusterErrAnal::FindNearestCluster(AliITSCI * clinfo,  Float_t dmax)
 459 {
 460   //
 461   //
 462   AliTrackReference * ref = &(clinfo->fRef);
 463   TObjArray * clusters = (TObjArray*)(fClusters->At(ref->GetTrack()));
 464   if (!clusters) return 0;
 465   Int_t nclusters = clusters->GetEntriesFast();
 466   Float_t maxd = dmax;
 467   AliITSclusterV2 * cluster=0;
 468   Float_t radius = ref->R();
 469   clinfo->fNClusters=0;
 470   for (Int_t icluster=0;icluster<nclusters;icluster++){
 471     AliITSclusterV2 * cl = (AliITSclusterV2*)clusters->At(icluster);
 472     if (!cl) continue;
 473     Float_t dz = cl->GetZ()-ref->Z();
 474     //
 475     if (TMath::Abs(dz)>dmax) continue;
 476     Int_t idet=cl->GetDetectorIndex();
 477     Double_t r=-1;
 478     Float_t  dy =0;
 479     //    for (Int_t i=cl->fLayer;i<=cl->fLayer-1;i++){
 480
 481     AliITStrackerMI::AliITSdetector & det = fTracker->GetDetector(cl->fLayer,idet);
 482       if (TMath::Abs(radius-det.GetR())<1.5) {
 483         Bool_t isOK=kFALSE;
 484         for (Int_t itrack=0;itrack<3;itrack++){
 485           if (cl->fTracks[itrack]== clinfo->fRef.GetTrack()) {
 486             clinfo->fNClusters++;
 487             isOK= kTRUE;
 488           }
 489         }
 490         if (!isOK) continue;
 491         r = det.GetR();
 492         Double_t phicl= det.GetPhi();
 493         Double_t ca=TMath::Cos(phicl), sa=TMath::Sin(phicl);
 494         Double_t x =  double(ref->X())*ca +double(ref->Y())*sa;
 495         Double_t y = -double(ref->X())*sa +double(ref->Y())*ca;
 496         if (TMath::Abs(x-r-0.012)>0.1) continue;
 497         //
 498         Double_t px = ref->Px(), py = ref->Py();
 499         Double_t lpx =  double(px)*ca +double(py)*sa;
 500         Double_t lpy =  -double(px)*sa +double(py)*ca;
 501         Double_t lphi =  TMath::ATan2(lpy,lpx);
 502
 503         Double_t theta  = ref->Pz()/ref->Pt();
 504         Double_t dx = x-r;
 505         dy = y       -cl->GetY()-dx*lphi;
 506         dz = ref->Z()-cl->GetZ()-dx*theta;
 507         //dx-=0.012;
 508         Double_t d = TMath::Sqrt(dy*dy+dz*dz);
 509         if (d<maxd){
 510           cluster=cl;
 511           maxd=d;
 512           clinfo->Update(fTracker);
 513           clinfo->fCl = *cluster;
 514           clinfo->fLx    = x;
 515           clinfo->fLy    = y;
 516           clinfo->fLayer = cl->fLayer;
 517           clinfo->fPhiCl = phicl;
 518           clinfo->fRCl   = det.GetR();
 519           clinfo->fPhiL  = lphi;
 520           clinfo->fDz    = dz;
 521           clinfo->fDy    = dy;
 522           cl->SetSigmaZ2(TMath::Abs(cl->GetSigmaZ2()));
 523           cl->SetSigmaY2(TMath::Abs(cl->GetSigmaY2()));
 524           clinfo->fPoolZ  = dz/TMath::Sqrt(cl->GetSigmaZ2());
 525           clinfo->fPoolY  = dy/TMath::Sqrt(cl->GetSigmaY2());
 526           clinfo->fPoolY2  = dy/TMath::Sqrt(cl->GetSigmaY2());
 527           clinfo->fErrY   = TMath::Sqrt(cl->GetSigmaY2());
 528           clinfo->fErrZ   = TMath::Sqrt(cl->GetSigmaZ2());
 529           Float_t erry,errz,ny,nz;
 530           clinfo->fErrType = GetError(cl->fLayer,cl,theta,lphi,clinfo->fExpectedQ,erry,errz);
 531           GetNTeor(cl->fLayer, theta,lphi,ny,nz);
 532           clinfo->fNormQ = cl->GetQ()/TMath::Sqrt(1+theta*theta+lpy*lpy/(lpx*lpx));
 533           clinfo->fTNy  = ny;
 534           clinfo->fTNz  = nz;
 535           clinfo->fErrY = erry;
 536           clinfo->fErrZ = errz;
 537           clinfo->fPoolY2 = dy/clinfo->fErrY;
 538           clinfo->fPoolZ2 = dz/clinfo->fErrZ;
 539           clinfo->fLabPos =-1;
 540           Double_t xyz[3];
 541           det.GetGlobalXYZ(cl,xyz);
 542           clinfo->fGx = xyz[0];
 543           clinfo->fGy = xyz[1];
 544           for (Int_t i=0;i<3;i++) if (cl->GetLabel(i)==clinfo->fRef.GetTrack()) clinfo->fLabPos = i;
 545         }
 546       }
 547     }
 548   //  }
 549   return cluster;
 550 }
 551
 552
 553
 554 Int_t AliITSClusterErrAnal::Analyze(Int_t trackmax) {
 555   //
 556   //
 557   // dummy cluster to be fill if not cluster info
 558   //
 559   AliITSCI * clinfo = new AliITSCI;
 560   AliITSclusterV2 dummy;
 561   TBranch *  branch = fTreeA->GetBranch("itscl");
 562   branch->SetAddress(&clinfo);
 563   Int_t nall =0;
 564   Int_t withcl =0;
 565   Int_t ntracks = fReferences->GetEntriesFast();
 566   ntracks = TMath::Min(trackmax,ntracks);
 567   for (Int_t itrack=0;itrack<ntracks;itrack++){
 568     TClonesArray * references = (TClonesArray*) fReferences->At(itrack);
 569     if (!references) continue;
 570     Int_t nreferences = references->GetEntriesFast();
 571     if (nreferences<4) continue;
 572     TObjArray * clarray = (TObjArray*)fClusters->At(itrack);
 573     for (Int_t iref=0;iref<nreferences;iref++){
 574       AliTrackReference* trackref = (AliTrackReference*)references->At(iref);
 575       if (trackref&& trackref->P()>0.01){
 576         TParticle * p = fStack->Particle(trackref->GetTrack());
 577         if (p) clinfo->fP = *p;
 578         clinfo->fRef = *trackref;
 579         clinfo->fStatus = 0;
 580         clinfo->Update(fTracker);
 581         nall++;
 582         AliITSclusterV2 * cluster = FindNearestCluster(clinfo,3.0);
 583         if (cluster){
 584           //clinfo->fCl = *cluster;
 585           clinfo->fStatus=1;
 586           withcl++;
 587         }else{
 588           clinfo->fCl= dummy;
 589         }
 590         clinfo->Update(fTracker);
 591         fTreeA->Fill();
 592       }
 593     }
 594   }
 595
 596   fFileA->cd();
 597   fTreeA->Write();
 598   fTreeB->Write();
 599   fFileA->Close();
 600   return 0;
 601 }
 602
 603
 604
 605 void AliITSClusterErrAnal::MakeSeeds(Double_t zv)
 606 {
 607   //
 608   //
 609   //
 610
 611   AliITStrackerMI::AliITSlayer & layer3 =   fTracker->GetLayer(3);
 612   AliITStrackerMI::AliITSlayer & layer2 =   fTracker->GetLayer(2);
 613   AliITStrackerMI::AliITSlayer & layer1 =   fTracker->GetLayer(1);
 614   AliITStrackerMI::AliITSlayer & layer0 =   fTracker->GetLayer(0);
 615   Int_t ncl3 =  layer3.GetNumberOfClusters();
 616   Int_t ncl2 =  layer2.GetNumberOfClusters();
 617   Int_t ncl1 =  layer1.GetNumberOfClusters();
 618   Int_t ncl0 =  layer0.GetNumberOfClusters();
 619   TH1F *hc = new TH1F("hc","hc",100,-0.01,0.01);
 620   TH1F *hrc = new TH1F("hrc","hrc",100,-1,1);
 621   TH1F *hl = new TH1F("hl","hl",100,-0.05,0.05);
 622   TH1F *hlr = new TH1F("hlr","hlr",1000,-30.,30.);
 623   TH2F *hcl = new TH2F("hcl","hcl",1000,-0.01,0.01,1000,-0.03,0.03);
 624   TH2F *hcc = new TH2F("hcc","hcc",1000,-0.01,0.01,1000,-0.03,0.03);
 625   TH1F *hr = new TH1F("hr","hr",100,0,2000.);
 626   //
 627   TH1F *hpoolc = new TH1F("hpoolc","hpoolc",100,-5,5);
 628   TH1F *hpooll = new TH1F("hpooll","hpooll",100,-5,5);
 629   //
 630   Double_t ratio23 = layer2.GetR()/layer3.GetR();
 631   Double_t ratio12 = layer1.GetR()/layer2.GetR();
 632   Double_t referenceR =  (layer0.GetR()+layer1.GetR()+layer2.GetR()+layer3.GetR())*0.25;
 633   //
 634   Double_t *clusterxyzr3 = new Double_t[ncl3*6];
 635   Double_t *clusterxyzr2 = new Double_t[ncl2*6];
 636   Double_t *clusterxyzr1 = new Double_t[ncl1*6];
 637   Double_t *clusterxyzr0 = new Double_t[ncl0*6];
 638   //
 639   // Get global coordinates 3
 640   for (Int_t icl3=0;icl3<ncl3;icl3++){
 641     AliITSclusterV2* cl3 = layer3.GetCluster(icl3);
 642     Double_t *xyz3 = &clusterxyzr3[6*icl3];
 643     AliITStrackerMI::AliITSdetector & det = layer3.GetDetector(cl3->GetDetectorIndex());
 644     det.GetGlobalXYZ(cl3,xyz3);
 645     xyz3[3] = TMath::Sqrt(xyz3[0]*xyz3[0]+xyz3[1]*xyz3[1]);
 646     Double_t fi3 = TMath::ATan2(xyz3[1],xyz3[0]);
 647     xyz3[4] = TMath::Cos(fi3);
 648     xyz3[5] = TMath::Sin(fi3);
 649   }
 650   //
 651   // Get global coordinates 2
 652   for (Int_t icl2=0;icl2<ncl2;icl2++){
 653     AliITSclusterV2* cl2 = layer2.GetCluster(icl2);
 654     Double_t *xyz2 = &clusterxyzr2[6*icl2];
 655     AliITStrackerMI::AliITSdetector & det = layer2.GetDetector(cl2->GetDetectorIndex());
 656     det.GetGlobalXYZ(cl2,xyz2);
 657     xyz2[3] = TMath::Sqrt(xyz2[0]*xyz2[0]+xyz2[1]*xyz2[1]);
 658     Double_t fi2 = TMath::ATan2(xyz2[1],xyz2[0]);
 659     xyz2[4] = TMath::Cos(fi2);
 660     xyz2[5] = TMath::Sin(fi2);
 661   }
 662   //
 663   // Get global coordinates 1
 664   for (Int_t icl1=0;icl1<ncl1;icl1++){
 665     AliITSclusterV2* cl1 = layer1.GetCluster(icl1);
 666     Double_t *xyz1 = &clusterxyzr1[6*icl1];
 667     AliITStrackerMI::AliITSdetector & det = layer1.GetDetector(cl1->GetDetectorIndex());
 668     det.GetGlobalXYZ(cl1,xyz1);
 669     xyz1[3] = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
 670     Double_t fi1 = TMath::ATan2(xyz1[1],xyz1[0]);
 671     xyz1[4] = TMath::Cos(fi1);
 672     xyz1[5] = TMath::Sin(fi1);
 673   }
 674   //
 675   // Get global coordinates 0
 676   for (Int_t icl0=0;icl0<ncl0;icl0++){
 677     AliITSclusterV2* cl0 = layer0.GetCluster(icl0);
 678     Double_t *xyz0 = &clusterxyzr0[6*icl0];
 679     AliITStrackerMI::AliITSdetector & det = layer0.GetDetector(cl0->GetDetectorIndex());
 680     det.GetGlobalXYZ(cl0,xyz0);
 681     xyz0[3] = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
 682     Double_t fi0 = TMath::ATan2(xyz0[1],xyz0[0]);
 683     xyz0[4] = TMath::Cos(fi0);
 684     xyz0[5] = TMath::Sin(fi0);
 685
 686   }
 687   //
 688   //
 689   Int_t ntracks=0;
 690   Int_t ntracks0=0;
 691   Int_t good0=0;
 692   Int_t accepted0=0;
 693   //
 694   Int_t nfoundtracks=0;
 695   Int_t good1 =0;;
 696   Int_t accepted1=0;
 697   Int_t accepted2=0;
 698   //
 699   Int_t *trackindex = new Int_t[ncl3*4*300];
 700   const AliITSclusterV2 **trackclusters = new const AliITSclusterV2*[ncl3*4*300];
 701   //
 702   for (Int_t icl3=0;icl3<ncl3;icl3++){
 703     AliITSclusterV2* cl3 = layer3.GetCluster(icl3);
 704     Double_t *xyz3 = &clusterxyzr3[icl3*6];
 705     Double_t &r3   = clusterxyzr3[icl3*6+3];
 706     ntracks=0;  //number of possible tracks corresponding to given cluster
 707     //
 708     //
 709     // find pairs in layer 2
 710     Float_t zc   = (xyz3[2]-zv)*ratio23+zv;
 711     Float_t lphi  = TMath::ATan2(xyz3[1],xyz3[0]);
 712     layer2.SelectClusters(zc-1,zc+1, layer2.GetR()*lphi-1, layer2.GetR()*lphi+1.);
 713     const AliITSclusterV2 * cl2=0;
 714     Int_t npairs2=0;
 715     Int_t ci2[100];
 716     while ( (cl2=layer2.GetNextCluster(ci2[npairs2]))!=0){
 717       Int_t index2 = ci2[npairs2];
 718       Double_t *xyz2 = &clusterxyzr2[index2*6];
 719       Double_t &r2   = clusterxyzr2[index2*6+3];
 720       //
 721       if (TMath::Abs((xyz3[2]-zv)*r2/r3-(xyz2[2]-zv))>0.3) continue;
 722       npairs2++;
 723     }
 724     //
 725     //
 726     // find pairs in layer 2
 727     for (Int_t jcl2 =0;jcl2<npairs2; jcl2++){
 728       Int_t index2 = ci2[jcl2];
 729       Double_t *xyz2 = &clusterxyzr2[index2*6];
 730       Float_t zc   = (xyz2[2]-zv)*ratio12+zv;
 731       Float_t lphi  = TMath::ATan2(xyz2[1],xyz2[0]);
 732       cl2 = layer2.GetCluster(index2);
 733       //
 734       //
 735       layer1.SelectClusters(zc-1,zc+1, layer1.GetR()*lphi-1, layer1.GetR()*lphi+1.);
 736       const AliITSclusterV2 * cl1=0;
 737       Int_t ci1;
 738       //
 739       Double_t lr2 = (xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1]);
 740       Double_t lr  = TMath::Sqrt(lr2);
 741       Double_t sinphi = (xyz3[1]-xyz2[1])/lr;
 742       Double_t cosphi = (xyz3[0]-xyz2[0])/lr;
 743       Double_t t1 = (xyz3[2]-xyz2[2])/lr;
 744       //
 745       while ( (cl1=layer1.GetNextCluster(ci1))!=0){
 746         if (cl1->GetQ()==0) continue;
 747         Int_t index1 = ci1;
 748         Double_t *xyz1 = &clusterxyzr1[index1*6];
 749         AliITSclusterV2* cl1 = layer1.GetCluster(index1);
 750
 751         if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) good0++;
 752         Double_t loc1[2];
 753         //
 754         loc1[0] =  (xyz1[0]-xyz2[0])*cosphi+(xyz1[1]-xyz2[1])*sinphi;
 755         loc1[1] = -(xyz1[0]-xyz2[0])*sinphi+(xyz1[1]-xyz2[1])*cosphi;
 756         //
 757         Double_t c1 = 2.*loc1[1]/(loc1[0]*(loc1[0]-lr));   //curvature approximated using 2 derivation
 758         Double_t c=c1;
 759         Double_t t2 = (xyz2[2]-xyz1[2])/TMath::Sqrt(loc1[0]*loc1[0]+loc1[1]*loc1[1]);
 760         if (TMath::Abs(c)>0.008) continue;
 761         if (TMath::Abs(t2-t1)>0.02) continue;
 762         //
 763         trackclusters[ntracks*4]    = 0;
 764         trackclusters[ntracks*4+1]  = cl1;
 765         trackclusters[ntracks*4+2]  = cl2;
 766         trackclusters[ntracks*4+3]  = cl3;
 767         trackindex[ntracks*4+0] = 0;
 768         trackindex[ntracks*4+1] = index1;
 769         trackindex[ntracks*4+2] = index2;
 770         trackindex[ntracks*4+3] = icl3;
 771         if (ntracks>=300*ncl3){
 772           break;
 773         }
 774         if (cl3->GetLabel(0)==cl2->GetLabel(0) && cl3->GetLabel(0)==cl1->GetLabel(0)) {
 775           //printf("%f\n",c);
 776
 777           //hc->Fill(c);
 778           //hl->Fill(t2-t1);
 779           //hr->Fill(1/c);
 780           //hcc->Fill(c,TMath::Abs(c1));
 781           //hcl->Fill(c1,t2-t1);
 782
 783           accepted0++;
 784         }
 785         ntracks++;
 786       }
 787     }   // end - find pairs in layer 2
 788     //
 789     ntracks0+=ntracks;
 790     //hcc->Draw();
 791     Int_t   ntracks2=0;
 792     Int_t   cindex2[10000*4];
 793     Float_t pools[10000];
 794
 795     for (Int_t itrack=0;itrack<ntracks;itrack++){
 796       const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3;
 797       cl3 = trackclusters[itrack*4+3];
 798       cl2 = trackclusters[itrack*4+2];
 799       cl1 = trackclusters[itrack*4+1];
 800       Int_t index1 = trackindex[itrack*4+1];
 801       Int_t index2 = trackindex[itrack*4+2];
 802       Int_t index3 = trackindex[itrack*4+3];
 803       //
 804       Double_t *xyz1 = &clusterxyzr1[index1*6],*xyz2 = &clusterxyzr2[index2*6],*xyz3=&clusterxyzr3[index3*6];
 805       Double_t &r1 = clusterxyzr1[index1*6+3],&r2 = clusterxyzr2[index2*6+3], &r3 = clusterxyzr3[index3*6+3];
 806
 807       Double_t r0 = layer0.GetR();
 808       Float_t zc   = xyz1[2] + (xyz2[2]-xyz1[2])*(r0-r1)/(r2-r1);
 809       Float_t lphi  = TMath::ATan2(xyz1[1],xyz1[0]);
 810
 811       layer0.SelectClusters(zc-1,zc+1, layer0.GetR()*lphi-0.3, layer0.GetR()*lphi+0.3);
 812       Int_t ci0;
 813       //
 814       //
 815       Double_t lr2 = (xyz2[0]-xyz1[0])*(xyz2[0]-xyz1[0])+(xyz2[1]-xyz1[1])*(xyz2[1]-xyz1[1]);
 816       Double_t lr  = TMath::Sqrt(lr2);
 817       Double_t sinphi = (xyz2[1]-xyz1[1])/lr;
 818       Double_t cosphi = (xyz2[0]-xyz1[0])/lr;
 819       //
 820       Double_t loc3[2]; // cluster3 in local coordindate
 821       loc3[0] =  (xyz3[0]-xyz1[0])*cosphi+(xyz3[1]-xyz1[1])*sinphi;
 822       loc3[1] = -(xyz3[0]-xyz1[0])*sinphi+(xyz3[1]-xyz1[1])*cosphi;
 823
 824       Double_t tan1 = (xyz3[2]-xyz2[2])/TMath::Sqrt((xyz3[0]-xyz2[0])*(xyz3[0]-xyz2[0])+(xyz3[1]-xyz2[1])*(xyz3[1]-xyz2[1]));
 825       Double_t c3 = 2.*loc3[1]/(loc3[0]*(loc3[0]-lr));   //curvature approximated using 2 derivation
 826       //
 827       while ( (cl0=layer0.GetNextCluster(ci0))!=0){
 828         if (cl0->GetQ()==0) continue;
 829         Double_t *xyz0 = &clusterxyzr0[ci0*6];
 830         if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0)
 831             && cl3->GetLabel(0)==cl1->GetLabel(0)&&cl3->GetLabel(0)==cl0->GetLabel(0)){
 832           good1++;
 833         }
 834         Double_t loc0[2]; // cluster3 in local coordindate
 835         loc0[0] =  (xyz0[0]-xyz1[0])*cosphi+(xyz0[1]-xyz1[1])*sinphi;
 836         loc0[1] = -(xyz0[0]-xyz1[0])*sinphi+(xyz0[1]-xyz1[1])*cosphi;
 837         //
 838         Double_t c0 = 2.*loc0[1]/(loc0[0]*(loc0[0]-lr));   //curvature approximated using 2 derivation
 839         Double_t tan2 = (xyz1[2]-xyz0[2])/TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+(xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]));
 840         //if (TMath::Abs(tan1-tan2)>0.05) continue;
 841         //if (TMath::Abs(c3-c0)>0.006) continue;
 842         Double_t poolphi = (c3-c0)/0.0011;
 843         poolphi*=poolphi;
 844         Double_t poolr = ((c3-c0)/(TMath::Abs(c3+c0)+0.001))/0.14;
 845         Double_t poolth = (tan1-tan2)/0.0094;
 846         poolth*=poolth;
 847         Double_t poolthr = ((tan1-tan2)/(TMath::Abs(c3+c0)+0.005))/0.83;
 848         if (poolr*poolr+poolthr*poolthr>40) continue;
 849         //if ( TMath::Abs((c3-c0)/(c3+c0))>0.6) continue;
 850         cindex2[ntracks2*4+0]   = ci0;
 851         cindex2[ntracks2*4+1]   = index1;
 852         cindex2[ntracks2*4+2]   = index2;
 853         cindex2[ntracks2*4+3]   = index3;
 854         pools[ntracks2] = poolr*poolr+poolthr*poolthr;
 855         //
 856         ntracks2++;
 857         if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0)
 858             && cl3->GetLabel(0)==cl1->GetLabel(0)){
 859           hl->Fill(tan1-tan2);
 860           hrc->Fill(poolr);
 861           hc->Fill((c3-c0));
 862           hlr->Fill(poolthr);
 863           hpoolc->Fill(poolr);
 864           hpooll->Fill(poolthr);
 865           accepted1++;
 866         }
 867       }
 868     }
 869     Int_t sindex[10000];
 870     TMath::Sort(ntracks2,pools,sindex,kFALSE);
 871     Float_t minpool = pools[sindex[0]];
 872     Int_t max = TMath::Min(ntracks2,20);
 873     for (Int_t itrack=0;itrack<ntracks2;itrack++){
 874       if (itrack>=max) continue;
 875       //if (pools[sindex[itrack]]>minpool+5) continue;
 876       const AliITSclusterV2 *cl0, *cl1,*cl2,*cl3;
 877       cl3 = layer3.GetCluster(cindex2[sindex[itrack]*4+3]);
 878       cl2 = layer2.GetCluster(cindex2[sindex[itrack]*4+2]);
 879       cl1 = layer1.GetCluster(cindex2[sindex[itrack]*4+1]);
 880       cl0 = layer0.GetCluster(cindex2[sindex[itrack]*4+0]);
 881       if (cl3->GetLabel(0)==cl0->GetLabel(0) && cl3->GetLabel(0)==cl2->GetLabel(0)
 882           && cl3->GetLabel(0)==cl1->GetLabel(0)){
 883         accepted2++;
 884       }
 885       nfoundtracks++;
 886     }
 887   }
 888   hc->Draw();
 889   printf("First pass \t%d\t%d\t%d\n",ntracks0,good0,accepted0);
 890   printf("Found tracks\t%d\t%d\t%d\t%d\n",nfoundtracks,good1, accepted1,accepted2);
 891
 892 }