more secure string operations
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibLaser.cxx
index acda597..e44f72c 100644 (file)
@@ -39,7 +39,7 @@
   gSystem->Load("libTPCcalib");
   TFile fcalib("CalibObjectsTrain2.root");
   AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)fcalib->Get("laserTPC");
-  laser->DumpMeanInfo(-0,0)
+  laser->DumpMeanInfo(run)
   TFile fmean("laserMean.root")
   //
   //  laser track clasification;
   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
   AliXRDPROOFtoolkit tool;
-  AliXRDPROOFtoolkit::FilterList("laser.txt","* driftvN",1) 
+  AliXRDPROOFtoolkit::FilterList("laserDebug.list","* driftvN",1) 
   TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
   chainDrift->Lookup();
-  TChain * chainDriftN = tool.MakeChainRandom("laser.txt.Good","driftvN",0,300);
+  TChain * chainDriftN = tool.MakeChainRandom("laserDebug.list.Good","driftvN",0,300);
   chainDriftN->Lookup();
 
 
 #include "TH2F.h"
 #include "TStatToolkit.h"
 #include "TROOT.h"
+#include "TDatabasePDG.h"
 
 
 #include "TTreeStream.h"
 #include "TTimeStamp.h"
 #include "AliDCSSensorArray.h"
 #include "AliDCSSensor.h"
+#include "AliGRPObject.h"
 
 using namespace std;
 
@@ -230,6 +232,12 @@ AliTPCcalibLaser::AliTPCcalibLaser():
   // Constructor
   //
   fTracksEsdParam.SetOwner(kTRUE);
+  for (Int_t i=0; i<336; i++) {
+    fFitZ[i]=0;
+    fCounter[i]=0;    //! counter of usage
+    fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
+    fClusterSatur[i]=0;   //!couter of saturated clusters in "sensitive are"
+  }
 }
 
 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
@@ -337,6 +345,12 @@ AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool
   // Constructor
   //
   fTracksEsdParam.SetOwner(kTRUE);
+  for (Int_t i=0; i<336; i++) {
+    fFitZ[i]=0;
+    fCounter[i]=0;    //! counter of usage
+    fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
+    fClusterSatur[i]=0;   //!couter of saturated clusters in "sensitive are"
+  }
 }
 
 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
@@ -440,6 +454,12 @@ AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
   //
   // copy constructor
   //
+  for (Int_t i=0; i<336; i++) {
+    fFitZ[i]=0;
+    fCounter[i]=0;    //! counter of usage
+    fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
+    fClusterSatur[i]=0;   //!couter of saturated clusters in "sensitive are"
+  }
 }
 
 
@@ -548,7 +568,8 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
   //
   // Loop over tracks and call  Process function
   //
