AliTPCcalibBase.h AliTPCcalibBase.cxx
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Mar 2009 19:03:46 +0000 (19:03 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Mar 2009 19:03:46 +0000 (19:03 +0000)
Adding trigger class as string

AliTPCPointCorrection.cxx AliTPCPointCorrection.h
Adding
R            distortion correction
RPhi         distortion correction
Quadrant  misalignment
Inter sector misalignment represented by 6 parameter to be still defined.

AliTPCcalibCosmic.cxx  AliTPCcalibCosmic.h
Adding MakeCombinedTrack       - merge 2 halves of the track
Removing setting of Gain map   - done automatically in the CookdEdx from the OCDB
Adding AliExternalComparison   - ndim THnsparse for performance study

AliTPCcalibTracks.cxx
Changing binning of QA histograms

AliTPCcalibCalib.cxx
Applying the R,RPhi and alignment correction before refit of the track
It is temporary sollution before we move the all correction to the AliTPCTransform
Corresponding OCDB  entry should be back compatible therefore I didn't want to commit it before correction class (AliTPCPointCorrection) is fixed

TPC/AliTPCPointCorrection.cxx
TPC/AliTPCPointCorrection.h
TPC/AliTPCcalibBase.cxx
TPC/AliTPCcalibBase.h
TPC/AliTPCcalibCalib.cxx
TPC/AliTPCcalibCosmic.cxx
TPC/AliTPCcalibCosmic.h
TPC/AliTPCcalibTracks.cxx

index 1b0fb3e..7541e64 100644 (file)
@@ -39,6 +39,7 @@
 #include "TVectorD.h"
 #include "AliLog.h"
 #include "AliTPCROC.h"
+#include "AliTPCClusterParam.h"
 #include "AliTPCPointCorrection.h"
 
 AliTPCPointCorrection* AliTPCPointCorrection::fgInstance = 0;
@@ -52,11 +53,27 @@ AliTPCPointCorrection::AliTPCPointCorrection():
   fParamOutRVersion(0),
   fErrorsOutR(38),
   fErrorsOutZ(38),
-  fParamOutZVersion(0)
+  fParamOutZVersion(0),
+  //
+  fXIO(0),
+  fXmiddle(0),
+  fXquadrant(0),
+  fArraySectorIntParam(36), // array of sector alignment parameters
+  fArraySectorIntCovar(36), // array of sector alignment covariances 
+  //
+  // Kalman filter for global alignment
+  //
+  fSectorParam(0),     // Kalman parameter   
+  fSectorCovar(0)     // Kalman covariance  
+
 {
   //
   // Default constructor
   //
+  AliTPCROC * roc = AliTPCROC::Instance();
+  fXquadrant = roc->GetPadRowRadii(36,53);
+  fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+  fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
 }
 
 AliTPCPointCorrection::AliTPCPointCorrection(const Text_t *name, const Text_t *title):
@@ -66,11 +83,30 @@ AliTPCPointCorrection::AliTPCPointCorrection(const Text_t *name, const Text_t *t
   fParamOutRVersion(0),
   fErrorsOutR(38),
   fErrorsOutZ(38),
-  fParamOutZVersion(0)
+  fParamOutZVersion(0),
+  //
+  // 
+  //
+  fXIO(0),
+  fXmiddle(0),
+  fXquadrant(0),
+  fArraySectorIntParam(36), // array of sector alignment parameters
+  fArraySectorIntCovar(36), // array of sector alignment covariances 
+  //
+  // Kalman filter for global alignment
+  //
+  fSectorParam(0),     // Kalman parameter   for A side
+  fSectorCovar(0)     // Kalman covariance  for A side 
+  
 {
   //
   // Non default constructor
   //
+  AliTPCROC * roc = AliTPCROC::Instance();
+  fXquadrant = roc->GetPadRowRadii(36,53);
+  fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+  fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
+
 }
 
 AliTPCPointCorrection::~AliTPCPointCorrection(){
@@ -151,7 +187,8 @@ Double_t AliTPCPointCorrection::CorrectionOutR0(Bool_t isGlobal, Bool_t type,  D
   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);  
   AliTPCROC * roc = AliTPCROC::Instance();
   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
-  Double_t dout = xout-radius;
+  Double_t dout = xout-radius;  
+  if (dout<0) return 0;
   //drift
   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
   //
@@ -171,9 +208,65 @@ Double_t AliTPCPointCorrection::CorrectionOutR0(Bool_t isGlobal, Bool_t type,  D
   result+=(*vec)[5]*eoutp*dr*dr;
   result+=(*vec)[6]*eoutm*dr*dr;
   return result;
+}
 
+Double_t AliTPCPointCorrection::RPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
+  //
+  // Calculates COG corection in RPHI direction
+  // cluster and track position  y is supposed to be corrected before for other effects
+  // (e.g ExB and alignemnt)
+  // Rphi distortion dependeds on the distance to the edge-pad, distance to the wire edge and
+  // relative distance to the center of the pad. Therefore the y position is trnsfromed to the 
+  // pad coordiante frame (correction offset (ExB alignemnt) substracted). 
+  //   
+  // Input parameters:
+  //
+  // sector - sector number - 0-71  - cl.GetDetector()
+  // padrow - padrow number -0-63 -IROC 0-95 OROC - cl->GetRow()
+  // pad    - mean pad number  - cl->GetPad()
+  // cy     - cluster y        - cl->GetY()
+  //  y     - track -or cluster y
+  //  qMax  - cluster max charge - cl-.GetMax()
+  //  threshold - clusterer threshold
+  //
+  AliTPCClusterParam * clparam = AliTPCClusterParam::Instance(); 
+  Int_t padtype=0;
+  if (sector>=36) padtype = (padrow>64)?2:1;
+  Double_t padwidth=(padtype==0)? 0.4:0.6;
+
+  Double_t sigma = clparam->GetRMS0(0,padtype,250-TMath::Abs(z),ky)/padwidth;
+  //
+  Int_t   npads   =  AliTPCROC::Instance()->GetNPads(sector,padrow);
+  Float_t yshift  =  TMath::Abs(cy)-TMath::Abs((-npads*0.5+pad)*padwidth);   // the clusters can be corrected before
+                                            // shift to undo correction
+  
+  y -= yshift*((y>0)?1.:-1.);                             // y in the sector coordinate
+  Double_t y0     = npads*0.5*padwidth;
+  Double_t dy     = (TMath::Abs(y0)-TMath::Abs(y))/padwidth-0.5;  // distance to  the center of pad0   
+  //
+  Double_t padcenter = TMath::Nint(dy);
+  Double_t sumw=0;
+  Double_t sumwi=0;
+  for (Double_t ip=padcenter-2.5; ip<=padcenter+2.5;ip++){
+    Double_t weightGaus = qMax*TMath::Exp(-(dy-ip)*(dy-ip)/(2*sigma*sigma));
+    Double_t ypad       = (ip-npads*0.5)*padwidth;
+    Double_t weightGain = GetEdgeQ0(sector,padrow,ypad);
+    //
+    Double_t weight = TMath::Nint(weightGaus*weightGain);
+    if (ip<0 &&weight> threshold) weight = threshold;  // this is done in cl finder
+    if (weight<0) continue;
+    sumw+=weight;
+    sumwi+=weight*(ip);    
+  }
+  Double_t result =0;
+  if (sumw>0) result = sumwi/sumw;
+  result = (result-dy)*padwidth;
+  return result;
 }
 
+
+
+
 Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  Double_t cx, Double_t cy, Double_t cz, Int_t sector){
   //
   // return dR correction - for correction version 0 
@@ -195,6 +288,7 @@ Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  D
   AliTPCROC * roc = AliTPCROC::Instance();
   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
   Double_t dout = xout-radius;
+  if (dout<0) return 0;
   //drift
   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
   //
@@ -217,5 +311,162 @@ Double_t AliTPCPointCorrection::CorrectionOutZ0(Bool_t isGlobal, Bool_t type,  D
 
 }
 
+Double_t  AliTPCPointCorrection::GetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
+  //
+  /* TF1 fexp("fexp","1-exp(-[0]*(x-[1]))",0,20)
+          | param [0] | param [1]
+     IROC | 4.71413   | 1.39558
+     OROC1| 2.11437   | 1.52643
+     OROC2| 2.15082   | 1.53537 
+  */
+
+  Double_t params[2]={100,0};
+  if (sector<36){
+    params[0]=4.71413; params[1]=1.39558;
+  }
+  if (sector>=36){
+    if (padrow<64) {  params[0]=2.114; params[1]=1.526;}
+    if (padrow>=64){  params[0]=2.15; params[1]=1.535;}
+  }
+  Double_t result= 1;
+  Double_t xrow  = AliTPCROC::Instance()->GetPadRowRadii(sector,padrow);
+  Double_t ymax  = TMath::Tan(TMath::Pi()/18.)*xrow;
+  Double_t dedge = ymax-TMath::Abs(y);
+  if (dedge-params[2]<0)             return 0;
+  if (dedge>10.*params[1]) return 1;
+  result = 1.-TMath::Exp(-params[0]*(dedge-params[2]));
+  return result;
+}
+
+Double_t AliTPCPointCorrection::SRPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,Float_t qMax, Float_t threshold){
+  return fgInstance->RPhiCOGCorrection(sector, padrow, pad, cy, y,z, ky, qMax, threshold);
+}
+
+Double_t AliTPCPointCorrection::SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y){
+  //
+  //
+  return fgInstance->GetEdgeQ0(sector, padrow, y);
+}
+
+void     AliTPCPointCorrection::AddCorrectionSector(TObjArray & sideAPar, TObjArray &sideCPar, TObjArray & sideACov, TObjArray &sideCCov, Bool_t reset){
+  //
+  //
+  //
+  for (Int_t isec=0;isec<36;isec++){
+    TMatrixD * mat1 = (TMatrixD*)(sideAPar.At(isec));
+    TMatrixD * mat1Covar = (TMatrixD*)(sideCPar.At(isec));
+    if (!mat1) continue;
+    TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
+    TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
+    if (!mat0) {
+      fArraySectorIntParam.AddAt(mat1->Clone(),isec); 
+      fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec); 
+      continue;
+    }
+    (*mat0Covar)=(*mat1Covar);      
+    if (reset)   (*mat0)=(*mat1);
+    if (!reset) (*mat0)+=(*mat1);
+  }
+
+  for (Int_t isec=0;isec<36;isec++){
+    TMatrixD * mat1 = (TMatrixD*)(sideCPar.At(isec));
+    TMatrixD * mat1Covar = (TMatrixD*)(sideCCov.At(isec));
+    if (!mat1) continue;
+    TMatrixD * mat0 = (TMatrixD*)(fArraySectorIntParam.At(isec));
+    TMatrixD * mat0Covar = (TMatrixD*)(fArraySectorIntCovar.At(isec));
+    if (!mat0) {
+      fArraySectorIntParam.AddAt(mat1->Clone(),isec); 
+      fArraySectorIntCovar.AddAt(mat1Covar->Clone(),isec); 
+      continue;
+    }
+    (*mat0Covar)=(*mat1Covar);      
+    if (reset)   (*mat0)=(*mat1);
+    if (!reset) (*mat0)+=(*mat1);
+  }
+}
+
+
+Double_t AliTPCPointCorrection::GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant){
+  //
+  // Get position correction for given sector
+  //
+
+  TMatrixD * param = (TMatrixD*)fArraySectorIntParam.At(sector%36);
+  if (!param) return 0;
+  if (quadrant<0){   //recalc quadrant if not specified
+    if (lx<fXIO) quadrant=0;
+    if(lx>fXIO){
+      if (lx<fXquadrant) {
+       if (ly<0) quadrant=1;
+       if (ly>0) quadrant=2;
+      }
+      if (lx>=fXquadrant) {
+       if (ly<0) quadrant=3;
+       if (ly>0) quadrant=4;
+      }
+    }    
+  }
+  Double_t a10 = (*param)(quadrant*6+0,0);
+  Double_t a20 = (*param)(quadrant*6+1,0);
+  Double_t a21 = (*param)(quadrant*6+2,0);
+  Double_t dx  = (*param)(quadrant*6+3,0);
+  Double_t dy  = (*param)(quadrant*6+4,0);
+  Double_t dz  = (*param)(quadrant*6+5,0);
+  Double_t deltaX = lx-fXmiddle;
+  if (coord==0) return dx;
+  if (coord==1) return dy+deltaX*a10;
+  if (coord==2) return dz+deltaX*a20+ly*a21;
+  if (coord==3) return a10;
+  if (coord==4) return a20;
+  if (coord==5) return a21;
+  //
+  return 0;
+}
+
+Double_t AliTPCPointCorrection::SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant){
+  //
+  //
+  //
+  if (!Instance()) return 0;
+  return Instance()->GetCorrectionSector(coord,sector,lx,ly,lz, quadrant);
+}
+
+
+
+Double_t AliTPCPointCorrection::GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
+  //
+  // Get position correction for given sector
+  //
+  if (!fSectorParam) return 0;
+  
+  Double_t a10 = (*fSectorParam)(sector*6+0,0);
+  Double_t a20 = (*fSectorParam)(sector*6+1,0);
+  Double_t a21 = (*fSectorParam)(sector*6+2,0);
+  Double_t dx  = (*fSectorParam)(sector*6+3,0);
+  Double_t dy  = (*fSectorParam)(sector*6+4,0);
+  Double_t dz  = (*fSectorParam)(sector*6+5,0);
+  Double_t deltaX = lx-fXIO;
+  //
+  if (coord==0) return dx;
+  if (coord==1) return dy+deltaX*a10;
+  if (coord==2) return dz+deltaX*a20+ly*a21;
+  if (coord==3) return a10;
+  if (coord==4) return a20;
+  if (coord==5) return a21;
+  return 0;
+}
+
+Double_t AliTPCPointCorrection::SGetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz){
+  //
+  //
+  //
+  if (!Instance()) return 0;
+  return Instance()->GetCorrection(coord,sector,lx,ly,lz);
+}
+
+
+
+
+
 
 
index cf7c510..8a54efa 100644 (file)
@@ -9,7 +9,6 @@
 #include "TObjArray.h"
 #include "TVectorD.h"
 
