Include stdlib.h needed on HP
[u/mrichter/AliRoot.git] / ITS / AliITSclusterTable.cxx
index fc14343ad5355914829d0de56a54e006cc728bbd..c2bc0eafe6dcf1afe683dcc0d60156053018d59b 100644 (file)
 // -Function FillArrayLabel fills an array wich contains, for each       //
 //  particle label, and for each layer, the information on clusters:     //
 //  0 if there is no cluster, 1 if there is a cluster with this label.   //
-//  This function is used to define trackable tracks.                    //   
+//  This function is used to define trackable tracks.                    //
+// -Function FillArrayCoorAngles fills arrays wich contains, for each    //
+//  layer, the global coordinates, errors on x,y,z and angles lambda     //
+//  and phi for each cluster                                             //
 ///////////////////////////////////////////////////////////////////////////
 
 
+#include <stdlib.h>
 #include "AliITSclusterTable.h"
 #include "AliITSclusterV2.h"
 #include "AliITSgeom.h"
@@ -35,7 +39,6 @@
 #include<TClonesArray.h>
 #include<TTree.h>
 ClassImp(AliITSclusterTable)
-
 //__________________________________________________________
 AliITSclusterTable::AliITSclusterTable(){
 // Default constructor
@@ -44,10 +47,19 @@ AliITSclusterTable::AliITSclusterTable(){
   fLbl=0;
   fGeom = 0;
   fTracker = 0;
+  fPhiList = 0; 
+  fLambdaList = 0;
+  fXList = 0;
+  fYList =0;
+  fZList =0;
+  fSxList = 0;
+  fSyList =0;
+  fSzList =0;
+  for(Int_t i=0;i<3;i++){fPrimaryVertex[i]=0;}
 }
 
 //__________________________________________________________
-AliITSclusterTable::AliITSclusterTable(AliITSgeom* geom, AliITStrackerSA* tracker) {
+AliITSclusterTable::AliITSclusterTable(AliITSgeom* geom, AliITStrackerSA* tracker,Double_t* primaryVertex) {
 // Standard constructor
   fDet=0;
   Int_t nm = geom->GetIndexMax();
@@ -56,8 +68,16 @@ AliITSclusterTable::AliITSclusterTable(AliITSgeom* geom, AliITStrackerSA* tracke
   fLbl=0;
   fGeom = geom;
   fTracker = tracker;
+  for(Int_t i=0;i<3;i++){fPrimaryVertex[i]=primaryVertex[i];}
+  fPhiList = 0;
+  fLambdaList = 0;
+  fXList = 0;
+  fYList = 0;
+  fZList =0;
+  fSxList = 0;
+  fSyList =0;
+  fSzList =0;
 }
-
 //______________________________________________________________________
 AliITSclusterTable::AliITSclusterTable(const AliITSclusterTable &tab) : 
                     TObject(tab) {
@@ -83,31 +103,40 @@ AliITSclusterTable::~AliITSclusterTable(){
     for (Int_t i=0; i<nm; i++){
       delete fDet[i];
     }
-    delete fDet;
+    delete[] fDet;
   }
-  /*  if(fLbl){
-    for (Int_t i=0; i<nm; i++){
-      delete fLbl[i];
-    }
-    delete fLbl;
-    }*/
-  if (fLbl)delete [] fLbl;
+  if (fLbl)delete [] fLbl; // memory leak!
   if (fNCl)delete [] fNCl;
+  for (Int_t i = 0; i < fGeom->GetNlayers(); i++) {
+    if (fPhiList) delete fPhiList[i];
+    if (fLambdaList) delete fLambdaList[i];
+    if (fXList) delete fXList[i];
+    if (fYList) delete fYList[i];
+    if (fZList) delete fZList[i];
+    if (fSxList) delete fSxList[i];
+    if (fSyList) delete fSyList[i];
+    if (fSzList) delete fSzList[i];
+  }
+  if (fPhiList) delete[] fPhiList;
+  if (fLambdaList) delete[] fLambdaList;
+  if (fXList) delete[] fXList;
+  if (fYList) delete[] fYList;
+  if (fZList) delete[] fZList;
+  if (fSxList) delete[] fSxList;
+  if (fSyList) delete[] fSyList;
+  if (fSzList) delete[] fSzList;
 }
 
 //__________________________________________________________
 
-void AliITSclusterTable::FillArray(TTree* clusterTree,Int_t evnumber){
+void AliITSclusterTable::FillArray(TTree* clusterTree){
   
   //
   Int_t nm = fGeom->GetIndexMax();
   fDet = new TArrayI*[nm];
-  fTracker->SetEventNumber(evnumber);
-  fTracker->LoadClusters(clusterTree);
 
   TArrayI** vect = new TArrayI*[fGeom->GetNlayers()];
-  
-  Int_t firstmod[fGeom->GetNlayers()+1];
+  Int_t * firstmod = new Int_t[fGeom->GetNlayers()+1];
   firstmod[fGeom->GetNlayers()]=fGeom->GetIndexMax();  // upper limit
   for(Int_t nlayer=0;nlayer<fGeom->GetNlayers();nlayer++){
     firstmod[nlayer] = fGeom->GetModuleIndex(nlayer+1,1,1);
@@ -124,7 +153,6 @@ void AliITSclusterTable::FillArray(TTree* clusterTree,Int_t evnumber){
   TClonesArray* clus = new TClonesArray("AliITSclusterV2",10000);
   brancht->SetAddress(&clus);
  
-  
   for(Int_t mod=0;mod<nm;mod++){
     Int_t nc=0;
     clusterTree->GetEvent(mod);
@@ -134,35 +162,34 @@ void AliITSclusterTable::FillArray(TTree* clusterTree,Int_t evnumber){
     Int_t nlr = FindIndex(fGeom->GetNlayers(),firstmod,mod);
     if(nlr<0){
       Fatal("FillArray","Wrong module number %d . Limits: %d , %d",mod,firstmod[0],firstmod[fGeom->GetNlayers()+1]);
-      return;
+      exit(1);
     }
     else {
       for(Int_t n=0;n<vect[nlr]->GetSize();n++){
        Int_t mm=vect[nlr]->At(n);
-       if(mm==mod) {fDet[mod]->AddAt(n,nc); nc+=1; }
+       if(nc>=fDet[mod]->GetSize()) fDet[mod]->Set(nc*2+10);
+       if(mm==mod) {(*fDet[mod])[nc]=n; nc+=1; }
       }
     }
   }
 
-
-  fTracker->UnloadClusters();
-
+  clus->Delete();
+  delete clus;
   for(Int_t n=0;n<fGeom->GetNlayers();n++)delete vect[n];
-  delete vect;
+  delete [] vect;
+  delete [] firstmod;
 }
 
 //_________________________________________________________________
-void AliITSclusterTable::FillArrayLabel(const Int_t numberofparticles,TTree* clusterTree,Int_t evnumber){
+void AliITSclusterTable::FillArrayLabel(Int_t numberofparticles){
   //
 
-
   fLbl = new TArrayI*[numberofparticles];
-  fTracker->SetEventNumber(evnumber);
-  fTracker->LoadClusters(clusterTree);
   const Int_t knm =fGeom->GetNlayers();
   for(Int_t nlab=0;nlab<numberofparticles;nlab++){
     fLbl[nlab] = new TArrayI(knm);
-    Int_t nn[knm]; for(Int_t i=0;i<knm;i++)nn[i]=0;
+    Int_t * nn = new Int_t[knm]; 
+    for(Int_t i=0;i<knm;i++)nn[i]=0;
     for(Int_t nlayer=0;nlayer<knm;nlayer++){
       Int_t ncl = fTracker->GetNumberOfClustersLayer(nlayer);
       while(ncl--){
@@ -176,9 +203,9 @@ void AliITSclusterTable::FillArrayLabel(const Int_t numberofparticles,TTree* clu
       }     
       fLbl[nlab]->AddAt(nn[nlayer],nlayer);
     }
+    delete [] nn;
   }
 
-  fTracker->UnloadClusters();
 }
 
 //_______________________________________________________________
@@ -211,3 +238,94 @@ Int_t AliITSclusterTable::FindIndex(Int_t ndim, Int_t *ptr, Int_t value){
   }
   return retval;
 }
+
+void AliITSclusterTable::FillArrayCoorAngles(){
+  //Fill arrays with phi,lambda and indices of clusters for each layer
+  fPhiList = new TArrayD*[fGeom->GetNlayers()];
+  fLambdaList = new TArrayD*[fGeom->GetNlayers()];
+  fXList = new TArrayF*[fGeom->GetNlayers()];
+  fYList = new TArrayF*[fGeom->GetNlayers()];
+  fZList = new TArrayF*[fGeom->GetNlayers()];
+  fSxList = new TArrayF*[fGeom->GetNlayers()];
+  fSyList = new TArrayF*[fGeom->GetNlayers()];
+  fSzList = new TArrayF*[fGeom->GetNlayers()];
+
+  Int_t * firstmod = new Int_t[fGeom->GetNlayers()+1];
+  firstmod[fGeom->GetNlayers()]=fGeom->GetIndexMax();  // upper limit
+
+  for(Int_t nlay=0;nlay<fGeom->GetNlayers();nlay++){
+    firstmod[nlay] = fGeom->GetModuleIndex(nlay+1,1,1);
+    Int_t ncl = fTracker->GetNumberOfClustersLayer(nlay);
+    fPhiList[nlay] = new TArrayD(ncl);
+    fLambdaList[nlay]=new TArrayD(ncl);
+    fXList[nlay]=new TArrayF(ncl);
+    fYList[nlay]=new TArrayF(ncl);
+    fZList[nlay]=new TArrayF(ncl);
+    fSxList[nlay]=new TArrayF(ncl);
+    fSyList[nlay]=new TArrayF(ncl);
+    fSzList[nlay]=new TArrayF(ncl);
+
+    for(Int_t j=0;j<ncl;j++){
+      AliITSclusterV2* cl = fTracker->GetClusterLayer(nlay,j);
+      Double_t phi=0;Double_t lambda=0;
+      Float_t x=0;Float_t y=0;Float_t z=0;
+      Float_t sx=0;Float_t sy=0;Float_t sz=0;
+      Int_t module = cl->GetDetectorIndex()+firstmod[nlay];
+      GetCoorAngles(cl,module,phi,lambda,x,y,z);
+      GetCoorErrors(cl,module,sx,sy,sz);
+      fPhiList[nlay]->AddAt(phi,j);
+      fLambdaList[nlay]->AddAt(lambda,j);
+      fXList[nlay]->AddAt(x,j);
+      fYList[nlay]->AddAt(y,j);
+      fZList[nlay]->AddAt(z,j);
+      fSxList[nlay]->AddAt(sx,j);
+      fSyList[nlay]->AddAt(sy,j);
+      fSzList[nlay]->AddAt(sz,j);
+      
+    }
+    
+  }
+  
+  delete [] firstmod;
+}
+void AliITSclusterTable::GetCoorAngles(AliITSclusterV2* cl,Int_t module,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z){
+  //Returns values of phi (azimuthal) and lambda angles for a given cluster
+  
+  Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
+  Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
+  Float_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);     
+
+  Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
+  Double_t phi1=TMath::Pi()/2+alpha;
+  if (lay==1) phi1+=TMath::Pi();
+
+  Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
+  Float_t r=tx*cp+ty*sp;
+
+  x= r*cp - cl->GetY()*sp;
+  y= r*sp + cl->GetY()*cp;
+  z=cl->GetZ();
+  
+  phi=TMath::ATan2(y,x);
+  lambda=TMath::ATan2(z-fPrimaryVertex[2],TMath::Sqrt((x-fPrimaryVertex[0])*(x-fPrimaryVertex[0])+(y-fPrimaryVertex[1])*(y-fPrimaryVertex[1])));
+}
+
+void AliITSclusterTable::GetCoorErrors(AliITSclusterV2* cl, Int_t module,Float_t &sx,Float_t &sy, Float_t &sz){
+
+  //returns x,y,z of cluster in global coordinates
+
+  Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
+  Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
+  Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
+  Double_t phi=TMath::Pi()/2+alpha;
+  if (lay==1) phi+=TMath::Pi();
+
+  Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
+
+  sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
+  sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
+  sz = TMath::Sqrt(cl->GetSigmaZ2());
+
+}