Warning removal +
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 22 Mar 2009 14:31:37 +0000 (14:31 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 22 Mar 2009 14:31:37 +0000 (14:31 +0000)
AliTPCcalibPID.h AliTPCcalibPID.cxx -
More ratio histograms - Ration dEdxfrac/dEdxfull

AliTPCcalibUnlinearity.h AliTPCcalibUnlinearity.cxx
Adding Quadrant alignment
Edge effect histograms

TPC/AliTPCcalibPID.cxx
TPC/AliTPCcalibPID.h
TPC/AliTPCcalibUnlinearity.cxx
TPC/AliTPCcalibUnlinearity.h

index 5d6f511..ccdb500 100644 (file)
@@ -25,13 +25,13 @@ The corrections which are done on the cluster level (angle,drift time etc.) are
 TFile f("CalibObjects.root")
 AliTPCcalibPID *cal = f->Get("TPCCalib")->FindObject("calibPID")
 
-cal->GetHistLongMediumRatio()->Projection(0)->Fit("gaus")
-cal->GetHistShortMediumRatio()->Projection(0)->Fit("gaus")
+cal->GetHistRatioMaxTot()->Projection(0)->Fit("gaus")
+cal->GetHistRatioTracklet()->Projection(0)->Fit("gaus")
 
 1.b) Update OCDB:
 .x $ALICE_ROOT/TPC/macros/ConfigOCDB.C
 AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
-(*parcl->fQpadMnorm)[ipad] = oldvalue*corrFactor
+*parcl->fQpadMnorm)[ipad] = oldvalue*corrFactor
 
 Int_t runNumber = 0;
 AliCDBMetaData *metaData= new AliCDBMetaData();
@@ -120,9 +120,11 @@ AliTPCcalibPID::AliTPCcalibPID()
    fLandau(0),
    fDeDxQmax(0),
    fDeDxQtot(0),