-
  
 class AliTPCPointCorrection:public TNamed {
 public:
@@ -32,6 +31,22 @@ public:
   Double_t CorrectionOutR0(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector);
   Double_t CorrectionOutZ0(Bool_t isGlobal, Bool_t type, Double_t cx, Double_t cy, Double_t cz, Int_t sector);
   
+  Double_t RPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky, Float_t qMax, Float_t threhsold);
+  Double_t SRPhiCOGCorrection(Int_t sector, Int_t padrow, Float_t pad, Float_t cy, Float_t y, Float_t z, Float_t ky,  Float_t qMax, Float_t threhsold);
+  //
+  Double_t GetEdgeQ0(Int_t sector, Int_t padrow, Float_t y);
+  static   Double_t SGetEdgeQ0(Int_t sector, Int_t padrow, Float_t y);
+  //
+  // IROC -OROC+Quadrant alignment
+  //
+  void     AddCorrectionSector(TObjArray & sideAPar, TObjArray &sideCPar, TObjArray & sideACov, TObjArray &sideCCov, Bool_t reset);
+  Double_t GetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant =-1); 
+  static Double_t SGetCorrectionSector(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz, Int_t quadrant=-1); 
+  //
+  // Global alignment
+  //
+  Double_t GetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz);
+  static Double_t SGetCorrection(Int_t coord, Int_t sector, Double_t lx, Double_t ly, Double_t lz);
 public: 
   //
   // Correction out