-  Int_t kMinTracks=20;
+  const Int_t  kMinTracks=20;
+  const Int_t  kMinClusters=40;
 
   fESD = event;
   if (!fESD) {
@@ -591,6 +612,7 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
     if (!friendTrack) continue;
     AliESDtrack *track=fESD->GetTrack(i);
+    if (!track) continue;
     Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
     if (binC>336) continue; //remove CE background
     TObject *calibObject=0;
@@ -598,7 +620,7 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
        break;
-    if (track&&seed) {
+    if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
       //filter CE tracks
       Int_t id = FindMirror(track,seed);
       if (id>=0) counter++;
@@ -608,8 +630,8 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) {
   fNtracks=counter;
   if (counter<kMinTracks) return;
 
-  FitDriftV();
-  FitDriftV(0.3);
+  //FitDriftV();
+  FitDriftV(0.2);
   if (!fFullCalib) return;
   static Bool_t init=kFALSE;
   if (!init){
@@ -926,7 +948,7 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
 
   // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
   // 1. Robust fit is used in the itteration number 0
-  // only fraction of laser uted
+  // only fraction of laser used
   // 2. Only the tracks close to the fit used in the second itteration
   /*
     Formulas:
@@ -959,8 +981,8 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
   const Float_t kSaturCut    = 0.05;    // remove saturated lasers - cut on fraction of saturated 
   const Float_t kDistCut     = 3.;      // distance sigma cut - 3 sigma
   const Float_t kDistCutAbs  = 1.;      // absolute cut 1 cm
-  const Float_t kMinClusters = 60.;      // minimal amount of the clusters
-  const Float_t kMinSignal   = 10.;      // minimal mean height of the signal
+  const Float_t kMinClusters = 40.;      // minimal amount of the clusters
+  const Float_t kMinSignal   = 2.5;      // minimal mean height of the signal
   const Float_t kChi2Cut     = 1.0;     // chi2 cut to accept drift fit
   //
   static TLinearFitter fdriftA(3,"hyp2");
@@ -1075,7 +1097,7 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       ltrp->GetXYZ(lxyz);
       ltrp->GetPxPyPz(lpxyz);
       Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
-      if (npointsAC>0) sz =TMath::Sqrt(chi2AC); 
+      //if (npointsAC>0) sz =TMath::Sqrt(chi2AC); 
       if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
       if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
 
@@ -1104,6 +1126,7 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
     if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
       //
       fdriftA.Eval();
+      //if (iter==0) fdriftA.FixParameter(2,0); //fix global y gradient
       npointsA= fdriftA.GetNpoints();
       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
       fdriftA.EvalRobust(kFraction[iter]);
@@ -1119,6 +1142,7 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
     }
     if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
       fdriftC.Eval();
+      //if (iter==0) fdriftC.FixParameter(2,0); //fix global y gradient
       npointsC= fdriftC.GetNpoints();
       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
       fdriftC.EvalRobust(kFraction[iter]);
@@ -1133,8 +1157,9 @@ Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
       }
     }
 
-    if (fdriftAC.GetNpoints()>minFraction*knLaser){
+    if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
       fdriftAC.Eval();
+      //if (iter==0) fdriftAC.FixParameter(2,0); //fix global y gradient
       npointsAC= fdriftAC.GetNpoints();
       chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
       fdriftAC.EvalRobust(kFraction[iter]);
@@ -1352,8 +1377,8 @@ Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
 
 
   AliExternalTrackParam param(*(track->GetOuterParam()));
-  AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
-  AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
+  AliTracker::PropagateTrackTo(&param,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
+  AliTracker::PropagateTrackTo(&param,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
   AliTPCLaserTrack ltr;
   AliTPCLaserTrack *ltrp=0x0;
   //  AliTPCLaserTrack *ltrpjw=0x0;
@@ -2111,7 +2136,7 @@ void AliTPCcalibLaser::RefitLaserJW(Int_t id){
        if (ltrp->GetSide()==0){
          if ((*fFitAside)[1]>0. || fUseFixedDriftV) { 
            // ignore global y dependence for now
-           Double_t zcorrected = 0;        
+           zcorrected = 0;         
            if(!fUseFixedDriftV) 
              zcorrected = (zcl + (*fFitAside)[0] - 
                            (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
@@ -2231,7 +2256,7 @@ void AliTPCcalibLaser::RefitLaserJW(Int_t id){
 
 
 
-void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
+void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
   //
   //  Dump information about laser beams
   //  isOK variable indicates usability of the beam  
@@ -2246,17 +2271,33 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
   const Float_t krmsCut1=0.16;
   const Float_t kmultiCut=2;
   const Float_t kcutP0=0.002;
+  AliMagF* magF=  dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
+  Double_t xyz[3]={90,0,10};         // tmp. global position 
+  Double_t bxyz[3]={90,0,10};        // tmp. mag field  integral - cylindrical
+  Double_t bgxyz[3]={90,0,10};       // tmp. mag field  integral - local
   //
   AliTPCcalibLaser *laser = this;
   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
   TF1 fg("fg","gaus");
-
+  AliTPCParam  * tpcparam    = 0;   
   // start set up for absolute residuals analysis
-  //AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
-  //  AliTPCParam  * tpcparam    = calib->GetParameters(); 
-  AliTPCParam  * tpcparam    = 0; 
+  //
+  AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
+  tpcparam    = calib->GetParameters(); 
   if (!tpcparam) tpcparam    = new AliTPCParamSR;
   tpcparam->Update();
+  AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
+  Float_t current=0;
+  Float_t bfield      = 0, bz=0;
+
+  if (grp){
+    Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
+    current = grp->GetL3Current((AliGRPObject::Stats)0);
+    bfield = polarity*5*current/30000.;
+    bz = polarity*5*current/30000.;
+    printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
+  }
+
   SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
   SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
   TLinearFitter lfabsyInner(2);
@@ -2289,6 +2330,7 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
      AliTPCLaserTrack::LoadTracks();
       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
     }
+    ltrp->UpdatePoints();
     pcstream->GetFile()->cd();
     if (hisphi)  hisphi->Write();
     if (hisphiP) hisphiP->Write();
@@ -2321,11 +2363,15 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
       hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
     Double_t gp31 = fg.GetParameter(1);
     Double_t gp32 = fg.GetParameter(2);
+    Double_t meanp3 = hisP3->GetMean();
+    Double_t rmsp3  = hisP3->GetRMS();
     //
     if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
       hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
     Double_t gp41 = fg.GetParameter(1);
     Double_t gp42 = fg.GetParameter(2);
+    Double_t meanp4 = hisP4->GetMean();
+    Double_t rmsp4  = hisP4->GetRMS();
     //
     Float_t meanS=hisS->GetMean();
     //
@@ -2535,6 +2581,49 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
     TVectorD vecEz(159);       //error z
     TVectorD vecPhi(159);      // local tangent
     TVectorD vecPhiR(159);     // local tangent
+    // magnetic field integrals
+    TVectorD vecIBR(159);        // radial
+    TVectorD vecIBRPhi(159);     // r-phi
+    TVectorD vecIBLX(159);       // local x
+    TVectorD vecIBLY(159);       // local y
+    TVectorD vecIBGX(159);       // local x
+    TVectorD vecIBGY(159);       // local y
+    TVectorD vecIBZ(159);        // z
+    //
+    for (Int_t irow=0;irow<159;irow++){
+      vecIBR[irow]=0;
+      vecIBRPhi[irow]=0;
+      vecIBLX[irow]=0;
+      vecIBLY[irow]=0;
+      vecIBGX[irow]=0;
+      vecIBGY[irow]=0;
+      vecIBZ[irow]=0;
+      Double_t gx    =(*(ltrp->fVecGX))[irow];
+      Double_t gy    =(*(ltrp->fVecGY))[irow];
+      Int_t    lsec  =TMath::Nint((*(ltrp->fVecSec))[irow]);
+      Double_t   ca  =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
+      Double_t   sa  =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
+      xyz[2]=(*(ltrp->fVecGZ))[irow];
+      xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
+      xyz[1]=TMath::ATan2(gy,gx);
+      Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
+      if (magF){
+       magF->GetTPCIntCyl(xyz,bxyz);
+       magF->GetTPCInt(gxyz,bgxyz);
+       vecIBR[irow]=bxyz[0];
+       vecIBRPhi[irow]=bxyz[1];
+       //
+       vecIBGX[irow]=bgxyz[0];
+       vecIBGY[irow]=bgxyz[1];
+       //
+       vecIBLX[irow]=  bgxyz[0]*ca+bgxyz[1]*sa;
+       vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
+       //
+
+       vecIBZ[irow]=bxyz[2];
+      }
+    }
+
 
     lfabsyInner.ClearPoints();    
     lfabszInner.ClearPoints();    
@@ -2561,7 +2650,7 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
     Float_t  lasTanPhiLocIn = -100;
     Float_t  lasTanPhiLocOut = -100;
 
-    if(histAbsY->GetEntries()>0) {
+    if(histAbsY && histAbsY->GetEntries()>0) {
       
       Double_t rotAngOut = 10;
       Double_t rotAngIn = 10;
@@ -2674,7 +2763,8 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
                                  TMath::Max(zprof->GetBinError(bin), 0.001));
          }
        }
-
+       // global position
+       
       }
        
       delete yprof; delete zprof;
@@ -2864,13 +2954,27 @@ void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
       //
       "gp31="<<gp31<<            //gaus mean - tgl
       "gp32="<<gp32<<            //gaus rms  -tgl
+      "meanp3="<<meanp3<<            //mean - tgl
+      "rmsp3="<<rmsp3<<            //rms  -tgl
       "gp41="<<gp41<<            //gaus mean - P4
       "gp42="<<gp42<<            //gaus rms  - P4
+      "meanp4="<<meanp4<<            //mean - P4
+      "rmsp4="<<rmsp4<<            //rms  - P4
       // Parameters from abs res analysis
       "SecIn="<<secInner<<              // inner sector
       "SecOut="<<secOuter<<             // outer sector
       "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
       "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
+      "ibr.="<<&vecIBR<<   // radial filed integral
+      "ibrphi.="<<&vecIBRPhi<<   // r=phifiled integral
+      "ibr.="<<&vecIBR<<   // radial filed integral
+      "ibz.="<<&vecIBZ<<   // z filed integral
+      //
+      "iblx.="<<&vecIBLX<<   // local bx  integral
+      "ibly.="<<&vecIBLY<<   // local by integral
+      "ibgx.="<<&vecIBGX<<   // global bx  integral
+      "ibgy.="<<&vecIBGY<<   // global by integral
+      //
       "X.="<<&vecX<<       // local x 
       "Y.="<<&vecY<<       // local y 
       "R.="<<&vecR<<       // radius 
@@ -2947,7 +3051,7 @@ void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
   //
   for (Int_t id=0; id<336; id++){
     // id =205;
-    sprintf(cut,"fId==%d&&%s",id,cutUser);
+    snprintf(cut,1000,"fId==%d&&%s",id,cutUser);
     Int_t entries = chain->Draw("bz",cut,"goff");
     if (entries<3) continue;
     AliTPCLaserTrack *ltrp = 0;
@@ -3358,7 +3462,8 @@ void   AliTPCcalibLaser::MakeFitHistos(){
   for (Int_t id=0; id<336;id++){
     TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
     TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
-    TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
+    //TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
+    TH2F *profy2 = 0;
     TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
     TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
     TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
@@ -3403,7 +3508,8 @@ void   AliTPCcalibLaser::MakeFitHistos(){
   //
   for (Int_t id=0; id<336;id++){
     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
-    TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
+    //TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
+    TH1F * hisP3 = 0;
     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
     
     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);