-   fDeDxRatio(0),
-   fDeDxShortMediumRatio(0),
-   fDeDxLongMediumRatio(0)
+   fDeDxRatioMaxTot(0),
+   fDeDxRatioQmax(0),
+   fDeDxRatioQtot(0),
+   fDeDxRatioTruncQtot(0),
+   fDeDxRatioTruncQmax(0)
 {  
   //
   // Empty default cosntructor
@@ -146,9 +148,11 @@ AliTPCcalibPID::AliTPCcalibPID(const Text_t *name, const Text_t *title)
    fLandau(0),
    fDeDxQmax(0),
    fDeDxQtot(0),
-   fDeDxRatio(0),
-   fDeDxShortMediumRatio(0),
-   fDeDxLongMediumRatio(0)
+   fDeDxRatioMaxTot(0),
+   fDeDxRatioQmax(0), 
+   fDeDxRatioQtot(0) ,
+   fDeDxRatioTruncQtot(0),
+   fDeDxRatioTruncQmax(0)
 {
   //
   //
@@ -161,30 +165,44 @@ AliTPCcalibPID::AliTPCcalibPID(const Text_t *name, const Text_t *title)
   fUpperTrunc = 0.6;
   //
   fUseShapeNorm = kTRUE;
-  fUsePosNorm = kFALSE;
+  fUsePosNorm = 0;
   fUsePadNorm = 1;
   //
   fIsCosmic  = kTRUE;
   //
-  //                  dE/dx,  z, phi, theta,    p,  bg, ncls
-  Int_t binsQA[7]    = {150, 10,  10,    10,   50,   1,  8};
-  Double_t xminQA[7] = {0.,  0,    0,     0, 0.01, 0.1, 60};
-  Double_t xmaxQA[7] = {10., 1,  0.6,   1.5,   50, 500, 160};
+  //                  dE/dx,  z, phi, theta,    p,  bg, ncls, tracklet type
+  Int_t binsQA[8]    = {150, 10,  10,    10,   50, 50,  8,    5};
+  Double_t xminQA[8] = {0.,  0,    0,     0, 0.01, 0.1, 60,   0};
+  Double_t xmaxQA[8] = {10., 1,  0.6,   1.5,   50, 500, 160,  5};
   //
-  Double_t xminRa[7] = {0.5, 0,    0,     0, 0.01, 0.1, 60};
-  Double_t xmaxRa[7] = { 2., 1,  0.6,   1.5,   50, 500, 160};
+  //
+  //
+  //                  dE/dx,  z, phi, theta, dEdx, dEdx*dl, ncls, tracklet type
+  Int_t    binsRA[9] = {100, 10,  10,    10,   25,  25,      4,    5};
+  Double_t xminRa[9] = {0.5, 0,    0,     0,  0.2, 0.2,     60,    0};
+  Double_t xmaxRa[9] = {4.0, 1,  0.6,   1.5,    5,   5,    160,    5};
   
   // z,sin(phi),tan(theta),p,betaGamma,ncls
-  fDeDxQmax  = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminQA,xmaxQA);
-  fDeDxQtot  = new THnSparseS("fDeDxQtot","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminQA,xmaxQA);
-  fDeDxRatio = new THnSparseS("fDeDxRatio","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
-  fDeDxShortMediumRatio = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
-  fDeDxLongMediumRatio  = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls); TPC signal Qmax (a.u.)",7,binsQA,xminRa,xmaxRa);
+  fDeDxQmax  = new THnSparseS("fDeDxQmax","Qmax;(z,sin(phi),tan(theta),p,betaGamma,ncls,type); TPC signal Qmax (a.u.)",8,binsQA,xminQA,xmaxQA);
+  fDeDxQtot  = new THnSparseS("fDeDxQtot","Q_{tot};(z,sin(phi),tan(theta),p,betaGamma,ncls,type); TPC signal Qmax (a.u.)",8,binsQA,xminQA,xmaxQA);
+  //
+  // ratio histograms
+  //
+  fDeDxRatioMaxTot = new THnSparseS("fDeDxRatioMaxTot","Q_{max}/Q_{tot};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type); TPC signal Qmax/Qtot (a.u.)",8,binsRA,xminRa,xmaxRa);
+  fDeDxRatioQmax = new THnSparseS("fDeDxRatioQmax","Q_{mtracklet}/Q_{mtrack};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Tracklet/Track (a.u.)",8,binsRA,xminRa,xmaxRa);
+  fDeDxRatioQtot = new THnSparseS("fDeDxRatioQtot","Q_{ttracklet}/Q_{ttrack};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Tracklet/Track (a.u.)",8,binsRA,xminRa,xmaxRa);
+  fDeDxRatioTruncQmax = new THnSparseS("fDeDxRatioTrunQmax","Q_{max}/Q_{maxtrunc};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Full/Trunc. (a.u.)",8,binsRA,xminRa,xmaxRa);
+  fDeDxRatioTruncQtot = new THnSparseS("fDeDxRatioTruncQtot","Q_{tot}/Q_{tottrunc};(z,sin(phi),tan(theta),dEdx,dEdx*dl,ncls,type,qtupe); TPC signal Full/Trunc (a.u.)",8,binsRA,xminRa,xmaxRa);
+
+
   BinLogX(fDeDxQmax,4); BinLogX(fDeDxQmax,5);
   BinLogX(fDeDxQtot,4); BinLogX(fDeDxQtot,5);
-  BinLogX(fDeDxRatio,4); BinLogX(fDeDxRatio,5);
-  BinLogX(fDeDxShortMediumRatio,4); BinLogX(fDeDxShortMediumRatio,5);
-  BinLogX(fDeDxLongMediumRatio,4); BinLogX(fDeDxLongMediumRatio,5);
+  //
+  BinLogX(fDeDxRatioMaxTot,4); BinLogX(fDeDxRatioMaxTot,5);
+  BinLogX(fDeDxRatioQmax,4); BinLogX(fDeDxRatioQmax,5);
+  BinLogX(fDeDxRatioQtot,4); BinLogX(fDeDxRatioQtot,5);
+  BinLogX(fDeDxRatioTruncQmax,4); BinLogX(fDeDxRatioTruncQmax,5);
+  BinLogX(fDeDxRatioTruncQtot,4); BinLogX(fDeDxRatioTruncQtot,5);
   //
   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",40,0,160);
@@ -200,6 +218,7 @@ AliTPCcalibPID::~AliTPCcalibPID(){
   //
   //
   //
+  
 }
 
 
@@ -212,7 +231,7 @@ void AliTPCcalibPID::Process(AliESDEvent *event) {
     Printf("ERROR: ESD not available");
     return;
   }  
-
+  const Int_t kMinClustersTracklet = 25;
   Int_t ntracks=event->GetNumberOfTracks(); 
   fHistNTracks->Fill(ntracks);
   
@@ -238,7 +257,7 @@ void AliTPCcalibPID::Process(AliESDEvent *event) {
     // calculate necessary track parameters
     //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
     Double_t meanP = trackIn->GetP();
-    Double_t d = trackIn->GetLinearD(0,0);
+    //Double_t d = trackIn->GetLinearD(0,0);
     Int_t NclsDeDx = track->GetTPCNcls();
 
     //if (meanP > 0.7 || meanP < 0.2) continue;
@@ -268,16 +287,16 @@ void AliTPCcalibPID::Process(AliESDEvent *event) {
       if (meanP > 30 && TMath::Abs(trackIn->GetSnp()) < 0.2) fClusters->Fill(track->GetTPCNcls());
       // calculate dEdx
       // (Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Bool_t posNorm, Int_t padNorm, Int_t returnVal)
-      Double_t TPCsignalTot    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+      //Double_t TPCsignalTot    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
       Double_t TPCsignalMax    = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
       //
-      Double_t TPCsignalShort  = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,63,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
-      Double_t TPCsignalMedium = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63,63+64,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
-      Double_t TPCsignalLong   = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63+64,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+      //Double_t TPCsignalShort  = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,63,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+      //Double_t TPCsignalMedium = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63,63+64,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+      //Double_t TPCsignalLong   = seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,63+64,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
       //
       //Double_t driftMismatch = seed->GetDriftTimeMismatch(0,0.95,0.5);
       Double_t driftMismatch = 0;
-      Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;        
+      //      Double_t drift = 1 - (TMath::Abs(trackIn->GetZ()) + TMath::Abs(trackOut->GetZ()))/500.;        
       
       // particle identification
       Double_t mass = 0.105658369;// muon
@@ -288,26 +307,113 @@ void AliTPCcalibPID::Process(AliESDEvent *event) {
        fLandau->Fill(0.1,0.6);
       }
       //var. of interest, z,sin(phi),tan(theta),p,betaGamma,ncls
-      Double_t snp = TMath::Abs(trackIn->GetSnp());
-      Double_t tgl = TMath::Abs(trackIn->GetTgl());
+      Double_t snpIn = TMath::Abs(trackIn->GetSnp());
+      Double_t snpOut = TMath::Abs(trackOut->GetSnp());
+      Double_t tglIn = TMath::Abs(trackIn->GetTgl());
+      Double_t tglOut = TMath::Abs(trackOut->GetTgl());
+      Double_t driftIn = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
+      Double_t driftOut = 1 - (TMath::Abs(trackIn->GetZ()))/250.;
       //
-      Double_t vecQmax[7] = {TPCsignalMax,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
-      Double_t vecQtot[7] = {TPCsignalTot,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
-      Double_t vecRatio[7] = {TPCsignalMax/TPCsignalTot,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
-      Double_t vecShortMediumRatio[7] = {TPCsignalShort/TPCsignalMedium,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
-      Double_t vecLongMediumRatio[7] = {TPCsignalLong/TPCsignalMedium,drift,snp,tgl,meanP,meanP/mass,NclsDeDx};
-          
-      fDeDxQmax->Fill(vecQmax); 
-      fDeDxQtot->Fill(vecQtot); 
-      fDeDxRatio->Fill(vecRatio); 
-      fDeDxShortMediumRatio->Fill(vecShortMediumRatio); 
-      fDeDxLongMediumRatio->Fill(vecLongMediumRatio);       
-      
-    }
-    
-  }
-  
-  
+      // dEdx in region
+      //
+      Float_t dEdxTot[5],dEdxTotFull[5];
+      Float_t dEdxMax[5],dEdxMaxFull[5];
+      Int_t   ncl[5];
+      for (Int_t itype=0; itype<5;itype++){
+       Int_t row0=0, row1 =159;
+       if (itype==1) {row0=0;      row1 = 63;};
+       if (itype==2) {row0= 64;    row1=63+64;}
+       if (itype==3) {row0= 63+64; row1=159;}
+       dEdxTot[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+       dEdxMax[itype]= (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+       // non trucated dEdx
+       dEdxTotFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,0,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+       dEdxMaxFull[itype]= (1/fMIP)*seed->CookdEdxNorm(0.0,0.99,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
+
+
+       ncl[itype]=TMath::Nint(seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,row0,row1,fUseShapeNorm,fUsePosNorm,fUsePadNorm,2));
+      }
+      //
+      //
+      //
+      Float_t wmeanTot=0, wmeanMax=0, sumW=0;
+      Double_t length[3] = {0.75,1,1.5};
+
+      for (Int_t ipad=0; ipad<3; ipad++){
+       if (ncl[1+ipad]<3) continue;
+       Double_t weight = Double_t(ncl[1+ipad])*TMath::Sqrt(length[ipad]);
+       wmeanTot+=dEdxTot[1+ipad]*weight;
+       wmeanMax+=dEdxMax[1+ipad]*weight;
+       sumW+=weight;
+      }
+      if (sumW>0){
+       dEdxTot[4]= wmeanTot/sumW;
+       dEdxMax[4]= wmeanMax/sumW;      
+      }
+      for (Int_t itype=0;itype<5;itype++){
+       //
+       Float_t snp=(TMath::Abs(snpIn)+TMath::Abs(snpOut))*0.5;
+       Float_t tgl=(TMath::Abs(tglIn)+TMath::Abs(tglOut))*0.5;
+       Float_t drift = (driftIn+driftOut)*0.5;
+       if (itype==1) {snp = TMath::Abs(snpIn); tgl = TMath::Abs(tglIn); drift= driftIn;};
+       if (itype==3) {snp = TMath::Abs(snpOut); tgl = TMath::Abs(tglOut);drift=driftOut;};
+       if (ncl[itype]<kMinClustersTracklet) continue;
+       Float_t deltaL = TMath::Sqrt(1+snp*snp+tgl*tgl);
+       //
+       Double_t vecQmax[8] = {dEdxMax[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
+       Double_t vecQtot[8] = {dEdxTot[itype],drift,snp,tgl,meanP,meanP/mass,NclsDeDx, itype};
+       //
+       //
+       //
+       Double_t ratioMaxTot           = (dEdxTot[itype]>0)  ? dEdxMax[itype]/dEdxTot[itype]:0;
+       Double_t ratioTrackletTot      = (dEdxTot[0]>0)      ? dEdxTot[itype]/dEdxTot[0]:0;
+       Double_t ratioTrackletMax      = (dEdxMax[0]>0)      ? dEdxMax[itype]/dEdxMax[0]:0;
+       Double_t ratioTrackletTruncTot = (dEdxTot[0]>0)      ? dEdxTotFull[itype]/dEdxTot[itype]:0;
+       Double_t ratioTrackletTruncMax = (dEdxMax[0]>0)      ? dEdxMaxFull[itype]/dEdxMax[itype]:0;
+
+       Double_t vecRatioMaxTot[8]      = {ratioMaxTot,      drift,snp,tgl,dEdxTot[0],  dEdxTot[0]*deltaL,NclsDeDx,itype};
+       Double_t vecRatioTrackletTot[8] = {ratioTrackletTot, drift,snp,tgl,dEdxTot[0],  dEdxTot[0]*deltaL,NclsDeDx,itype};      
+       Double_t vecRatioTrackletMax[8] = {ratioTrackletMax, drift,snp,tgl,dEdxMax[0],  dEdxMax[0]*deltaL,NclsDeDx,itype};      
+       Double_t vecRatioTrackletTruncTot[8] = {ratioTrackletTruncTot, drift,snp,tgl,dEdxTot[0],  dEdxTot[0]*deltaL,NclsDeDx,itype};    
+       Double_t vecRatioTrackletTruncMax[8] = {ratioTrackletTruncMax, drift,snp,tgl,dEdxMax[0],  dEdxMax[0]*deltaL,NclsDeDx,itype};    
+       fDeDxQmax->Fill(vecQmax); 
+       fDeDxQtot->Fill(vecQtot); 
+       fDeDxRatioMaxTot->Fill(vecRatioMaxTot); 
+       fDeDxRatioQmax->Fill(vecRatioTrackletTot); 
+       fDeDxRatioQtot->Fill(vecRatioTrackletMax); 
+       fDeDxRatioTruncQmax->Fill(vecRatioTrackletTruncTot); 
+       fDeDxRatioTruncQtot->Fill(vecRatioTrackletTruncMax); 
+       //
+       TTreeSRedirector * cstream =  GetDebugStreamer();
+       if (cstream){
+         TVectorD vQT(9,vecQtot);
+         TVectorD vQM(9,vecQmax);
+         TVectorD vQMT(9, vecRatioMaxTot);
+         TVectorD vQRT(9,vecRatioTrackletTot);
+         TVectorD vQRM(9,vecRatioTrackletMax);
+         TVectorD vQRTT(9,vecRatioTrackletTruncTot);
+         TVectorD vQRTM(9,vecRatioTrackletTruncMax);
+         (*cstream) << "dEdx" <<
+           "run="<<fRun<<              //  run number
+           "event="<<fEvent<<          //  event number
+           "time="<<fTime<<            //  time stamp of event
+           "trigger="<<fTrigger<<      //  trigger
+           "mag="<<fMagF<<             //  magnetic field            
+           "track.="<<seed<<           //  original seed
+           //
+           "vQT.="<<&vQT<<             // trunc mean - total charge
+           "vQM.="<<&vQM<<             // trunc mean - max charge 
+           //
+           "vQMT.="<<&vQMT<<           // ratio max/tot      
+           "vQRT.="<<&vQRT<<           // ratio tracklet/track - total charge
+           "vQRM.="<<&vQRM<<           // ratio tracklet/track - max charge
+           "vQRTT.="<<&vQRTT<<         // ratio trunc/full     - total charge
+           "vQRTM.="<<&vQRTM<<         // ratio trunc/full     - max charge
+           "\n";
+       }
+      }   
+    }    
+  }  
 }
 
 
@@ -338,9 +444,11 @@ Long64_t AliTPCcalibPID::Merge(TCollection *li) {
     //
     if (cal->GetHistQmax()) fDeDxQmax->Add(cal->GetHistQmax());
     if (cal->GetHistQtot()) fDeDxQtot->Add(cal->GetHistQtot());
-    if (cal->GetHistRatio()) fDeDxRatio->Add(cal->GetHistRatio());
-    if (cal->GetHistShortMediumRatio()) fDeDxShortMediumRatio->Add(cal->GetHistShortMediumRatio());
-    if (cal->GetHistLongMediumRatio()) fDeDxLongMediumRatio->Add(cal->GetHistLongMediumRatio());
+    if (cal->GetHistRatioMaxTot()) fDeDxRatioMaxTot->Add(cal->GetHistRatioMaxTot());
+    if (cal->GetHistRatioQmax()) fDeDxRatioQmax->Add(cal->GetHistRatioQmax());
+    if (cal->GetHistRatioQtot()) fDeDxRatioQtot->Add(cal->GetHistRatioQtot());
+    if (cal->GetHistRatioTruncQmax()) fDeDxRatioTruncQmax->Add(cal->GetHistRatioTruncQmax());
+    if (cal->GetHistRatioTruncQtot()) fDeDxRatioTruncQtot->Add(cal->GetHistRatioTruncQtot());
   }
   
   return 0;
@@ -349,17 +457,6 @@ Long64_t AliTPCcalibPID::Merge(TCollection *li) {
 
 
 
-
-void AliTPCcalibPID::MakeReport() {
-  //
-  // 1. standard dEdx picture
-  TCanvas *cDeDx = new TCanvas("cDeDx","cDeDx");
-  GetHistQmax()->Projection(0,4)->Draw("colz");
-  return;
-}
-
-
-
 void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
 
   // Method for the correct logarithmic binning of histograms
@@ -382,3 +479,377 @@ void AliTPCcalibPID::BinLogX(THnSparse *h, Int_t axisDim) {
   
 }
 
+
+
+
+
+void AliTPCcalibPID::MakeReport(const char *outputpath) {
+  //
+  // Make a standard QA plots
+  //
+  for (Int_t ipad=0;ipad<4;ipad++){
+    DrawRatioTot(ipad,outputpath);
+    DrawRatioMax(ipad,outputpath);
+  }
+  DrawRatiodEdx(0.5,3,outputpath);
+  DrawResolBGQtot(140,160,outputpath);
+  DrawResolBGQmax(140,160,outputpath);
+  return;
+}
+
+void  AliTPCcalibPID::DrawRatioTot(Int_t ipad, const char* outputpath){
+  //
+  // Draw default ratio histogram for given pad type
+  // ipad - 0 - short pads
+  //        1 - medium pads
+  //        2 - long pads
+  //
+  //  Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+  Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+  AliTPCcalibPID * pid = this;
+  TCanvas *canvas= new TCanvas(Form("QtotRatio%d)",ipad),Form("QtotRatio%d)",ipad));
+  canvas->Divide(3,2);
+  pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+  canvas->cd(1);
+  TH1 * his0 = pid->GetHistRatioQtot()->Projection(0);
+  his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
+  his0->SetYTitle("");  
+  his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  his0->Draw();
+  //
+  canvas->cd(2);
+  TH1 * hprofz = (TH1*) (pid->GetHistRatioQtot()->Projection(0,1)->ProfileX());
+  hprofz->SetXTitle("drift length");
+  hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofz->SetMarkerStyle(kmimarkers[0]);
+  hprofz->SetMaximum(1.1);
+  hprofz->SetMinimum(0.9);
+  hprofz->Draw();
+  //
+  canvas->cd(3);
+  TH1 * hprofphi = (TH1*) (pid->GetHistRatioQtot()->Projection(0,2)->ProfileX());
+  hprofphi->SetXTitle("sin(#phi)");
+  hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));    
+  hprofphi->SetMarkerStyle(kmimarkers[0]);  
+  hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofphi->SetMaximum(1.1);
+  hprofphi->SetMinimum(0.9);
+  hprofphi->Draw();
+  //
+  canvas->cd(4);
+  TH1 * hproftheta = (TH1*) (pid->GetHistRatioQtot()->Projection(0,3)->ProfileX());
+  hproftheta->SetXTitle("tan(#theta)");
+  hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hproftheta->SetMarkerStyle(kmimarkers[0]);
+  hproftheta->SetMaximum(1.1);
+  hproftheta->SetMinimum(0.9);
+  hproftheta->Draw();
+  //
+  canvas->cd(5);
+  TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQtot()->Projection(0,4)->ProfileX());
+  hprofdedx->SetXTitle("dEdx_{track}");
+  hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofdedx->SetMarkerStyle(kmimarkers[0]);
+  hprofdedx->SetMaximum(1.1);
+  hprofdedx->SetMinimum(0.9);
+  hprofdedx->Draw();
+
+  canvas->cd(6);
+  TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
+  hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
+  hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofdedxL->SetMarkerStyle(kmimarkers[0]);
+  hprofdedxL->SetMaximum(1.1);
+  hprofdedxL->SetMinimum(0.9);
+  hprofdedxL->Draw();
+
+
+
+  canvas->SaveAs(Form("%s/QtotRatioType%d.eps",outputpath,ipad));
+  canvas->SaveAs(Form("%s/QtotRatioType%d.gif",outputpath,ipad));
+}
+
+void  AliTPCcalibPID::DrawRatioMax(Int_t ipad, const char* outputpath){
+  //
+  // Draw default ration histogram for given pad type
+  // ipad - 0 - short pads
+  //        1 - medium pads
+  //        2 - long pads
+  //
+  //  Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+  Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+  AliTPCcalibPID * pid = this;
+  TCanvas *canvas= new TCanvas(Form("QmaxRatio%d)",ipad),Form("QmaxRatio%d)",ipad));
+  canvas->Divide(3,2);
+  pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+  canvas->cd(1);
+  TH1 * his0 = pid->GetHistRatioQmax()->Projection(0);
+  his0->SetXTitle("dEdx_{tracklet}/dEdx_{track}");
+  his0->SetYTitle("");  
+  his0->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  his0->Draw();
+  //
+  canvas->cd(2);
+  TH1 * hprofz = (TH1*) (pid->GetHistRatioQmax()->Projection(0,1)->ProfileX());
+  hprofz->SetXTitle("drift length");
+  hprofz->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofz->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofz->SetMarkerStyle(kmimarkers[0]);
+  hprofz->SetMaximum(1.1);
+  hprofz->SetMinimum(0.9);
+  hprofz->Draw();
+  //
+  canvas->cd(3);
+  TH1 * hprofphi = (TH1*) (pid->GetHistRatioQmax()->Projection(0,2)->ProfileX());
+  hprofphi->SetXTitle("sin(#phi)");
+  hprofphi->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));    
+  hprofphi->SetMarkerStyle(kmimarkers[0]);  
+  hprofphi->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofphi->SetMaximum(1.1);
+  hprofphi->SetMinimum(0.9);
+  hprofphi->Draw();
+  //
+  canvas->cd(4);
+  TH1 * hproftheta = (TH1*) (pid->GetHistRatioQmax()->Projection(0,3)->ProfileX());
+  hproftheta->SetXTitle("tan(#theta)");
+  hproftheta->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hproftheta->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hproftheta->SetMarkerStyle(kmimarkers[0]);
+  hproftheta->SetMaximum(1.1);
+  hproftheta->SetMinimum(0.9);
+  hproftheta->Draw();
+
+  canvas->cd(5);
+  TH1 * hprofdedx = (TH1*) (pid->GetHistRatioQmax()->Projection(0,4)->ProfileX());
+  hprofdedx->SetXTitle("dEdx_{track}");
+  hprofdedx->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofdedx->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofdedx->SetMarkerStyle(kmimarkers[0]);
+  hprofdedx->SetMaximum(1.1);
+  hprofdedx->SetMinimum(0.9);
+  hprofdedx->Draw();
+
+  canvas->cd(6);
+  TH1 * hprofdedxL = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
+  hprofdedxL->SetXTitle("dEdx_{track}#Delta_{x}");
+  hprofdedxL->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+  hprofdedxL->SetTitle(Form("dEdx_{tracklet}/dEdx_{track} type %d",ipad));
+  hprofdedxL->SetMarkerStyle(kmimarkers[0]);
+  hprofdedxL->SetMaximum(1.1);
+  hprofdedxL->SetMinimum(0.9);
+  hprofdedxL->Draw();
+
+
+  canvas->SaveAs(Form("%s/QmaxRatioType%d.eps",outputpath,ipad));
+  canvas->SaveAs(Form("%s/QmaxRatioType%d.gif",outputpath,ipad));
+}
+
+
+
+void AliTPCcalibPID::DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath){
+  //
+  //
+  //
+  //
+  //  Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
+  Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+  AliTPCcalibPID * pid = this;
+  TCanvas *canvas= new TCanvas("QRatiodEdx","QRatiodEdx");
+  canvas->Divide(2,4);
+  canvas->SetLogx(kTRUE);
+  TH1 * hprofP=0;
+  for (Int_t ipad=0;ipad<4;ipad++){
+    pid->GetHistRatioQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+    pid->GetHistRatioQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+    pid->GetHistRatioQmax()->GetAxis(5)->SetRangeUser(demin,demax);
+    pid->GetHistRatioQtot()->GetAxis(5)->SetRangeUser(demin,demax);
+
+    canvas->cd(ipad*2+1)->SetLogx(kFALSE);
+    hprofP  = (TH1*) (pid->GetHistRatioQmax()->Projection(0,5)->ProfileX());
+    hprofP->SetXTitle("dEdx_{track}");
+    hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+    hprofP->SetTitle(Form("Q_{max} dEdx_{tracklet}/dEdx_{track} type %d",0));
+    hprofP->SetMarkerStyle(kmimarkers[0]);     
+    hprofP->SetMaximum(1.1);
+    hprofP->SetMinimum(0.9);
+    hprofP->Draw();
+    // pad  Tot
+    canvas->cd(ipad*2+2)->SetLogx(kFALSE);
+    hprofP  = (TH1*) (pid->GetHistRatioQtot()->Projection(0,5)->ProfileX());
+    hprofP->SetXTitle("dEdx_{track}");
+    hprofP->SetYTitle("dEdx_{tracklet}/dEdx_{track}");
+    hprofP->SetTitle(Form("Q_{tot} dEdx_{tracklet}/dEdx_{track} type %d",0));
+    hprofP->SetMarkerStyle(kmimarkers[0]);     
+    hprofP->SetMaximum(1.1);
+    hprofP->SetMinimum(0.9);
+    hprofP->Draw();
+  }
+  //
+  //
+  canvas->SaveAs(Form("%s/QratiodEdx%d.eps",outputpath));
+  canvas->SaveAs(Form("%s/QratiodEdx%d.gif",outputpath));
+}
+
+
+
+void AliTPCcalibPID::DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, const char *outputpath){
+  //
+  // 
+  //
+  AliTPCcalibPID * pid = this;
+  Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+
+  TObjArray arr;
+  TH2 * his   =0; 
+  TH1 * hmean =0;
+  TH1 * hsigma=0;
+  //
+  // set cut
+  pid->GetHistQtot()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
+  pid->GetHistQtot()->GetAxis(5)->SetRangeUser(1,10000);
+  TCanvas *canvas= new TCanvas("dEdxResolQ_{Tot}","dEdxResolQ_{Tot}");
+  canvas->Divide(2,3);
+  //
+  //
+  //
+  for (Int_t ipad=0;ipad<5;ipad++){
+    canvas->cd(1+ipad)->SetLogx(kTRUE);
+    if (ipad<4) pid->GetHistQtot()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+    if (ipad==4) pid->GetHistQtot()->GetAxis(7)->SetRange(1,1);
+    his = (TH2*)(pid->GetHistQtot()->Projection(0,5));
+    his->FitSlicesY(0,0,-1,0,"QNR",&arr);
+    hmean = (TH1*)arr.At(1);
+    hsigma = (TH1*)arr.At(2)->Clone();    
+    hsigma->Divide(hmean);
+    arr.SetOwner(kTRUE);
+    arr.Delete();
+    delete his;
+
+    hsigma->SetXTitle("#beta#gamma");
+    hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
+    hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
+    hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Pad %d",ipad));
+    if (ipad==4) {
+      hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
+      hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{tot} Full"));
+    }
+    hsigma->SetMarkerStyle(kmimarkers[0]);
+    hsigma->SetMaximum(0.15);
+    hsigma->SetMinimum(0.0);
+    hsigma->Draw();
+  }
+  canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
+  canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
+}
+
+void AliTPCcalibPID::DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, const char *outputpath){
+  //
+  // Int_t minClusters=140;  Int_t maxClusters=200; const char *outputpath="picPID06"
+  //
+  Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
+  AliTPCcalibPID * pid = this;
+  TObjArray arr;
+  TH2 * his   =0; 
+  TH1 * hmean =0;
+  TH1 * hsigma=0;
+  //
+  // set cut
+  pid->GetHistQmax()->GetAxis(6)->SetRangeUser(minClusters,maxClusters);
+  pid->GetHistQmax()->GetAxis(5)->SetRangeUser(1,10000);
+  TCanvas *canvas= new TCanvas("dEdxResolQ_{Max}","dEdxResolQ_{Max}");
+  canvas->Divide(2,3);
+  //
+  //
+  //
+  for (Int_t ipad=0;ipad<5;ipad++){
+    canvas->cd(1+ipad)->SetLogx(kTRUE);
+    if (ipad<4) pid->GetHistQmax()->GetAxis(7)->SetRange(ipad+2,ipad+2);
+    if (ipad==4) pid->GetHistQmax()->GetAxis(7)->SetRange(1,1);
+    his = (TH2*)(pid->GetHistQmax()->Projection(0,5));
+    his->FitSlicesY(0,0,-1,0,"QNR",&arr);
+    hmean = (TH1*)arr.At(1);
+    hsigma = (TH1*)arr.At(2)->Clone();
+    hsigma->Divide(hmean);
+    arr.SetOwner(kTRUE);
+    arr.Delete();
+    delete his;
+    hsigma->SetXTitle("#beta#gamma");
+    hsigma->SetYTitle("#sigma_{dEdx}/dEdx");
+    hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
+    hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Pad %d",ipad));
+    if (ipad==4){
+      hsigma->SetTitle(Form("#sigma_{dEdx}/dEdx_{max} Full"));
+      hsigma->SetName(Form("#sigma_{dEdx}/dEdx_{max} Full"));
+    }
+    hsigma->SetMarkerStyle(kmimarkers[0]);
+    hsigma->SetMaximum(0.15);
+    hsigma->SetMinimum(0.0);
+    hsigma->Draw();
+  }
+  canvas->SaveAs(Form("%s/dEdxResolMax.eps",outputpath));
+  canvas->SaveAs(Form("%s/dEdxResolMax.gif",outputpath));
+}
+
+
+
+void AliTPCcalibPID::DumpTree(THnSparse * hndim, const char * outname){
+  //
+  // make a tree
+  // the tree will be written to the file - outname
+  //
+  TTreeSRedirector *pcstream = new TTreeSRedirector(outname);
+  //
+  //
+  Double_t position[10];
+  Double_t value; 
+  Int_t *bins = new Int_t[10];
+  //
+  //
+  Int_t ndim = hndim->GetNdimensions();
+  //
+  for (Long64_t i = 0; i < hndim->GetNbins(); ++i) {
+    value = hndim->GetBinContent(i, bins);
+    for (Int_t idim = 0; idim < ndim; idim++) {
+      position[idim] = hndim->GetAxis(idim)->GetBinCenter(bins[idim]);
+    }      
+    Double_t ty = TMath::Tan(TMath::ASin(position[2]));
+    Double_t tz = position[3]*TMath::Sqrt(1+ty*ty);
+    //
+    (*pcstream)<<"Dump"<<
+      "bincont="<<value<<      // bin content
+      "val="<<position[0]<<    // parameter difference
+      "dr="<<position[1]<<     //drift
+      "ty="<<ty<<              //phi
+      "tz="<<tz<<              //theta      
+      "p2="<<position[2]<<      //p2
+      "p3="<<position[3]<<     //p3
+      "p="<<position[4]<<      //p
+      "bg="<<position[5]<<     //bg
+      "ncl="<<position[6]<<    //ncl
+      "type="<<position[7]<<   //tracklet
+      "tot="<<position[8]<<    //is tot 
+      "\n";
+  }
+  delete pcstream;
+}
+
+
+void AliTPCcalibPID::DumpTrees(){
+  //
+  // dump the content of histogram to the tree
+  //
+  AliTPCcalibPID *pid = this;
+  DumpTree(pid->GetHistQtot(),"dumpQtot.root");
+  DumpTree(pid->GetHistQmax(),"dumpQmax.root");
+  DumpTree(pid->GetHistRatioQtot(),"dumpRatioQtot.root");
+  DumpTree(pid->GetHistRatioQmax(),"dumpRatioQmax.root");
+}
+
+
+
index 0dee910..dbbc689 100644 (file)
@@ -27,8 +27,12 @@ public:
   virtual void           Process(AliESDEvent *event);
   virtual Long64_t       Merge(TCollection *li);
   virtual void           Analyze();