@@ -43,11 +58,33 @@ public:
   TObjArray   fErrorsOutZ;       // Parameters  for z      distortion  - outer field cage
   Int_t       fParamOutZVersion;  // version of the parameterization
   //
+  //  Edge rfi
+  // 
+  //
+  // Alignment part
+  //
+  Double_t fXIO;               // OROC-IROC boundary
+  Double_t fXmiddle;           // center of the TPC sector local X
+  Double_t fXquadrant;         // x quadrant
+  //
+  // IROC OROC alignment
+  //
+  TObjArray fArraySectorIntParam; // array of sector alignment parameters
+  TObjArray fArraySectorIntCovar; // array of sector alignment covariances 
+  //
+  // Kalman filter for global alignment
+  //
+  TMatrixD  *fSectorParam;     // Kalman parameter   
+  TMatrixD  *fSectorCovar;     // Kalman covariance  
+  //
+  //
+  //
 private:
+
   AliTPCPointCorrection(const AliTPCPointCorrection&); 
   AliTPCPointCorrection& operator=(const AliTPCPointCorrection&); 
   static AliTPCPointCorrection*   fgInstance; //! Instance of this class (singleton implementation)
-  ClassDef(AliTPCPointCorrection, 1); 
+  ClassDef(AliTPCPointCorrection, 3); 
 };
 
 #endif
index b77b3ca..6c1a930 100644 (file)
@@ -63,6 +63,7 @@ AliTPCcalibBase::AliTPCcalibBase():
     fTriggerMaskAccept(-1),   //trigger mask - accept trigger
     fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
     fRejectLaser(kTRUE),                 //flag- reject laser
+    fTriggerClass(),
     fDebugLevel(0)
 {
   //
@@ -83,6 +84,7 @@ AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
   fTriggerMaskAccept(-1),   //trigger mask - accept trigger
   fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
   fRejectLaser(kTRUE),                 //flag- reject laser
+  fTriggerClass(),
   fDebugLevel(0)
 {
   //
@@ -166,7 +168,9 @@ void    AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
   fTime    = event->GetTimeStamp();
   fTrigger = event->GetTriggerMask();
   fMagF    = event->GetMagneticField();
+  fTriggerClass = event->GetFiredTriggerClasses().Data();
   fHasLaser = HasLaser(event); 
+  
 }
 
 
index 78e2d62..44cb2e1 100644 (file)
@@ -9,6 +9,7 @@
 ////
 
 #include "TNamed.h"
+#include "TObjString.h"
 class AliTPCseed;
 class AliESDEvent;
 class AliESDtrack;
@@ -54,6 +55,7 @@ protected:
   Int_t   fTriggerMaskAccept;           //trigger mask - accept
   Bool_t  fHasLaser;                    //flag the laser is overlayed with given event
   Bool_t  fRejectLaser;                 //flag- reject laser
+  TObjString fTriggerClass;                // trigger class
 private:
   Int_t  fDebugLevel;                   //  debug level
 
index c08b82b..e308d66 100644 (file)
@@ -60,6 +60,7 @@
 #include "AliTPCTransform.h"
 #include "AliTPCclusterMI.h"
 #include "AliTPCseed.h"
+#include "AliTPCPointCorrection.h"
 
 ClassImp(AliTPCcalibCalib)
 
@@ -149,10 +150,15 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
   //
 
   //
+  // 0 - Setup transform object
+  //
+  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+  transform->SetCurrentRun(fRun);
+  transform->SetCurrentTimeStamp((UInt_t)fTime);
+  //
   // First apply calibration
   //
-
-  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform();
+  AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
   for (Int_t irow=0;irow<159;irow++) {
     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
     if (!cluster) continue; 
@@ -169,8 +175,38 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
     Float_t dz =AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
     //x[1]-=dy;
     //x[2]-=dz;
+    //
+    // Apply sector alignment
+    //
+    Double_t dxq = AliTPCPointCorrection::SGetCorrectionSector(0,cluster->GetDetector()%36,cluster->GetX(),
+                                                              cluster->GetY(),cluster->GetZ());
+    Double_t dyq = AliTPCPointCorrection::SGetCorrectionSector(1,cluster->GetDetector()%36,cluster->GetX(),
+                                                              cluster->GetY(),cluster->GetZ());
+    Double_t dzq = AliTPCPointCorrection::SGetCorrectionSector(2,cluster->GetDetector()%36,cluster->GetX(),
+                                                              cluster->GetY(),cluster->GetZ());
+    if (kTRUE){
+      x[0]-=dxq;
+      x[1]-=dyq;
+      x[2]-=dzq;
+    }
+    //
+    // Apply r-phi correction  - To be done on track level- knowing the track angle !!!
+    //
+    Double_t corrclY =  
+      corr->RPhiCOGCorrection(cluster->GetDetector(),cluster->GetRow(), cluster->GetPad(),
+                                 cluster->GetY(),cluster->GetY(), cluster->GetZ(), 0., cluster->GetMax(),2.5);
+    // R correction
+    Double_t corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,cluster->GetX(),cluster->GetY(),cluster->GetZ(),cluster->GetDetector());
+
+    if (kTRUE){
+      if (cluster->GetY()>0) x[1]+=corrclY;  // rphi correction
+      if (cluster->GetY()<0) x[1]-=corrclY;  // rphi correction
+      x[0]+=corrR;                           // radial correction
+    }
 
     //
+    //
+    //
     cluster->SetX(x[0]);
     cluster->SetY(x[1]);
     cluster->SetZ(x[2]);
@@ -182,15 +218,23 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
          "event="<<fEvent<<          //  event number
          "time="<<fTime<<            //  time stamp of event
          "trigger="<<fTrigger<<      //  trigger
+         "triggerClass="<<&fTriggerClass<<      //  trigger      
          "mag="<<fMagF<<             //  magnetic field
          "cl0.="<<&cl0<<
          "cl.="<<cluster<<
          "cy="<<dy<<
          "cz="<<dz<<
+         "cR="<<corrR<<
+         "dxq="<<dxq<<
+         "dyq="<<dyq<<
+         "dzq="<<dzq<<
          "\n";
       }
     }
   }
+  //
+  //
+  //
   Int_t ncl = seed->GetNumberOfClusters();
   Double_t covar[15];
   for (Int_t i=0;i<15;i++) covar[i]=0;
@@ -207,16 +251,43 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
   AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();
 
   AliExternalTrackParam trackIn  = *trackOutOld;
-  AliExternalTrackParam trackOut = *trackInOld;
   trackIn.ResetCovariance(200.);
-  trackOut.ResetCovariance(200.);
   trackIn.AddCovariance(covar);
-  trackOut.AddCovariance(covar);
   Double_t xyz[3];
   Int_t nclIn=0,nclOut=0;
   //
+  // Refit in
+  //
+
+  for (Int_t irow=159; irow>0; irow--){
+    AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+    if (!cl) continue;
+    if (cl->GetX()<80) continue;
+    Int_t sector = cl->GetDetector();
+    Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
+    if (TMath::Abs(dalpha)>0.01)
+      trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
+    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
+    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
+    AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    trackIn.GetXYZ(xyz);
+    Double_t bz = AliTracker::GetBz(xyz);
+
+    if (!trackIn.PropagateTo(r[0],bz)) continue;
+    if (RejectCluster(cl,&trackIn)) continue;
+    nclIn++;
+    trackIn.Update(&r[1],cov);    
+  }
+  //
+  AliExternalTrackParam trackOut = trackIn;
+  trackOut.ResetCovariance(200.);
+  trackOut.AddCovariance(covar);
+  //
   // Refit out
   //
+  //Bool_t lastEdge=kFALSE;
   for (Int_t irow=0; irow<160; irow++){
     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
     if (!cl) continue;
@@ -234,14 +305,22 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
     cov[2]*=cov[2];
     trackOut.GetXYZ(xyz);
     Double_t bz = AliTracker::GetBz(xyz);
-    if (trackOut.PropagateTo(r[0],bz)) nclOut++;
+    if (!trackOut.PropagateTo(r[0],bz)) continue;
     if (RejectCluster(cl,&trackOut)) continue;
-    trackOut.Update(&r[1],cov);    
+    nclOut++;
+    trackOut.Update(&r[1],cov);        
+    //if (cl->GetType()<0) lastEdge=kTRUE;
+    //if (cl->GetType()>=0) lastEdge=kFALSE;    
   }
   //
-  // Refit in
   //
-
+  //
+  nclIn=0;
+  trackIn  = trackOut;
+  trackIn.ResetCovariance(10.);
+  //
+  // Refit in one more time
+  //
   for (Int_t irow=159; irow>0; irow--){
     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
     if (!cl) continue;
@@ -258,10 +337,13 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
     trackIn.GetXYZ(xyz);
     Double_t bz = AliTracker::GetBz(xyz);
 
-    if (trackIn.PropagateTo(r[0],bz)) nclIn++;
+    if (!trackIn.PropagateTo(r[0],bz)) continue;
     if (RejectCluster(cl,&trackIn)) continue;
+    nclIn++;
     trackIn.Update(&r[1],cov);    
   }
+
+
   trackIn.Rotate(trackInOld->GetAlpha());
   trackOut.Rotate(trackOutOld->GetAlpha());
   //
@@ -282,6 +364,7 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
        "event="<<fEvent<<          //  event number
        "time="<<fTime<<            //  time stamp of event
        "trigger="<<fTrigger<<      //  trigger