-  void                   MakeReport();
-  //
+  void                   MakeReport(const char * outputpath);
+  void                   DrawRatioTot(Int_t ipad, const char* outputpath);
+  void                   DrawRatioMax(Int_t ipad, const char* outputpath);
+  void                   DrawRatiodEdx(Float_t demin, Float_t demax, const char* outputpath);
+  void                   DrawResolBGQtot(Int_t minClusters, Int_t maxClusters, const char *outputpath); //
+  void                   DrawResolBGQmax(Int_t minClusters, Int_t maxClusters, const char *outputpath);
   //
   TH1F   *          GetHistNTracks(){return fHistNTracks;};
   TH1F   *          GetHistClusters(){return fClusters;};
@@ -37,22 +41,26 @@ public:
   //
   THnSparseS *      GetHistQmax(){return fDeDxQmax;};
   THnSparseS *      GetHistQtot(){return fDeDxQtot;};
-  THnSparseS *      GetHistRatio(){return fDeDxRatio;};
-  THnSparseS *      GetHistShortMediumRatio(){return fDeDxShortMediumRatio;};
-  THnSparseS *      GetHistLongMediumRatio(){return fDeDxLongMediumRatio;};
+  THnSparseS *      GetHistRatioMaxTot(){return fDeDxRatioMaxTot;};
+  THnSparseS *      GetHistRatioQmax(){return fDeDxRatioQmax;};
+  THnSparseS *      GetHistRatioQtot(){return fDeDxRatioQtot;};
+  THnSparseS *      GetHistRatioTruncQmax(){return fDeDxRatioTruncQmax;};
+  THnSparseS *      GetHistRatioTruncQtot(){return fDeDxRatioTruncQtot;};
   //
   void SetMIPvalue(Float_t mip){fMIP = mip;};
   void SetLowerTrunc(Float_t lowerTrunc){fLowerTrunc = lowerTrunc;};
   void SetUpperTrunc(Float_t upperTrunc){fUpperTrunc = upperTrunc;};
   void SetUseShapeNorm(Bool_t useShapeNorm){fUseShapeNorm = useShapeNorm;};