+       "triggerClass="<<&fTriggerClass<<      //  trigger
        "mag="<<fMagF<<             //  magnetic field
        "nclIn="<<nclIn<<
        "nclOut="<<nclOut<<
@@ -302,6 +385,7 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
   AliExternalTrackParam *t = &trackIn;
   track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
   seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
+  seed->SetNumberOfClusters((nclIn+nclOut)*0.5);
   return kTRUE;
 }
 
@@ -312,12 +396,22 @@ Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackPara
   // check the acceptance of cluster
   // Cut on edge effects
   //
+  Float_t kEdgeCut=2.5;
+  Float_t kSigmaCut=6;
+
   Bool_t isReject = kFALSE;
   Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
   Float_t dist  = edgeY - TMath::Abs(cl->GetY());
   if (param)  dist  = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
-  if (dist<3) isReject=kTRUE;
-  if (cl->GetType()<0) isReject=kTRUE;
+  if (dist<kEdgeCut) isReject=kTRUE;
+
+  Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation    
+  AliTPCseed::GetError(cl, param,cov[0],cov[2]);
+  Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
+  Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
+  //
+  if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
+  
   return isReject;
 }
 
index 4169f11..610aa5f 100644 (file)
@@ -56,6 +56,7 @@
 #include "AliESDfriend.h"
 #include "AliESDInputHandler.h"
 #include "AliAnalysisManager.h"
+#include "AliExternalComparison.h"
 
 #include "AliTracker.h"
 #include "AliMagF.h"
 #include "AliLog.h"
 
 #include "AliTPCcalibCosmic.h"
-#include "AliExternalComparison.h"
 #include "TTreeStream.h"
 #include "AliTPCTracklet.h"
+#include "AliESDcosmic.h"
+
 
 ClassImp(AliTPCcalibCosmic)
 
 
 AliTPCcalibCosmic::AliTPCcalibCosmic() 
   :AliTPCcalibBase(),
-   fGainMap(0),
+   fComp(0),
    fHistNTracks(0),
    fClusters(0),
    fModules(0),
@@ -92,7 +94,7 @@ AliTPCcalibCosmic::AliTPCcalibCosmic()
 
 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 
   :AliTPCcalibBase(),
-   fGainMap(0),
+   fComp(0),
    fHistNTracks(0),
    fClusters(0),
    fModules(0),
@@ -124,6 +126,7 @@ AliTPCcalibCosmic::~AliTPCcalibCosmic(){
   //
   //
   //
+  delete fComp;
 }
 
 
@@ -151,7 +154,13 @@ void AliTPCcalibCosmic::Process(AliESDEvent *event) {
   Int_t ntracks=event->GetNumberOfTracks(); 
   fHistNTracks->Fill(ntracks);
   if (ntracks==0) return;
-
+  AliESDcosmic cosmicESD;    
+  TTreeSRedirector * cstream =  GetDebugStreamer();
+  cosmicESD.SetDebugStreamer(cstream);
+  cosmicESD.ProcessEvent(event);
+  if (cstream) cosmicESD.DumpToTree();
+      
+  
 }
 
 
@@ -205,11 +214,11 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 
    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
-     fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
+     fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
      //
-     if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
+     if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
      //
-     if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap)>300) {
+     if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
        TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
        if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
        if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
@@ -246,14 +255,14 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
       if (! seed0) continue;
       if (! seed1) continue;
-      Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
-      Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
+      Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
+      Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
       //
-      Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
-      Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
+      Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
+      Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
       //
-      Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
-      Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
+      Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
+      Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
       //
       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
       Float_t d0  = track0->GetLinearD(0,0);
@@ -298,6 +307,17 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
       Bool_t isPair = IsPair(&param0,&param1);
       //
       if (isPair) FillAcordeHist(track0);
+      if (fComp){
+       Bool_t acceptComp =   fComp->AcceptPair(&param0,&param1);
+       if (acceptComp) fComp->Process(&param0,&param1);
+      }
+      //
+      // combined track params 
+      //
+      AliExternalTrackParam *par0U=MakeCombinedTrack(&param0,&param1);
+      AliExternalTrackParam *par1U=MakeCombinedTrack(&param1,&param0);
+
+
       //
       if (fStreamLevel>0){
        TTreeSRedirector * cstream =  GetDebugStreamer();
@@ -316,6 +336,7 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
            "event="<<fEvent<<          //  event number
            "time="<<fTime<<            //  time stamp of event
            "trigger="<<fTrigger<<      //  trigger
+           "triggerClass="<<&fTriggerClass<<      //  trigger
            "mag="<<fMagF<<             //  magnetic field
            "dir="<<dir<<               //  direction
            "OK="<<isPair<<             //  will be accepted
@@ -332,6 +353,8 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
            "Ip1.="<<ip1<<              //  inner param - lower
            "Op0.="<<op0<<              //  outer param - upper
            "Op1.="<<op1<<              //  outer param - lower
+           "Up0.="<<par0U<<           //  combined track 0
+           "Up1.="<<par1U<<           //  combined track 1
            //
            "v00="<<dvertex0[0]<<       //  distance using kalman
            "v01="<<dvertex0[1]<<       // 
@@ -374,7 +397,9 @@ void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
            "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
            "\n";
        }
-      }      
+      }
+      delete par0U;
+      delete par1U;
     }
   }  
 }    
@@ -432,7 +457,7 @@ Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
     fDeDxMIP->Add(cal->GetHistMIP());
   
   }
-  
+  //if (fComp && cal->fComp) fComp->Add(cal->fComp);
   return 0;
   
 }
@@ -519,19 +544,6 @@ void AliTPCcalibCosmic::BinLogX(TH1 *h) {
   
 }
 
-AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
-{
-  //
-  // Invert paramerameter  - not covariance yet
-  //
-  AliExternalTrackParam *output = new AliExternalTrackParam(*input);
-  Double_t * param = (Double_t*)output->GetParameter();
-  param[0]*=-1;
-  param[3]*=-1;
-  param[4]*=-1;
-  //
-  return output;
-}
 
 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
   //
@@ -539,10 +551,12 @@ AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam
   //
   AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
   par1R->Rotate(track0->GetAlpha());
+  par1R->PropagateTo(track0->GetX(),AliTracker::GetBz()); 
   //
   //
   Double_t * param = (Double_t*)par1R->GetParameter();
   Double_t * covar = (Double_t*)par1R->GetCovariance();