-  void SetUsePosNorm(Bool_t usePosNorm){fUsePosNorm = usePosNorm;};
+  void SetUsePosNorm(Int_t usePosNorm){fUsePosNorm = usePosNorm;};
   void SetPadNorm(Int_t padNorm){fUsePadNorm = padNorm;};
   void SetIsCosmic(Bool_t isCosmic){fIsCosmic = isCosmic;};
   //
   //
   static void       BinLogX(THnSparse * h, Int_t axisDim);   // method for correct histogram binning
-
-
+  void DumpTree(THnSparse * hndim, const char * outname);
+  void DumpTrees();
+  void     Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
+  void     Process(AliTPCseed *track){return AliTPCcalibBase::Process(track);}
 private:
   //
   // parameter specifications
@@ -61,7 +69,7 @@ private:
   Float_t fLowerTrunc;
   Float_t fUpperTrunc;
   Bool_t  fUseShapeNorm;
-  Bool_t  fUsePosNorm;
+  Int_t  fUsePosNorm;
   Int_t   fUsePadNorm;
   //
   Bool_t  fIsCosmic;
@@ -75,9 +83,15 @@ private:
   //
   THnSparseS * fDeDxQmax;               //  histogram which shows dEdx (Qmax) as a function of z,sin(phi),tan(theta),p,betaGamma
   THnSparseS * fDeDxQtot;               //  histogram which shows dEdx (Qtot) as a function of z,sin(phi),tan(theta),p,betaGamma
-  THnSparseS * fDeDxRatio;              //  histogram which shows dEdx ratio (Qmax/Qtot) as a function of z,sin(phi),tan(theta),p,betaGamma
-  THnSparseS * fDeDxShortMediumRatio;   //  histogram which shows dEdx ratio (QmaxShort/QmaxMedium) as a function of z,sin(phi),tan(theta),p,betaGamma
-  THnSparseS * fDeDxLongMediumRatio;    //  histogram which shows dEdx ratio (QmaxLong/QmaxMedium) as a function of z,sin(phi),tan(theta),p,betaGamma
+  //
+  // ratio histograms
+  //
+  THnSparseS * fDeDxRatioMaxTot;              //  histogram which shows dEdx ratio (Qmax/Qtot) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+  THnSparseS * fDeDxRatioQmax;   // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+  THnSparseS * fDeDxRatioQtot;   // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+  THnSparseS * fDeDxRatioTruncQtot;   // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+  THnSparseS * fDeDxRatioTruncQmax;   // dEdx ratio (tracklet/track) as a function of z,sin(phi),tan(theta),dEdx,dEdx*dl
+
   //
   AliTPCcalibPID(const AliTPCcalibPID&); 
   AliTPCcalibPID& operator=(const AliTPCcalibPID&); 