+
   param[0]*=1;  //OK
   param[1]*=1;  //OK
   param[2]*=1;  //?
@@ -552,15 +566,23 @@ AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam
   covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
   //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
   covar[13]*=-1.;
-  par1R->PropagateTo(track0->GetX(),0); // bz shold be set -
-  //if (1){
-  //  printf("Print param\n");
-  //  track1->Print();
-  //  par1R->Print();
-  //}
   return par1R;
 }
 
+AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
+  //
+  // Make combined track
+  //
+  //
+  AliExternalTrackParam * par1T = MakeTrack(track0,track1);
+  AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
+  //
+  UpdateTrack(*par0U,*par1T);
+  delete par1T;
+  return par0U;
+}
+
+
 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
   //
   // Update track 1 with track 2
@@ -588,19 +610,21 @@ void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExte
   //
   // copy data to the matrix
   for (Int_t ipar=0; ipar<5; ipar++){
-    vecXk(ipar,0) = param1[ipar];
-    vecZk(ipar,0) = param2[ipar];
     for (Int_t jpar=0; jpar<5; jpar++){
       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
+      matHk(ipar,jpar)=0;
+      mat1(ipar,jpar)=0;
     }
+    vecXk(ipar,0) = param1[ipar];
+    vecZk(ipar,0) = param2[ipar];
+    matHk(ipar,ipar)=1;
+    mat1(ipar,ipar)=0;
   }
   //
   //
   //
   //
-  matHk(0,0)=1; matHk(1,1)= 1; matHk(2,2)= 1;  
-  matHk(3,3)= 1;    matHk(4,4)= 1;           // vector to measurement
   //
   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
   matHkT=matHk.T(); matHk.T();
@@ -608,7 +632,6 @@ void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExte
   matSk.Invert();
   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
   vecXk += matKk*vecYk;                      //  updated vector 
-  mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
   covXk2 = (mat1-(matKk*matHk));
   covOut =  covXk2*covXk; 
   //
@@ -703,405 +726,3 @@ void AliTPCcalibCosmic::ProcessTree(TTree * chainTracklet, AliExternalComparison
 
 
 
-
-
-
-/*
-
-void Init(){
-  
-.x ~/UliStyle.C
-.x ~/rootlogon.C
-gSystem->Load("libSTAT.so");
-gSystem->Load("libANALYSIS");
-gSystem->Load("libTPCcalib");
-gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
-
-gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
-AliXRDPROOFtoolkit tool; 
-TChain * chainCosmic = tool.MakeChain("cosmic.txt","Track0",0,1000000);
-chainCosmic->Lookup();
-
-TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
-TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
-TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<20");   // OK
-TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
-TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>50");
-TCut cutM("cutM","abs(mag)>0.01");
-TCut cutA=cutT+cutD+cutPt+cutN+cutP1+"trigger!=16";
-
-TCut cuthpt("abs(Tr0.fP[4])+abs(Tr1.fP[4])<0.2");
-TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
-
-//
-chainCosmic->Draw(">>listELP",cutA,"entryList");
-//TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
-//TEntryList *elist = (TEntryList*)gProof->GetOutputList()->At(1);
-chainCosmic->SetEntryList(elist);
-//
-chainCosmic->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
-TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
-chainCosmic->SetEntryList(elistV40Z100);
-
-//
-// Aliases
-//
-chainCosmic->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
-chainCosmic->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
-chainCosmic->SetAlias("signy","(-1+(Tr0.fP[0]>0)*2)");
-
-chainCosmic->SetAlias("dy","Tr0.fP[0]+Tr1.fP[0]");
-chainCosmic->SetAlias("dz","Tr0.fP[1]-Tr1.fP[1]");
-chainCosmic->SetAlias("d1pt","Tr0.fP[4]+Tr1.fP[4]");
-chainCosmic->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
-chainCosmic->SetAlias("dphi","(Tr0.fAlpha-Tr1.fAlpha-pi)");
-
-chainCosmic->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
-chainCosmic->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
-chainCosmic->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
-
-
-
-chainCosmic->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdyA->FitSlicesY();
-hisdyA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdyA_2->SetYTitle("#sigma_{y}(cm)");
-hisdyA_2->SetTitle("Cosmic - Y matching");
-hisdyA_2->SetMaximum(0.5);
-
-
-chainCosmic->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdyC->FitSlicesY();
-hisdyC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdyC_2->SetYTitle("#sigma_{y}(cm)");
-hisdyC_2->SetTitle("Cosmic - Y matching");
-hisdyC_2->SetMaximum(1);
-hisdyC_2->SetMinimum(0);
-hisdyC_2->SetMarkerStyle(22);
-hisdyA_2->SetMarkerStyle(21);
-hisdyC_2->SetMarkerSize(1.5);
-hisdyA_2->SetMarkerSize(1.5);
-hisdyC_2->Draw();
-hisdyA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/ymatching.gif")
-
-chainCosmic->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdzA->FitSlicesY();
-hisdzA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdzA_2->SetYTitle("#sigma_{z}(cm)");
-hisdzA_2->SetTitle("Cosmic - Z matching - A side ");
-hisdzA_2->SetMaximum(0.5);
-
-chainCosmic->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdzC->FitSlicesY();
-hisdzC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdzC_2->SetYTitle("#sigma_{z}(cm)");
-hisdzC_2->SetTitle("Cosmic - Z matching");
-hisdzC_2->SetMaximum(0.5);
-hisdzC_2->SetMarkerStyle(22);
-hisdzA_2->SetMarkerStyle(21);
-hisdzC_2->SetMarkerSize(1.5);
-hisdzA_2->SetMarkerSize(1.5);
-
-hisdzC_2->Draw();
-hisdzA_2->Draw("same");
-
-
-//
-// PICTURE 1/pt
-//
-chainCosmic->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptA(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0"+cutM);
-hisd1ptA->FitSlicesY();
-hisd1ptA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisd1ptA_2->SetYTitle("#sigma_{z}(cm)");
-hisd1ptA_2->SetTitle("Cosmic - Z matching - A side ");
-hisd1ptA_2->SetMaximum(0.5);
-
-chainCosmic->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptC(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0"+cutM);
-hisd1ptC->FitSlicesY();
-hisd1ptC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisd1ptC_2->SetYTitle("#sigma_{1/pt}(1/GeV)");
-hisd1ptC_2->SetTitle("Cosmic - 1/pt matching");
-hisd1ptC_2->SetMaximum(0.05);
-hisd1ptC_2->SetMarkerStyle(22);
-hisd1ptA_2->SetMarkerStyle(21);
-hisd1ptC_2->SetMarkerSize(1.5);
-hisd1ptA_2->SetMarkerSize(1.5);
-
-hisd1ptC_2->Draw();
-hisd1ptA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/1ptmatching.gif")
-
-//
-// Theta
-//
-chainCosmic->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdthetaA->FitSlicesY();
-hisdthetaA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdthetaA_2->SetYTitle("#sigma_{#theta}(cm)");
-hisdthetaA_2->SetTitle("Cosmic - Z matching - A side ");
-hisdthetaA_2->SetMaximum(0.5);
-
-chainCosmic->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdthetaC->FitSlicesY();
-hisdthetaC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdthetaC_2->SetYTitle("#sigma_{#theta}(rad)");
-hisdthetaC_2->SetTitle("Cosmic - Theta matching");
-hisdthetaC_2->SetMaximum(0.01);
-hisdthetaC_2->SetMinimum(0.0);
-hisdthetaC_2->SetMarkerStyle(22);
-hisdthetaA_2->SetMarkerStyle(21);
-hisdthetaC_2->SetMarkerSize(1.5);
-hisdthetaA_2->SetMarkerSize(1.5);
-
-hisdthetaC_2->Draw();
-hisdthetaA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/thetamatching.gif")
-//
-// Phi
-//
-chainCosmic->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
-hisdphiA->FitSlicesY();
-hisdphiA_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdphiA_2->SetYTitle("#sigma_{#phi}(rad)");
-hisdphiA_2->SetTitle("Cosmic - Z matching - A side ");
-hisdphiA_2->SetMaximum(0.5);
-
-chainCosmic->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
-hisdphiC->FitSlicesY();
-hisdphiC_2->SetXTitle("#sqrt{1/p_{t}}");
-hisdphiC_2->SetYTitle("#sigma_{#phi}(rad)");
-hisdphiC_2->SetTitle("Cosmic - Phi matching");
-hisdphiC_2->SetMaximum(0.01);
-hisdphiC_2->SetMinimum(0.0);
-hisdphiC_2->SetMarkerStyle(22);
-hisdphiA_2->SetMarkerStyle(21);
-hisdphiC_2->SetMarkerSize(1.5);
-hisdphiA_2->SetMarkerSize(1.5);
-
-hisdphiC_2->Draw();
-hisdphiA_2->Draw("same");
-gPad->SaveAs("~/Calibration/Cosmic/pic/phimatching.gif")
-
-
-
-}
-
-
-*/
-
-
-/*
-void MatchTheta(){
-
-TStatToolkit toolkit;
-Double_t chi2=0;
-Int_t    npoints=0;
-TVectorD fitParamA0;
-TVectorD fitParamA1;
-TVectorD fitParamC0;
-TVectorD fitParamC1;
-TMatrixD covMatrix;
-
-
-TString fstring="";
-// 
-fstring+="mtheta++";
-fstring+="ca++";
-fstring+="sa++";
-fstring+="ca*mtheta++";
-fstring+="sa*mtheta++";
-//
-fstring+="side++";
-fstring+="side*mtheta++";
-fstring+="side*ca++";
-fstring+="side*sa++";
-fstring+="side*ca*mtheta++";
-fstring+="side*sa*mtheta++";
-
-
-TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt&&!crossI&&!crossO", chi2,npoints,fitParamA0,covMatrix,0.8);
-chainCosmic->SetAlias("dtheta0",strTheta0.Data());
-strTheta0->Tokenize("+")->Print();
-
-
-//fstring+="mtheta++";
-//fstring+="mtheta^2++";
-//fstring+="ca*mtheta^2++";
-//fstring+="sa*mtheta^2++";
-
-
-
-}
-
-*/
-
-
-
-
-/*
- void PosCorrection()
-
-
- TStatToolkit toolkit;
- Double_t chi2=0;
- Int_t    npoints=0;
- TVectorD fitParam;
- TMatrixD covMatrix;
- //Theta
-chainCosmic->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
-chainCosmic->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
-chainCosmic->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
-//
-//
-TString strFit="";
-//
-strFit+="sign++";                              //1
-strFit+="Tr0.fP[3]++";                         //2
-// 
-strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //3
-strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //4
-//
-strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //5
-strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //6  
-//
-strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++";         //7
-strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++";         //8
- //                                        
- TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
- chainCosmic->SetAlias("corTheta",thetaParam->Data());
- chainCosmic->Draw("dthe:Tr0.fP[1]","","",50000);
-
-
-
-*/
-
-
-
-/*
-
-void AliTPCcalibCosmic::dEdxCorrection(){
-  TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
-  TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
-  TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
-  TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
-  TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
-  TCut cutA=cutT+cutD+cutPt+cutN+cutS;
-
-
- .x ~/UliStyle.C
-  gSystem->Load("libANALYSIS");
-  gSystem->Load("libTPCcalib");
-  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
-  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
- AliXRDPROOFtoolkit tool; 
-  TChain * chainCosmic = tool.MakeChain("cosmic.txt","Track0",0,1000000);
-  chainCosmic->Lookup();
-  
-  chainCosmic->Draw(">>listEL",cutA,"entryList");
-  TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
-  chainCosmic->SetEntryList(elist);
-
-  .x ~/rootlogon.C
-   gSystem->Load("libSTAT.so");
-   TStatToolkit toolkit;
-  Double_t chi2=0;
-  Int_t    npoints=0;
-  TVectorD fitParam;
-  TMatrixD covMatrix;
-  
-  chainCosmic->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
-  
-  TString strFit;
-  strFit+="(Tr0.fP[1]/250)++";
-  strFit+="(Tr0.fP[1]/250)^2++";
-  strFit+="(Tr0.fP[3])++";
-  strFit+="(Tr0.fP[3])^2++";
-
-  TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix) 
-
-
-
-*/
-
-
-/*
-gSystem->Load("libANALYSIS");
-gSystem->Load("libSTAT");
-gSystem->Load("libTPCcalib");
-
-TStatToolkit toolkit;
-Double_t chi2;
-TVectorD fitParam;
-TMatrixD covMatrix;
-Int_t npoints;
-//
-TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03");  // OK
-TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5");     // OK
-TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
-TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
-TCut cutA=cutT+cutD+cutPt+cutN;
-
-
-
-TTree * chainCosmic = Track0;
-
-
-chainCosmic->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
-//
-chainCosmic->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
-chainCosmic->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
-chainCosmic->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
-chainCosmic->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
-
-TString fstring="";  
-fstring+="dr1++";
-fstring+="dr2++";
-fstring+="dr4++";
-fstring+="dr5++";
-//
-fstring+="dr1*dr2++";
-fstring+="dr1*dr4++";
-fstring+="dr1*dr5++";
-fstring+="dr2*dr4++";
-fstring+="dr2*dr5++";
-fstring+="dr4*dr5++";
-
-
-
-TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
-  
-chainCosmic->SetAlias("corQT",strqdedx->Data());
-
-*/
-
-
-/*
-  chainCosmic->SetProof(kTRUE);
-  chainCosmic->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE):Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)",""+cutA,"",100000);
-
-
-chainCosmic->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)/Seed1.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)>>his(100,0.5,1.5)","min(Orig0.fTPCncls,Orig1.fTPCncls)>130"+cutA,"",50000);
-
-*/
-
-
-/*
-chainCosmic->Draw("Tr0.fP[1]-Tr1.fP[1]:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0&&abs(mag)>0.1"+cutA); 
-
-TGraph *grdzA = (TGraph*)gProof->GetOutputList()->At(1)->Clone();
-
-
-
-*/
-
-
-
-
index ffca65a..35c3e97 100644 (file)
@@ -30,12 +30,15 @@ public:
   //
   void              FindPairs(AliESDEvent *event);
   Bool_t            IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1);