index 1748e3b..297fac8 100644 (file)
   dy =          +dr(r,z)*tan(phi)
 
    
-
-  //  
+  .x ~/NimStyle.C
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libTPCcalib");
+  TFile fcalib("CalibObjects.root");
+  TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+  AliTPCcalibUnlinearity * calibUnlin = ( AliTPCcalibUnlinearity *)array->FindObject("calibUnlinearity");
+  //
+  
 */
 
 
@@ -76,7 +82,6 @@
 #include "AliTPCTracklet.h"
 #include "TTimeStamp.h"
 #include "AliTPCcalibDB.h"
-#include "AliTPCcalibLaser.h"
 #include "AliDCSSensorArray.h"
 #include "AliDCSSensor.h"
 #include "TLinearFitter.h"
 
 
 #include "AliTPCcalibUnlinearity.h"
+#include "AliTPCPointCorrection.h"
 
-AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::fgInstance = 0;
 
 ClassImp(AliTPCcalibUnlinearity)
 
 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
   AliTPCcalibBase(),
-  fDiffHistoLine(0),
-  fDiffHistoPar(0),
+  fInit(kFALSE),
+  fDiffHistoLineY(0),
+  fDiffHistoParY(0),
+  fDiffHistoLineZ(0),
+  fDiffHistoParZ(0),
   fFittersOutR(38),
   fFittersOutZ(38),
   fParamsOutR(38),
   fParamsOutZ(38),
   fErrorsOutR(38),
-  fErrorsOutZ(38)
+  fErrorsOutZ(38),
+  fDistRPHIPlus(74),
+  fDistRPHIMinus(74),
+  fFitterQuadrantY(16*38),
+  fFitterQuadrantPhi(16*38)
 {
   //
   // Defualt constructor
@@ -108,42 +120,39 @@ AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
 
 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
   AliTPCcalibBase(name,title),
-  fDiffHistoLine(0),
-  fDiffHistoPar(0),  
+  fInit(kFALSE),
+  fDiffHistoLineY(0),
+  fDiffHistoParY(0),  
+  fDiffHistoLineZ(0),
+  fDiffHistoParZ(0),  
   fFittersOutR(38),
   fFittersOutZ(38),
   fParamsOutR(38),
   fParamsOutZ(38),
   fErrorsOutR(38),
-  fErrorsOutZ(38)
+  fErrorsOutZ(38),
+  fDistRPHIPlus(74),
+  fDistRPHIMinus(74),
+  fFitterQuadrantY(16*38),
+  fFitterQuadrantPhi(16*38)
 {
   //
   // Non default constructor
   //
-  MakeFitters();
+  MakeHisto();
 }
 
 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
   //
   //
   //
-  if (fDiffHistoLine) delete fDiffHistoLine;
-  if (fDiffHistoPar)  delete fDiffHistoPar;
+  if (fDiffHistoLineZ) delete fDiffHistoLineY;
+  if (fDiffHistoParZ)  delete fDiffHistoParY;
+  if (fDiffHistoLineZ) delete fDiffHistoLineZ;
+  if (fDiffHistoParZ)  delete fDiffHistoParZ;
 }
 
 
-AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::Instance()
-{
-  //
-  // Singleton implementation
-  // Returns an instance of this class, it is created if neccessary
-  //
-  if (fgInstance == 0){
-    fgInstance = new AliTPCcalibUnlinearity();
-  }
-  return fgInstance;
-}
-
 
 
 
@@ -151,9 +160,11 @@ void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
   //
   // 
   //  
+  if (!fInit) {
+    Init();   // work around for PROOF - initialize firs time used
+  }
   const Int_t  kdrow=10;
   const Int_t kMinCluster=35;
-  if (!fDiffHistoLine) MakeHisto();    
   if (!seed) return;
   if (TMath::Abs(fMagF)>0.01) return;   // use only data without field
   Int_t counter[72];
@@ -167,6 +178,9 @@ void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
   //
   for (Int_t is=0; is<72;is++){
     if (counter[is]>kMinCluster) ProcessDiff(seed, is);
+    if (counter[is]>kMinCluster) {
+      AlignOROC(seed, is);
+    }
   }
 }
 
@@ -177,6 +191,7 @@ void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
   //
   const  Double_t kXmean=134;
   const Int_t     kdrow=10;
+  const Float_t kSagitaCut = 1;
   static TLinearFitter fy1(2,"pol1");
   static TLinearFitter fy2(3,"pol2");
   static TLinearFitter fz1(2,"pol1");
@@ -214,7 +229,6 @@ void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
     AliTPCclusterMI *c=track->GetClusterPointer(irow);
     if (!c) continue;
     if (c->GetDetector()!=isec) continue;
-    if (c->GetType()<0) continue;
     Double_t dx = c->GetX()-kXmean;
     Double_t y1 = py1[0]+py1[1]*dx;
     Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
@@ -222,29 +236,24 @@ void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
     Double_t z1 = pz1[0]+pz1[1]*dx;
     Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
     //
+    Double_t edgeMinus= -y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
+    Double_t edgePlus =  y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
+    Int_t    npads = AliTPCROC::Instance()->GetNPads(isec,c->GetRow());
+    Float_t    dpad  = TMath::Min(c->GetPad(), npads-1- c->GetPad());
+    Float_t dy =c->GetY()-y1;
+    //
+    // Corrections 
+    //    
+    AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
+    Double_t corrclY=0, corrtrY=0, corrR=0;
+    corrclY =  corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+                                      c->GetY(),c->GetY(), c->GetZ(), py1[1],c->GetMax(),2.5); 
+    corrtrY =  corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+                                      c->GetY(),y1, c->GetZ(), py1[1], c->GetMax(),2.5);     
+    corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,c->GetX(),c->GetY(),c->GetZ(),isec);
+    //
     //
-    Double_t x[10];
-    x[0]=isec;
-    x[1]=irow;
-    x[2]=c->GetY();
-    x[3]=c->GetZ();
-    x[4]=py1[1];
-    x[5]=pz1[1];
-    x[6]=py2[2]*150*150/4; // sagita 150 cm
     //
-    x[7]=c->GetY()-y1;
-    x[8]=c->GetZ()-z1;
-    fDiffHistoLine->Fill(x);
-    x[7]=c->GetY()-y2;
-    x[8]=c->GetZ()-z2;
-    fDiffHistoPar->Fill(x);   
-
-    if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(pz2[2]*150*150/4)<1){
-      // Apply sagita cut
-      AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
-              py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
-    }
-            
     if (fStreamLevel>1){
       TTreeSRedirector *cstream = GetDebugStreamer();
       if (cstream){
@@ -254,7 +263,14 @@ void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
          "time="<<fTime<<            //  time stamp of event
          "trigger="<<fTrigger<<      //  trigger
          "mag="<<fMagF<<             //  magnetic field          
-         "isec="<<isec<<
+         "isec="<<isec<<             // current sector
+         "npads="<<npads<<           // number of pads at given sector
+         "dpad="<<dpad<<             // distance to edge pad
+         //
+         "crR="<<corrR<<              // radial correction
+         "cclY="<<corrclY<<            // edge effect correction cl
+         "ctrY="<<corrtrY<<            // edge effect correction using track
+         //
          "Cl.="<<c<<
          "y1="<<y1<<
          "y2="<<y2<<
@@ -264,13 +280,269 @@ void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
          "py2.="<<&py2<<
          "pz1.="<<&pz1<<
          "pz2.="<<&pz2<<
+         "eP="<<edgePlus<<
+         "eM="<<edgeMinus<<
          "\n";
       }
     }
+    if (TMath::Min(edgeMinus,edgePlus)<6){
+      AddPointRPHI(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
+      TTreeSRedirector *cstream = GetDebugStreamer();
+      if (TMath::Abs(dy)<2 && cstream){
+       (*cstream)<<"rphi"<<
+         "isec="<<isec<<             // current sector
+         "npads="<<npads<<           // number of pads at given sector
+         "dpad="<<dpad<<             // distance to edge pad
+         "eP="<<edgePlus<<           // distance to edge 
+         "eM="<<edgeMinus<<          // distance to edge minus
+         //
+         "crR="<<corrR<<              // radial correction
+         "cclY="<<corrclY<<            // edge effect correction cl
+         "ctrY="<<corrtrY<<            // edge effect correction using track
+         //
+         "dy="<<dy<<                 // dy
+         "Cl.="<<c<<
+         "y1="<<y1<<
+         "y2="<<y2<<
+         "z1="<<z1<<
+         "z2="<<z2<<
+         "py1.="<<&py1<<
+         "pz1.="<<&pz1<<
+         "\n";
+      }
+    }
+    if (c->GetType()<0) continue;  // don't use edge rphi cluster
+    //
+    //
+    Double_t x[10];
+    x[0]=c->GetDetector();
+    x[1]=c->GetRow();
+    x[2]=c->GetY()/c->GetX();
+    x[3]=c->GetZ();
+    //
+    // apply sagita cut
+    //
+    if (TMath::Abs(py2[2]*150*150/4)<kSagitaCut && 
+       TMath::Abs(pz2[2]*150*150/4)<kSagitaCut){
+      //
+      x[4]=py1[1];
+      x[5]=c->GetY()-y1;
+      fDiffHistoLineY->Fill(x);
+      x[5]=c->GetY()-y2;
+      //fDiffHistoParY->Fill(x);
+      //
+      x[4]=pz1[1];
+      x[5]=c->GetY()-z1;
+      fDiffHistoLineZ->Fill(x);
+      x[5]=c->GetY()-z2;
+      //fDiffHistoParZ->Fill(x);      
+      // Apply sagita cut
+      AddPoint(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
+    }
+            
   }
 }
 
 
+void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
+  //
+  // fit the tracklet in OROC sepatately in Quadrant
+  //
+  //
+  if (isec<36) return; 
+  const Int_t kMinClusterF=35;
+  const Int_t kMinClusterQ=10;
+  //
+  const  Int_t     kdrow1 =8;       // rows to skip at the end      
+  const  Int_t     kdrow0 =2;        // rows to skip at beginning  
+  const  Float_t   kedgey=3;
+  const  Float_t   kMaxDist=0.5;
+  const  Float_t   kMaxCorrY=0.1;
+  const  Float_t   kPRFWidth = 0.6;   //cut 2 sigma of PRF
+  //
+  //
+  AliTPCROC * roc = AliTPCROC::Instance();
+  AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
+  const  Double_t kXmean=roc->GetPadRowRadii(isec,53);
+  //
+  // full track fit parameters
+  // 
+  TLinearFitter fyf(2,"pol1");
+  TLinearFitter fzf(2,"pol1");
+  TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
+  Int_t nf=0;
+  //
+  // make full fit as reference
+  //
+  for (Int_t iter=0; iter<2; iter++){
+    fyf.ClearPoints();
+    for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+      AliTPCclusterMI *c=track->GetClusterPointer(irow);
+      if (!c) continue;
+      if (c->GetDetector()!=isec) continue;
+      if (c->GetRow()<kdrow0) continue;
+      Double_t dx = c->GetX()-kXmean;
+      Double_t x[2]={dx, dx*dx};
+      if (iter==0 &&c->GetType()<0) continue;
+      if (iter==1){    
+       Double_t yfit  =  pyf[0]+pyf[1]*dx;
+       Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+       if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+       if (dedge<kedgey) continue;
+       Double_t corrtrY =  
+         corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+                                 c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+       if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+      }
+      fyf.AddPoint(x,c->GetY(),0.1);
+      fzf.AddPoint(x,c->GetZ(),0.1);
+    }
+    nf = fyf.GetNpoints();
+    if (nf<kMinClusterF) return;   // not enough points - skip 
+    fyf.Eval(); 
+    fyf.GetParameters(pyf); 
+    fyf.GetErrors(peyf);
+    fzf.Eval(); 
+    fzf.GetParameters(pzf); 
+    fzf.GetErrors(pezf);
+  }
+  //
+  // Make Fitters and params for 3 fitters
+  //
+  TLinearFitter *fitters[4];
+  Int_t npoints[4];
+  TVectorD params[4];
+  TVectorD errors[4];
+  Double_t chi2C[4];
+  for (Int_t i=0;i<4;i++) {
+    fitters[i] = new TLinearFitter(2,"pol1");
+    npoints[i]=0;
+    params[i].ResizeTo(2);
+    errors[i].ResizeTo(2);
+  }
+  //
+  // Update fitters
+  //
+  for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+    AliTPCclusterMI *c=track->GetClusterPointer(irow);
+    if (!c) continue;
+    if (c->GetDetector()!=isec) continue;
+    if (c->GetRow()<kdrow0) continue;      
+    Double_t dx = c->GetX()-kXmean;
+    Double_t x[2]={dx, dx*dx};
+    Double_t yfit  =  pyf[0]+pyf[1]*dx;
+    Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+    if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+    if (dedge<kedgey) continue;
+    Double_t corrtrY =  
+      corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+                             c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+    if (TMath::Abs(corrtrY)>kMaxCorrY) continue;  
+    if (dx<0){
+      if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
+      if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
+    }
+    if (dx>0){
+      if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
+      if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
+    }
+  }
+
+  for (Int_t i=0;i<4;i++) {
+    npoints[i] = fitters[i]->GetNpoints();
+    if (npoints[i]>=kMinClusterQ){
+      fitters[i]->Eval();
+      Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
+      chi2C[i]=chi2Fac;
+      fitters[i]->GetParameters(params[i]);
+      fitters[i]->GetErrors(errors[i]);
+      // renormalize errors
+      errors[i][0]*=chi2Fac;
+      errors[i][1]*=chi2Fac;
+    }
+  }
+  //
+  // Fill fitters
+  //
+  for (Int_t i0=0;i0<4;i0++){
+    for (Int_t i1=0;i1<4;i1++){
+      if (i0==i1) continue;
+      if(npoints[i0]<kMinClusterQ) continue;
+      if(npoints[i1]<kMinClusterQ) continue;
+      Int_t index0=i0*4+i1;
+      Double_t xy[1];
+      Double_t xfi[1];
+      xy[0] = pyf[1];   
+      xfi[0] = pyf[1];
+      //
+      Double_t dy = params[i1][0]-params[i0][0];
+      Double_t sy = TMath::Sqrt(errors[i1][0]*errors[i1][0]+errors[i0][0]*errors[i0][0]);
+      Double_t dphi = params[i1][1]-params[i0][1];
+      Double_t sphi = TMath::Sqrt(errors[i1][1]*errors[i1][1]+errors[i0][1]*errors[i0][1]);
+      //
+      Int_t sector = isec-36;
+      //
+      if (TMath::Abs(dy/(sy+0.1))>5.) continue;          // 5 sigma cut
+      if (TMath::Abs(dphi/(sphi+0.004))>5.) continue;    // 5 sigma cut 
+      TLinearFitter * fitterY,*fitterPhi;
+      fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*sector+index0);
+      if (TMath::Abs(dy/(sy+0.1))<5.)  fitterY->AddPoint(xy,dy,sy);
+      fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*36+index0);
+      if (TMath::Abs(dy/(sy+0.1))<5.)  fitterY->AddPoint(xy,dy,sy);
+      //
+      fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*sector+index0);
+      if (TMath::Abs(dphi/(sphi+0.001))<5.)  fitterPhi->AddPoint(xfi,dphi,sphi);
+      fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*36+index0);
+      if (TMath::Abs(dphi/(sphi+0.001))<5.)  fitterPhi->AddPoint(xfi,dphi,sphi);
+    }
+  }
+  //
+  // Dump debug information
+  //
+  if (fStreamLevel>0){    
+    TTreeSRedirector *cstream = GetDebugStreamer();      
+    Int_t sector = isec-36;
+    if (cstream){      
+      for (Int_t i0=0;i0<4;i0++){
+       for (Int_t i1=i0;i1<4;i1++){
+         if (i0==i1) continue;
+         if(npoints[i0]<kMinClusterQ) continue;
+         if(npoints[i1]<kMinClusterQ) continue;
+         (*cstream)<<"fitQ"<<
+           "run="<<fRun<<              //  run number
+           "event="<<fEvent<<          //  event number
+           "time="<<fTime<<            //  time stamp of event
+           "trigger="<<fTrigger<<      //  trigger
+           "mag="<<fMagF<<             //  magnetic field        
+           "sec="<<sector<<            // current sector from 0 to 36
+           "isec="<<isec<<             //  current sector
+           // full fit
+           "nf="<<nf<<                 //  total number of points
+           "pyf.="<<&pyf<<             //  full OROC fit y
+           "pzf.="<<&pzf<<             //  full OROC fit z
+           // quadrant fit
+           "i0="<<i0<<                 // quadrant number
+           "i1="<<i1<<                 
+           "n0="<<npoints[i0]<<        // number of points
+           "n1="<<npoints[i1]<<
+           "p0.="<<&params[i0]<<       // parameters
+           "p1.="<<&params[i1]<< 
+           "e0.="<<&errors[i0]<<       // errors
+           "e1.="<<&errors[i1]<<
+           "chi0="<<chi2C[i0]<<       // chi2s
+           "chi1="<<chi2C[i1]<<
+
+           "\n";
+       }    
+      }
+    }
+  }
+}
+
+
+
+
+
 void AliTPCcalibUnlinearity::MakeHisto(){
   //
   //
@@ -282,38 +554,53 @@ void AliTPCcalibUnlinearity::MakeHisto(){
   //
   axisName[0]="sector";
   xmin[0]=0; xmax[0]=72; nbins[0]=72;
-  //
+  //          
   axisName[1]="row";
-  xmin[1]=0; xmax[1]=159; nbins[1]=159;
-  //
-  axisName[2]="ly";
-  xmin[2]=-50; xmax[2]=50; nbins[2]=10;
-  //
-  axisName[3]="lz";
-  xmin[3]=-250; xmax[3]=250; nbins[3]=50;
-  //
-  axisName[4]="p2";
-  xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
-  //
-  axisName[5]="p3";
-  xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
+  xmin[1]=0; xmax[1]=96; nbins[1]=96;
   //
-  axisName[6]="p4";
-  xmin[6]=-2; xmax[6]=2; nbins[6]=20;
-  //
-  axisName[7]="dy";
-  xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
+  axisName[2]="tphi";            // tan phi - ly/lx 
+  xmin[2]=-0.17; xmax[2]=0.17; nbins[2]=5;
+  //  
+  axisName[3]="lz";              // global z        
+  xmin[3]=-250; xmax[3]=250; nbins[3]=10;
+  //
+  axisName[4]="k";              // tangent - in angle
+  xmin[4]=-0.5; xmax[4]=0.5; nbins[4]=10;
+  //
+  //
+  axisName[5]="delta";          // delta
+  xmin[5]=-0.5; xmax[5]=0.5; nbins[5]=100;
+  //
+  //
+  fDiffHistoLineY = new THnSparseF("hnDistHistoLineY","hnDistHistoLineY",6,nbins,xmin,xmax);
+  fDiffHistoParY = new THnSparseF("hnDistHistoParY","hnDistHistoParY",6,nbins,xmin,xmax);
+  fDiffHistoLineZ = new THnSparseF("hnDistHistoLineZ","hnDistHistoLineZ",6,nbins,xmin,xmax);
+  fDiffHistoParZ = new THnSparseF("hnDistHistoParZ","hnDistHistoParZ",6,nbins,xmin,xmax);
+ //  for (Int_t i=0; i<8;i++){
+//     fDiffHistoLineY->GetAxis(i)->SetName(axisName[i].Data());
+//     fDiffHistoParY->GetAxis(i)->SetName(axisName[i].Data());
+//     fDiffHistoLineY->GetAxis(i)->SetTitle(axisName[i].Data());
+//     fDiffHistoParY->GetAxis(i)->SetTitle(axisName[i].Data());
+//     fDiffHistoLineZ->GetAxis(i)->SetName(axisName[i].Data());
+//     fDiffHistoParZ->GetAxis(i)->SetName(axisName[i].Data());
+//     fDiffHistoLineZ->GetAxis(i)->SetTitle(axisName[i].Data());
+//     fDiffHistoParZ->GetAxis(i)->SetTitle(axisName[i].Data());
+//   }
+  
   //
-  axisName[8]="dz";
-  xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
+  // 
   //
-  fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
-  fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
-  for (Int_t i=0; i<9;i++){
-    fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
-    fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
-    fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
-    fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
+  char hname[1000];
+  TH2F * his=0;
+  for (Int_t isec=0;isec<74;isec++){
+    sprintf(hname,"DeltaRPhiPlus_Sector%d",isec);
+    his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
+    his->SetDirectory(0);
+    fDistRPHIPlus.AddAt(his,isec);
+    sprintf(hname,"DeltaRPhiMinus_Sector%d",isec);
+    his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
+    his->SetDirectory(0);
+    fDistRPHIMinus.AddAt(his,isec);
   }
 }
 
@@ -329,13 +616,77 @@ void AliTPCcalibUnlinearity::Terminate(){
   AliTPCcalibBase::Terminate();
 }
 
+Long64_t AliTPCcalibUnlinearity::Merge(TCollection *list) {
+  //
+  // merge results
+  //
+  TIterator* iter = list->MakeIterator();
+  AliTPCcalibUnlinearity* cal = 0;
+  //
+  while ((cal = (AliTPCcalibUnlinearity*)iter->Next())) {
+    if (!cal->InheritsFrom(AliTPCcalibUnlinearity::Class())) {
+      return -1;
+    }    
+    Add(cal);
+  }
+  return 0;  
+}
+
+void    AliTPCcalibUnlinearity::Add(AliTPCcalibUnlinearity * calib){
+  //
+  //
+  //
+  if (!fInit) Init();
+  if (fDiffHistoLineY && calib->fDiffHistoLineY){
+    fDiffHistoLineY->Add(calib->fDiffHistoLineY);
+  }
+  if (fDiffHistoParY && calib->fDiffHistoParY){
+    fDiffHistoParY->Add(calib->fDiffHistoParY);
+  }
+  if (fDiffHistoLineZ && calib->fDiffHistoLineZ){
+    fDiffHistoLineZ->Add(calib->fDiffHistoLineZ);
+  }
+  if (fDiffHistoParZ && calib->fDiffHistoParZ){
+    fDiffHistoParZ->Add(calib->fDiffHistoParZ);
+  }
+  for (Int_t isec=0;isec<38;isec++){
+    TLinearFitter *f0r = (TLinearFitter*)fFittersOutR.At(isec);
+    TLinearFitter *f1r = (TLinearFitter*)(calib->fFittersOutR.At(isec));
+    if (f0r&&f1r) f0r->Add(f1r);
+    TLinearFitter *f0z = (TLinearFitter*)fFittersOutZ.At(isec);
+    TLinearFitter *f1z = (TLinearFitter*)(calib->fFittersOutZ.At(isec));
+    if (f0z&&f1z) f0z->Add(f1z);
+  }
+
+  for (Int_t isec=0;isec<16*38;isec++){
+    TLinearFitter *f0y = (TLinearFitter*)fFitterQuadrantY.At(isec);
+    TLinearFitter *f1y = (TLinearFitter*)(calib->fFitterQuadrantY.At(isec));
+    if (f0y&&f1y) f0y->Add(f1y);
+    TLinearFitter *f0phi = (TLinearFitter*)fFitterQuadrantPhi.At(isec);
+    TLinearFitter *f1phi = (TLinearFitter*)(calib->fFitterQuadrantPhi.At(isec));
+    if (f0phi&&f1phi) f0phi->Add(f1phi);
+  }
+  
+  
+  for (Int_t isec=0;isec<74;isec++){
+    TH2F * h0p= (TH2F*)(fDistRPHIPlus.At(isec));
+    TH2F * h1p= (TH2F*)(calib->fDistRPHIPlus.At(isec));
+    if (h0p&&h1p) h0p->Add(h1p);
+    TH2F * h0m= (TH2F*)(fDistRPHIMinus.At(isec));
+    TH2F * h1m= (TH2F*)(calib->fDistRPHIMinus.At(isec));
+    if (h0m&&h1m) h0m->Add(h1m);
+  }
+
+}
+
+
 
-void AliTPCcalibUnlinearity::DumpTree(){
+void AliTPCcalibUnlinearity::DumpTree(const char *fname){
   //
   //
   // 
-  if (fStreamLevel==0) return;
-  TTreeSRedirector *cstream = GetDebugStreamer();
+  TTreeSRedirector *cstream =new TTreeSRedirector(fname);
   if (!cstream) return;
   //  
   THnSparse *his=0;  
@@ -344,7 +695,7 @@ void AliTPCcalibUnlinearity::DumpTree(){
   Int_t *bins = new Int_t[10];
   //
   for (Int_t ihis=0; ihis<2; ihis++){
-    his =  (ihis==0)? fDiffHistoLine:fDiffHistoPar;
+    his =  (ihis==0)? fDiffHistoLineY:fDiffHistoParY;
     if (!his) continue;
     //
     Int_t ndim = his->GetNdimensions();
@@ -364,10 +715,11 @@ void AliTPCcalibUnlinearity::DumpTree(){
        "\n";      
     }
   }
+  delete cstream;
 }
 
 
-void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints){
+void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t ky, Double_t kz, Int_t npoints){
   //
   //
   // Simple distortion fit in outer sectors
@@ -376,24 +728,17 @@ void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Doub
   //       - Decay length 10 and 5 cm
   //       - Decay amplitude - Z dependent 
   //
-  /*
-  chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
-  chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
-  chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
-  chainUnlin->SetAlias("dzt","dz*side");
-  chainUnlin->SetAlias("dr","(1-abs(lz/250))");
-  chainUnlin->SetAlias("dout","(159-row)*1.5");
-  chainUnlin->SetAlias("din","row*0.75");
-  chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
-  chainUnlin->SetAlias("eout5","exp(-(dout)/5.)");  
-  */
-
   Double_t side   = (-1+((sector%36)<18)*2); // A side +1   C side -1
-  Double_t dzt    = dz*side;
-  Double_t dout   = (159-row)*1.5;  // distance to the last pad row - (valid only for long pads)
-  if (dout>30) return; // use only edge pad rows
-  dr-=0.5;             // dr shifted to the middle - reduce numerical instabilities
-
+  Double_t dzt    = (cz-tz)*side;
+  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;
+  if (dout>30) return;
+  //drift
+  Double_t dr   = 0.5 - TMath::Abs(cz/250.);
+  Double_t dy = cy-ty;
   Double_t eout10 = TMath::Exp(-dout/10.);
   Double_t eout5 = TMath::Exp(-dout/5.);
   //
@@ -410,19 +755,19 @@ void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Doub
   xxxZ[4]=dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
   xxxZ[5]=dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
   //
-  xxxR[0]=p2*eoutp;                //fstring+="eoutp++";  
-  xxxR[1]=p2*eoutm;               //fstring+="eoutm++";  
-  xxxR[2]=p2*dr*eoutp;            //fstring+="dr*eoutp++";  
-  xxxR[3]=p2*dr*eoutm;            //fstring+="dr*eoutm++";  
-  xxxR[4]=p2*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
-  xxxR[5]=p2*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
+  xxxR[0]=ky*eoutp;                //fstring+="eoutp++";  
+  xxxR[1]=ky*eoutm;               //fstring+="eoutm++";  
+  xxxR[2]=ky*dr*eoutp;            //fstring+="dr*eoutp++";  
+  xxxR[3]=ky*dr*eoutm;            //fstring+="dr*eoutm++";  
+  xxxR[4]=ky*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
+  xxxR[5]=ky*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
   //
-  xxxRZ[0]=p3*eoutp;                //fstring+="eoutp++";  
-  xxxRZ[1]=p3*eoutm;               //fstring+="eoutm++";  
-  xxxRZ[2]=p3*dr*eoutp;            //fstring+="dr*eoutp++";  
-  xxxRZ[3]=p3*dr*eoutm;            //fstring+="dr*eoutm++";  
-  xxxRZ[4]=p3*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
-  xxxRZ[5]=p3*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
+  xxxRZ[0]=kz*eoutp;                //fstring+="eoutp++";  
+  xxxRZ[1]=kz*eoutm;               //fstring+="eoutm++";  
+  xxxRZ[2]=kz*dr*eoutp;            //fstring+="dr*eoutp++";  
+  xxxRZ[3]=kz*dr*eoutm;            //fstring+="dr*eoutm++";  
+  xxxRZ[4]=kz*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
+  xxxRZ[5]=kz*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
   //
   TLinearFitter * fitter=0;
   Double_t err=0.1;
@@ -450,12 +795,36 @@ void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Doub
     }
   }
 }
+void AliTPCcalibUnlinearity::AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t /*ky*/, Double_t /*kz*/, Int_t npoints){
+  //
+  //
+  //
+  Float_t kMaxDz = 0.5;  // cut on z in cm
+  Double_t dy = cy-ty;
+  Double_t dz = cz-tz;
+  if (TMath::Abs(dz)>kMaxDz) return;
+  Double_t edgeMinus= -ty+cx*TMath::Tan(TMath::Pi()/18.);
+  Double_t edgePlus =  ty+cx*TMath::Tan(TMath::Pi()/18.);
+  //
+  TH2F * his =0;
+  his = (TH2F*)fDistRPHIPlus.At(sector);
+  his->Fill(edgePlus,dy,npoints);
+  his = (TH2F*)((sector<36)? fDistRPHIPlus.At(72):fDistRPHIPlus.At(73));
+  his->Fill(edgePlus,dy,npoints);
+  his = (TH2F*)fDistRPHIMinus.At(sector);
+  his->Fill(edgeMinus,dy,npoints);
+  his = (TH2F*)((sector<36)? fDistRPHIMinus.At(72):fDistRPHIMinus.At(73));
+  his->Fill(edgeMinus,dy,npoints);
+}
+
 
 
-void AliTPCcalibUnlinearity::MakeFitters(){
+void AliTPCcalibUnlinearity::Init(){
   //
   //
   //
+  //
+  MakeHisto();
   // Make outer fitters
   TLinearFitter * fitter = 0;
   for (Int_t ifit=0; ifit<38; ifit++){
@@ -466,6 +835,15 @@ void AliTPCcalibUnlinearity::MakeFitters(){
     fitter->StoreData(kFALSE);
     fFittersOutZ.AddAt(fitter,ifit);
   }
+  for (Int_t ifit=0; ifit<16*38;ifit++){
+    fitter = new TLinearFitter(2,"hyp1");
+    fitter->StoreData(kFALSE);
+    fFitterQuadrantY.AddAt(fitter,ifit);
+    fitter = new TLinearFitter(2,"hyp1");
+    fitter->StoreData(kFALSE);
+    fFitterQuadrantPhi.AddAt(fitter,ifit);    
+  }  
+  fInit= kTRUE;
 }
 
 void AliTPCcalibUnlinearity::EvalFitters(){
@@ -495,97 +873,6 @@ void AliTPCcalibUnlinearity::EvalFitters(){
   }
 }
 
-Double_t AliTPCcalibUnlinearity::GetDr(Int_t sector, Double_t dout, Double_t dr){
-  //
-  //
-  //
-  TVectorD * vec = GetParamOutR(sector);
-  if (!vec) return 0;
-  dr-=0.5;        // dr=0 at the middle drift length
-  Double_t eout10 = TMath::Exp(-dout/10.);
-  Double_t eout5 = TMath::Exp(-dout/5.);                   
-  Double_t eoutp  = (eout10+eout5)*0.5;
-  Double_t eoutm  = (eout10-eout5)*0.5;
-
-  Double_t result=0;
-  result+=(*vec)[1]*eoutp;
-  result+=(*vec)[2]*eoutm;
-  result+=(*vec)[3]*eoutp*dr;
-  result+=(*vec)[4]*eoutm*dr;
-  result+=(*vec)[5]*eoutp*dr*dr;
-  result+=(*vec)[6]*eoutm*dr*dr;
-  return result;
-}
-
-
-Double_t AliTPCcalibUnlinearity::GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
-  //
-  //
-  //
-  Double_t alpha    = TMath::ATan2(gy,gx);
-  if (alpha<0)  alpha+=TMath::Pi()*2;
-  Int_t lsector = Int_t(18*alpha/(TMath::Pi()*2));
-  alpha = (lsector+0.5)*TMath::Pi()/9.;
-  //
-  Double_t r[3];
-  r[0]=gx; r[1]=gy;
-  Double_t cs=TMath::Cos(-alpha), sn=TMath::Sin(-alpha), x=r[0];
-  r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
-  //
-  AliTPCROC * roc = AliTPCROC::Instance();
-  Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
-  Double_t dout = xout-r[0];
-  if (dout<0) return 0;
-  Double_t dr   = 1-TMath::Abs(gz/250.);
-  if (gz<0) lsector+=18;
-  if (stype==0) lsector = (gz>0) ? 36:37;
-  if (stype<0) return lsector;  // 
-  return GetDr(lsector,dout,dr);
-}
-
-
-Double_t AliTPCcalibUnlinearity::GetDz(Int_t sector, Double_t dout, Double_t dr){
-  //
-  //
-  //
-  TVectorD * vec = GetParamOutZ(sector);
-  if (!vec) return 0;
-  dr-=0.5; // 0 at the middle
-  Double_t eout10 = TMath::Exp(-dout/10.);
-  Double_t eout5 = TMath::Exp(-dout/5.);
-  Double_t eoutp  = (eout10+eout5)*0.5;
-  Double_t eoutm  = (eout10-eout5)*0.5;
-
-  
-  Double_t result=0;
-  result+=(*vec)[1]*eoutp;
-  result+=(*vec)[2]*eoutm;
-  result+=(*vec)[3]*eoutp*dr;
-  result+=(*vec)[4]*eoutm*dr;
-  result+=(*vec)[5]*eoutp*dr*dr;
-  result+=(*vec)[6]*eoutm*dr*dr;
-  return result;
-}
-
-
-Double_t      AliTPCcalibUnlinearity::SGetDr(Int_t sector, Double_t dout, Double_t dr){
-  //
-  //
-  //
-  return fgInstance->GetDr(sector,dout,dr);
-}
-Double_t      AliTPCcalibUnlinearity::SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
-  //
-  //
-  //
-  return fgInstance->GetGDr(stype,gx,gy,gz);
-}
-Double_t      AliTPCcalibUnlinearity::SGetDz(Int_t sector, Double_t dout, Double_t dr){
-  //
-  //
-  //
-  return fgInstance->GetDz(sector,dout,dr);
-}
 
 
 