-  void              SetGainMap(AliTPCCalPad *GainMap){fGainMap = GainMap;};
   static void       CalculateBetheParams(TH2F *hist, Double_t * initialParam);
   static Double_t   CalculateMIPvalue(TH1F * hist);
-  AliExternalTrackParam *Invert(AliExternalTrackParam *input);
   AliExternalTrackParam *MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1);
+  AliExternalTrackParam *MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1);
+
   void UpdateTrack(AliExternalTrackParam &track0, const AliExternalTrackParam &track1);
+  void SetComparison(AliExternalComparison * comp) { fComp=comp;}
+  AliExternalComparison *GetComparison() { return fComp;}
+  
   //
   TH1F   *          GetHistNTracks(){return fHistNTracks;};
   TH1F   *          GetHistClusters(){return fClusters;};
@@ -55,7 +58,7 @@ private:
 
   void              FillAcordeHist(AliESDtrack *upperTrack);
 
-  AliTPCCalPad *fGainMap;         //  gain map from Krypton calibration
+  AliExternalComparison * fComp;  //  comparison histogram
   TH1F  *fHistNTracks;            //  histogram showing number of ESD tracks per event
   TH1F  *fClusters;               //  histogram showing the number of clusters per track
   TH2F  *fModules;                //  2d histogram of tracks which are propagated to the ACORDE scintillator array
index 14c86d4..ee6f472 100644 (file)
@@ -350,22 +350,22 @@ AliTPCcalibTracks::AliTPCcalibTracks(const Text_t *name, const Text_t *title, Al
       
       // amplitude
       sprintf(chname,"Amp_Sector%d",i);
-      his1 = new TH1F(chname,chname,250,0,500);         // valgrind 
+      his1 = new TH1F(chname,chname,100,0,32);         // valgrind 
       his1->SetXTitle("Max Amplitude (ADC)");
       fArrayAmp->AddAt(his1,i);
       sprintf(chname,"Amp_Sector%d",i+36);
-      his1 = new TH1F(chname,chname,200,0,600);         // valgrind 3   13,408,208 bytes in 229 blocks are still reachable
+      his1 = new TH1F(chname,chname,100,0,32);         // valgrind 3   13,408,208 bytes in 229 blocks are still reachable
       his1->SetXTitle("Max Amplitude (ADC)");
       fArrayAmp->AddAt(his1,i+36);
       
       // driftlength
       sprintf(chname, "driftlengt vs. charge, ROC %i", i);
-      prof1 = new TProfile(chname, chname, 500, 0, 250);
+      prof1 = new TProfile(chname, chname, 25, 0, 250);
       prof1->SetYTitle("Charge");
       prof1->SetXTitle("Driftlength");
       fArrayChargeVsDriftlength->AddAt(prof1,i);
       sprintf(chname, "driftlengt vs. charge, ROC %i", i+36);
-      prof1 = new TProfile(chname, chname, 500, 0, 250);
+      prof1 = new TProfile(chname, chname, 25, 0, 250);
       prof1->SetYTitle("Charge");
       prof1->SetXTitle("Driftlength");
       fArrayChargeVsDriftlength->AddAt(prof1,i+36);
@@ -889,16 +889,17 @@ void  AliTPCcalibTracks::FillResolutionHistoLocal(AliTPCseed * track){
       
       TProfile *profAmpRow =  (TProfile*)fArrayAmpRow->At(sector);
       if ( irow >= kFirstLargePad) max /= kLargePadSize;
-      profAmpRow->Fill( (Double_t)cluster0->GetRow(), max );
+      Double_t smax = TMath::Sqrt(max);
+      profAmpRow->Fill( (Double_t)cluster0->GetRow(), smax );
       TH1F *hisAmp =  (TH1F*)fArrayAmp->At(sector);
-      hisAmp->Fill(max);
+      hisAmp->Fill(smax);
       
       // remove the following two lines one day:
       TProfile *profDriftLength = (TProfile*)fArrayChargeVsDriftlength->At(sector);
-      profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
+      profDriftLength->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
       
       TProfile *profDriftLengthTmp = (TProfile*)(fcalPadRegionChargeVsDriftlength->GetObject(isegment, padSize));
-      profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), max );
+      profDriftLengthTmp->Fill( 250.-(Double_t)TMath::Abs(cluster0->GetZ()), smax );
       
       fHclusterPerPadrow->Fill(irow);   // fill histogram showing clusters per padrow
       Int_t ipad = TMath::Nint(cluster0->GetPad());