@@ -619,10 +906,12 @@ void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
     if (i%10000==0) printf("%d\n",(Int_t)i);
     Int_t row= cl->GetRow();
     if (cl->GetDetector()>36) row+=64;
+    calib->AddPointRPHI(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
+                     ty,tz,(*vy)[1],(*vz)[1],1);
+
     if (cl->GetType()<0) continue; 
     Double_t dy = cl->GetY()-ty;
     Double_t dz = cl->GetZ()-tz;
-    Double_t dr = 1-TMath::Abs(cl->GetZ()/250);
     //
     //
     if (TMath::Abs(dy)>0.25) continue;
@@ -631,8 +920,10 @@ void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
     if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
     if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
       // Apply sagita cut
-    calib->AddPoint(cl->GetDetector(), row, dz, dy,
-                   (*vy)[1],(*vz)[1], dr, 1);
+    if (cl->GetType()>=0) {
+      calib->AddPoint(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
+                     ty,tz,(*vy)[1],(*vz)[1],1);
+    }
   }
   }
 }
index d447c31..382eccf 100644 (file)
@@ -20,6 +20,7 @@ class TList;
 class AliESDEvent;
 class AliESDtrack;
 class TLinearFitter;
+class AliTPCClusterParam;
 
  
 class AliTPCcalibUnlinearity:public AliTPCcalibBase {
@@ -31,38 +32,38 @@ public:
   virtual void Process(AliTPCseed *track);
   virtual void Analyze(){return;}
   virtual void Terminate();
-  virtual Long64_t Merge(TCollection* list){return 0;}
+  virtual Long64_t Merge(TCollection* list);
+  void    Add(AliTPCcalibUnlinearity * calib);
   //
   void ProcessTree(TTree * tree, Long64_t nmaxPoints);
-  void AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints=1);
+  void AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t ky, Double_t kz, Int_t npoints=1);
+  void AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t ky, Double_t kz, Int_t npoints=1);
   //
   void MakeHisto();
   void ProcessDiff(AliTPCseed *track, Int_t isec);
-  void DumpTree();
-  void MakeFitters();
+  void AlignOROC(AliTPCseed *track, Int_t isec);
+  //
+  void DumpTree(const char *fname="unlinResidual.root");
+  void Init();
   void EvalFitters();
   TLinearFitter * GetFitterOutR(Int_t sector) {return (TLinearFitter*)fFittersOutR.At(sector);}
   TLinearFitter * GetFitterOutZ(Int_t sector) {return (TLinearFitter*)fFittersOutZ.At(sector);}
   TVectorD * GetParamOutR(Int_t sector) {return (TVectorD*)fParamsOutR.At(sector);}
   TVectorD * GetParamOutZ(Int_t sector) {return (TVectorD*)fParamsOutZ.At(sector);}
   //
-  Double_t      GetDr(Int_t sector, Double_t dout, Double_t dr);
-  Double_t      GetDz(Int_t sector, Double_t dout, Double_t dr);
-  Double_t      GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz);
-  //
-  static Double_t      SGetDr(Int_t sector, Double_t dout, Double_t dr);
-  static Double_t      SGetDz(Int_t sector, Double_t dout, Double_t dr);
-  static Double_t      SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz);
-
-  static AliTPCcalibUnlinearity* Instance();
-  void SetInstance(AliTPCcalibUnlinearity*param){fgInstance = param;}
   //TMatrixD * GetNormCovariance(Int_t sector, Int_t type);
   //TMatrixD * GetNormCovariance(Int_t sector, Int_t type);
   static void MakeQPosNormAll(TTree * chainres, AliTPCClusterParam * param, Int_t maxPoints);
+  void     Process(AliESDEvent *event) {AliTPCcalibBase::Process(event);};
+  void     Process(AliESDtrack *track, Int_t runNo=-1){AliTPCcalibBase::Process(track,runNo);};
 public:
-  THnSparse * fDiffHistoLine;    // matrix with cluster residuals - linear fit
-  THnSparse * fDiffHistoPar;     // matrix with cluster residuals - parabolic fit
+  Bool_t      fInit;              // initialization flag
+  THnSparse * fDiffHistoLineY;    // matrix with cluster residuals - linear fit
+  THnSparse * fDiffHistoParY;     // matrix with cluster residuals - parabolic fit
+  THnSparse * fDiffHistoLineZ;    // matrix with cluster residuals - linear fit
+  THnSparse * fDiffHistoParZ;     // matrix with cluster residuals - parabolic fit
   //
+  // Outer residula fit
   //
   TObjArray   fFittersOutR;      // Unlinearity fitters for radial distortion  - outer field cage
   TObjArray   fFittersOutZ;      // Unlinearity fitters for z      distortion  - outer field cage
@@ -71,13 +72,20 @@ public:
   TObjArray   fErrorsOutR;       // Parameters  for radial distortion  - outer field cage
   TObjArray   fErrorsOutZ;      // Parameters  for z      distortion  - outer field cage
   //
-
+  // R-phi residual histogram
+  //
+  TObjArray   fDistRPHIPlus;     // edge effect histograms  - plus  direction 
+  TObjArray   fDistRPHIMinus;    // edge effect histograms  - minus direction
+  //
+  // Quadrant fitters
+  //
+   TObjArray   fFitterQuadrantY;        //qudrant misalignemnt fit Y
+   TObjArray   fFitterQuadrantPhi;      //qudrant misalignemnt fit Phi
 private:
   AliTPCcalibUnlinearity(const AliTPCcalibUnlinearity&); 
   AliTPCcalibUnlinearity& operator=(const AliTPCcalibUnlinearity&); 
- static AliTPCcalibUnlinearity*   fgInstance; //! Instance of this class (singleton implementation)
  
-  ClassDef(AliTPCcalibUnlinearity, 1); 
+  ClassDef(AliTPCcalibUnlinearity, 4); 
 };
 
 #endif