2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
20 Basic calibration and QA class for the TPC gain calibration based on tracks from BEAM events.
23 Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch
27 #include "Riostream.h"
40 #include "TGraphErrors.h"
42 #include "AliTPCcalibDB.h"
43 #include "AliTPCclusterMI.h"
44 #include "AliTPCClusterParam.h"
45 #include "AliTPCseed.h"
46 #include "AliESDVertex.h"
47 #include "AliESDEvent.h"
48 #include "AliESDfriend.h"
49 #include "AliESDInputHandler.h"
50 #include "AliAnalysisManager.h"
51 #include "AliTPCParam.h"
53 #include "AliComplexCluster.h"
54 #include "AliTPCclusterMI.h"
58 #include "AliTPCcalibGainMult.h"
60 #include "TTreeStream.h"
61 #include "TDatabasePDG.h"
62 #include "AliKFParticle.h"
63 #include "AliKFVertex.h"
65 #include "AliESDkink.h"
66 #include "AliRecoParam.h"
67 #include "AliTracker.h"
68 #include "AliTPCTransform.h"
69 #include "AliTPCROC.h"
70 #include "TStatToolkit.h"
72 ClassImp(AliTPCcalibGainMult)
74 Double_t AliTPCcalibGainMult::fgMergeEntriesCut=10000000.;
76 AliTPCcalibGainMult::AliTPCcalibGainMult()
84 fCutRequireITSrefit(0),
105 // Empty default cosntructor
107 AliDebug(5,"Default Constructor");
111 AliTPCcalibGainMult::AliTPCcalibGainMult(const Text_t *name, const Text_t *title)
119 fCutRequireITSrefit(0),
126 fHistClusterShape(0),
146 fLowerTrunc = 0.02; // IMPORTANT CHANGE --> REMOVE HARDWIRED TRUNCATION FROM TRACKER
148 fUseMax = kTRUE; // IMPORTANT CHANGE FOR PbPb; standard: kFALSE;
150 // define track cuts and default BB parameters for interpolation around mean
154 fCutRequireITSrefit = kFALSE;
157 // default values for MIP window selection
158 fMinMomentumMIP = 0.4;
159 fMaxMomentumMIP = 0.6;
160 fAlephParameters[0] = 0.07657; // the following parameters work for most of the periods and are therefore default
161 fAlephParameters[1] = 10.6654; // but they can be overwritten in the train setup of cpass0
162 fAlephParameters[2] = 2.51466e-14;
163 fAlephParameters[3] = 2.05379;
164 fAlephParameters[4] = 1.84288;
166 // basic QA histograms - mainly for debugging purposes
168 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5);
169 fHistClusterShape = new TH1F("fHistClusterShape","cluster shape; rms meas. / rms exp.;",300,0,3);
170 fHistQA = new TH3F("fHistQA","dEdx; momentum p (GeV); TPC signal (a.u.); pad",500,0.1,20.,500,0.,500,6,-0.5,5.5);
171 AliTPCcalibBase::BinLogX(fHistQA);
174 // MIP, sect, pad (short,med,long,full,oroc), run, ncl
175 Int_t binsGainSec[5] = { 145, 72, 4, 10000000, 65};
176 Double_t xminGainSec[5] = { 10., -0.5, -0.5, -0.5, -0.5};
177 Double_t xmaxGainSec[5] = {300., 71.5, 3.5, 9999999.5, 64.5};
178 TString axisNameSec[5]={"Q","sector","pad type","run", "ncl"};
179 TString axisTitleSec[5]={"Q (a.u)","sector","pad type","run","ncl"};
181 fHistGainSector = new THnSparseF("fHistGainSector","0:MIP, 1:sect, 2:pad, 3:run, 4:ncl", 5, binsGainSec, xminGainSec, xmaxGainSec);
182 for (Int_t iaxis=0; iaxis<5;iaxis++){
183 fHistGainSector->GetAxis(iaxis)->SetName(axisNameSec[iaxis]);
184 fHistGainSector->GetAxis(iaxis)->SetTitle(axisTitleSec[iaxis]);
187 // pad region equalization
189 Int_t binsPadEqual[5] = { 400, 400, 4, 5, 20};
190 Double_t xminPadEqual[5] = { 0.0, 0.0, -0.5, 0, -1.};
191 Double_t xmaxPadEqual[5] = { 2.0, 2.0, 3.5, 13000, +1};
192 TString axisNamePadEqual[5] = {"dEdxRatioMax","dEdxRatioTot","padType","mult","dipAngle"};
193 TString axisTitlePadEqual[5] = {"dEdx_padRegion/mean_dEdx Qmax", "dEdx_padRegion/mean_dEdx Qtot","padType","mult","tan(lambda)"};
195 fHistPadEqual = new THnSparseF("fHistPadEqual","0:dEdx_pad/dEdx_mean, 1:pad, 2:mult, 3:drift", 5, binsPadEqual, xminPadEqual, xmaxPadEqual);
196 for (Int_t iaxis=0; iaxis<5;iaxis++){
197 fHistPadEqual->GetAxis(iaxis)->SetName(axisNamePadEqual[iaxis]);
198 fHistPadEqual->GetAxis(iaxis)->SetTitle(axisTitlePadEqual[iaxis]);
201 // multiplicity dependence
202 // MIP Qmax, MIP Qtot, z, pad, vtx. contribut., ncl
203 Int_t binsGainMult[6] = { 145, 145, 25, 4, 100, 80};
204 Double_t xminGainMult[6] = { 10., 10., 0, -0.5, 0, -0.5};
205 Double_t xmaxGainMult[6] = {300., 300., 250, 3.5, 13000, 159.5};
206 TString axisNameMult[6]={"Qmax","Qtot","drift","padtype""multiplicity","ncl"};
207 TString axisTitleMult[6]={"Qmax (a.u)","Qtot (a.u.)","driftlenght l (cm)","Pad Type","multiplicity","ncl"};
209 fHistGainMult = new THnSparseF("fHistGainMult","MIP Qmax, MIP Qtot, z, type, vtx. contribut.", 6, binsGainMult, xminGainMult, xmaxGainMult);
210 for (Int_t iaxis=0; iaxis<6;iaxis++){
211 fHistGainMult->GetAxis(iaxis)->SetName(axisNameMult[iaxis]);
212 fHistGainMult->GetAxis(iaxis)->SetTitle(axisTitleMult[iaxis]);
215 // dip-angle (theta or eta) and inclination angle (local phi) dependence -- absolute calibration
217 // (0.) weighted dE/dx, (1.) 0:Qtot - 1:Qmax, (2.) tgl, (3.) 1./pT
218 Int_t binsTopology[4] = {145, 2, 20, 20};
219 Double_t xminTopology[4] = { 10, -0.5, -1, 0};
220 Double_t xmaxTopology[4] = { 300, 1.5, +1, 5};
221 TString axisNameTopology[4] = {"dEdx", "QtotQmax", "tgl", "1/pT"};
222 TString axisTitleTopology[4] = {"dEdx", "QtotQmax", "tgl", "1/pT"};
224 fHistTopology = new THnF("fHistTopology", "dEdx,QtotQmax,tgl, 1/pT", 4, binsTopology, xminTopology, xmaxTopology);
225 for (Int_t iaxis=0; iaxis<4;iaxis++){
226 fHistTopology->GetAxis(iaxis)->SetName(axisNameTopology[iaxis]);
227 fHistTopology->GetAxis(iaxis)->SetTitle(axisTitleTopology[iaxis]);
230 // MI suggestion for all dEdx histograms we shpuld keep log scale - to have the same relative resolution
232 // e.g: I want to enable - AliTPCcalibBase::BinLogX(fHistTopology->GetAxis(0));
236 // dedx maps - bigger granulatity in phi -
237 // to construct the dedx sector/phi map
238 Int_t binsGainMap[4] = { 100, 90, 10, 6};
239 Double_t xminGainMap[4] = { 0.3, -TMath::Pi(), 0, 0};
240 Double_t xmaxGainMap[4] = { 2, TMath::Pi(), 1, 6};
241 TString axisNameMap[4] = {"Q_Qexp","phi", "1/Qexp","Pad Type"};
242 TString axisTitleMap[4] = {"Q/Q_{exp}","#phi (a.u.)","1/Q_{exp}","Pad Type"};
244 fHistdEdxMap = new THnSparseF("fHistdEdxMap","fHistdEdxMap", 4, binsGainMap, xminGainMap, xmaxGainMap);
245 for (Int_t iaxis=0; iaxis<4;iaxis++){
246 fHistdEdxMap->GetAxis(iaxis)->SetName(axisNameMap[iaxis]);
247 fHistdEdxMap->GetAxis(iaxis)->SetTitle(axisTitleMap[iaxis]);
253 Int_t binsGainMax[6] = { 100, 10, 10, 10, 5, 3};
254 Double_t xminGainMax[6] = { 0.5, 0, 0, 0, 0, 0};
255 Double_t xmaxGainMax[6] = { 1.5, 1, 1.0, 1.0, 3000, 3};
256 TString axisNameMax[6] = {"Q_Qexp","1/Qexp", "phi","theta","mult", "Pad Type"};
257 TString axisTitleMax[6] = {"Q/Q_{exp}","1/Qexp", "#phi","#theta","mult","Pad Type"};
259 fHistdEdxMax = new THnSparseF("fHistdEdxMax","fHistdEdxMax", 6, binsGainMax, xminGainMax, xmaxGainMax);
260 fHistdEdxTot = new THnSparseF("fHistdEdxTot","fHistdEdxTot", 6, binsGainMax, xminGainMax, xmaxGainMax);
261 for (Int_t iaxis=0; iaxis<6;iaxis++){
262 fHistdEdxMax->GetAxis(iaxis)->SetName(axisNameMax[iaxis]);
263 fHistdEdxMax->GetAxis(iaxis)->SetTitle(axisTitleMax[iaxis]);
264 fHistdEdxTot->GetAxis(iaxis)->SetName(axisNameMax[iaxis]);
265 fHistdEdxTot->GetAxis(iaxis)->SetTitle(axisTitleMax[iaxis]);
268 AliDebug(5,"Non Default Constructor");
272 AliTPCcalibGainMult::~AliTPCcalibGainMult(){
276 delete fHistNTracks; // histogram showing number of ESD tracks per event
277 delete fHistClusterShape; // histogram to check the cluster shape
278 delete fHistQA; // dE/dx histogram showing the final spectrum
280 delete fHistGainSector; // histogram which shows MIP peak for each of the 3x36 sectors (pad region)
281 delete fHistPadEqual; // histogram for the equalization of the gain in the different pad regions -> pass0
282 delete fHistGainMult; // histogram which shows decrease of MIP signal as a function
283 delete fHistTopology;
289 if (fBBParam) delete fBBParam;
294 void AliTPCcalibGainMult::Process(AliESDEvent *event) {
296 // Main function of the class
297 // 1. Select Identified particles - for identified particles the flag in the PID matrix is stored
298 // 1.a) ProcessV0s - select Electron (gamma coversion) and pion canditates (from K0s)
299 // 1.b) ProcessTOF - select - Proton, kaon and pions candidates
300 // AS THE TOF not calibrated yet in Pass0 - we are calibrating the TOF T0 in this function
301 // 1.c) ProcessCosmic - select cosmic mumn candidates - too few entries - not significant for the calibration
302 // 1.d) ProcessKinks - select Kaon and pion candidates. From our experience (see Kink debug streamer), the angular cut for kink daughter is not sufficient - big contamination - delta rays, hadronic interaction (proton)
303 // - NOT USED for the
305 // 2. Loop over tracks
306 // 2.a DumpTrack() - for identified particles dump the track and dEdx information into the tree (for later fitting)
307 // 3. Actual fitting for the moment macro
310 // Criteria for the track selection
312 // const Int_t kMinNCL=80; // minimal number of cluster - tracks accepted for the dedx calibration
313 //const Double_t kMaxEta=0.8; // maximal eta fo the track to be accepted
314 //const Double_t kMaxDCAR=10; // maximal DCA R of the track
315 //const Double_t kMaxDCAZ=5; // maximal DCA Z of the track
316 // const Double_t kMIPPt=0.525; // MIP pt
319 Printf("ERROR: ESD not available");
322 fCurrentEvent=event ;
323 fMagF = event->GetMagneticField();
324 Int_t ntracks=event->GetNumberOfTracks();
325 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
327 //Printf("ERROR: esdFriend not available");
331 if (!(esdFriend->TestSkipBit())) fPIDMatrix= new TMatrixD(ntracks,5);
332 fHistNTracks->Fill(ntracks);
333 // ProcessCosmic(event); // usually not enogh statistic
335 if (esdFriend->TestSkipBit()) {
339 //ProcessV0s(event); //
340 //ProcessTOF(event); //
341 //ProcessKinks(event); // not relyable
343 UInt_t runNumber = event->GetRunNumber();
344 Int_t nContributors = event->GetNumberOfTracks();
348 for (Int_t i=0;i<ntracks;++i) {
350 AliESDtrack *track = event->GetTrack(i);
351 if (!track) continue;
353 AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam();
354 if (!trackIn) continue;
356 // calculate necessary track parameters
357 Double_t meanP = trackIn->GetP();
358 Int_t ncls = track->GetTPCNcls();
359 Int_t nCrossedRows = track->GetTPCCrossedRows();
361 // correction factor of dE/dx in MIP window
362 Float_t corrFactorMip = AliExternalTrackParam::BetheBlochAleph(meanP/0.13957,
367 fAlephParameters[4]);
368 if (TMath::Abs(corrFactorMip) < 1e-10) continue;
370 if (nCrossedRows < fCutCrossRows) continue;
371 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
372 if (TMath::Abs(trackIn->Eta()) > fCutEtaWindow) continue;
374 UInt_t status = track->GetStatus();
375 if ((status&AliESDtrack::kTPCrefit)==0) continue;
376 if ((status&AliESDtrack::kITSrefit)==0 && fCutRequireITSrefit) continue; // ITS cluster
377 Float_t dca[2], cov[3];
378 track->GetImpactParameters(dca,cov);
379 Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0]);
380 if (TMath::Abs(dca[0]) > fCutMaxDcaXY || TMath::Abs(dca[0]) < 0.0000001) continue; // cut in xy
381 if (((status&AliESDtrack::kITSrefit) == 1 && TMath::Abs(dca[1]) > 3.) || TMath::Abs(dca[1]) > fCutMaxDcaZ ) continue;
384 // fill Alexander QA histogram
386 if (primVtxDCA < 3 && track->GetNcls(0) > 3 && track->GetKinkIndex(0) == 0 && ncls > 100) fHistQA->Fill(meanP, track->GetTPCsignal(), 5);
389 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
390 if (!friendTrack) continue;
391 TObject *calibObject;
392 AliTPCseed *seed = 0;
393 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
394 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
396 //if (seed) DumpTrack(track, friendTrack,seed,i); // MI implementation for the identified particles
398 if (seed) { // seed the container with track parameters and the clusters
400 const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut(); // tack at the outer radius of TPC
401 if (!trackIn) continue;
402 if (!trackOut) continue;
403 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
404 Double_t dipAngleTgl = trackIn->GetTgl();
406 for (Int_t irow =0; irow<160;irow++) {
407 AliTPCTrackerPoint * point = seed->GetTrackPoint(irow);
408 if (point==0) continue;
409 AliTPCclusterMI * cl = seed->GetClusterPointer(irow);
412 Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
413 fHistClusterShape->Fill(rsigmay);
419 Double_t signalShortMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,0,62);
420 Double_t signalMedMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,63,126);
421 Double_t signalLongMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,127,159);
422 Double_t signalMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1);
423 Double_t signalArrayMax[4] = {signalShortMax, signalMedMax, signalLongMax, signalMax};
425 Double_t signalShortTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,0,62);
426 Double_t signalMedTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,63,126);
427 Double_t signalLongTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,127,159);
429 Double_t signalTot = 0;
431 Double_t signalArrayTot[4] = {signalShortTot, signalMedTot, signalLongTot, signalTot};
433 Double_t mipSignalShort = fUseMax ? signalShortMax : signalShortTot;
434 Double_t mipSignalMed = fUseMax ? signalMedMax : signalMedTot;
435 Double_t mipSignalLong = fUseMax ? signalLongMax : signalLongTot;
436 Double_t mipSignalOroc = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax,63,159);
437 Double_t signal = fUseMax ? signalMax : signalTot;
439 fHistQA->Fill(meanP, mipSignalShort, 0);
440 fHistQA->Fill(meanP, mipSignalMed, 1);
441 fHistQA->Fill(meanP, mipSignalLong, 2);
442 fHistQA->Fill(meanP, signal, 3);
443 fHistQA->Fill(meanP, mipSignalOroc, 4);
445 // normalize pad regions to their common mean
447 Float_t meanMax = (63./159)*signalArrayMax[0] + (64./159)*signalArrayMax[1] + (32./159)*signalArrayMax[2];
448 Float_t meanTot = (63./159)*signalArrayTot[0] + (64./159)*signalArrayTot[1] + (32./159)*signalArrayTot[2];
449 //MY SUGGESTION COMMENT NEXT LINE
450 // if (meanMax < 1e-5 || meanTot < 1e-5) continue;
451 //AND INTRODUCE NEW LINE
453 const Double_t kMinAmp=0.001;
454 if (signalArrayMax[0]<=kMinAmp) continue;
455 if (signalArrayMax[1]<=kMinAmp) continue;
456 if (signalArrayMax[2]<=kMinAmp) continue;
457 if (signalArrayTot[0]<=kMinAmp) continue;
458 if (signalArrayTot[1]<=kMinAmp) continue;
459 if (signalArrayTot[2]<=kMinAmp) continue;
461 for(Int_t ipad = 0; ipad < 4; ipad ++) {
462 // "dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength"
463 Double_t vecPadEqual[5] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, static_cast<Double_t>(ipad), static_cast<Double_t>(nContributors), dipAngleTgl};
464 if (fMinMomentumMIP > meanP && meanP < fMaxMomentumMIP) fHistPadEqual->Fill(vecPadEqual);
468 if (meanP < fMaxMomentumMIP && meanP > fMinMomentumMIP) {
469 Double_t vecMult[6] = {seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1)/corrFactorMip,
470 seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1)/corrFactorMip,
473 static_cast<Double_t>(nContributors),
474 static_cast<Double_t>(ncls)};
476 fHistGainMult->Fill(vecMult);
477 vecMult[0]=mipSignalShort/corrFactorMip; vecMult[1]=mipSignalShort/corrFactorMip; vecMult[3]=0;
478 fHistGainMult->Fill(vecMult);
479 vecMult[0]=mipSignalMed/corrFactorMip; vecMult[1]=mipSignalMed/corrFactorMip; vecMult[3]=1;
480 fHistGainMult->Fill(vecMult);
481 vecMult[0]=mipSignalLong/corrFactorMip; vecMult[1]=mipSignalLong/corrFactorMip; vecMult[3]=2;
482 fHistGainMult->Fill(vecMult);
484 // topology histogram (absolute)
485 // (0.) weighted dE/dx, (1.) 0:Qtot - 1:Qmax, (2.) tgl, (3.) 1./pT
486 Double_t vecTopolTot[4] = {meanTot, 0, dipAngleTgl, TMath::Abs(track->GetSigned1Pt())};
487 Double_t vecTopolMax[4] = {meanMax, 1, dipAngleTgl, TMath::Abs(track->GetSigned1Pt())};
488 fHistTopology->Fill(vecTopolTot);
489 fHistTopology->Fill(vecTopolMax);
493 if (fMinMomentumMIP < meanP || meanP > fMaxMomentumMIP) continue; // only MIP pions
494 //if (meanP > 0.5 || meanP < 0.4) continue; // only MIP pions
496 // for each track, we look at the three different pad regions, split it into tracklets, check if the sector does not change and fill the histogram
498 Bool_t isNotSplit[3] = {kTRUE, kTRUE, kTRUE}; // short, medium, long (true if the track is not split between two chambers)
500 Double_t sector[4] = {-1, -1, -1, -1}; // sector number short, medium, long, all
501 Int_t ncl[3] = {0,0,0};
504 for (Int_t irow=0; irow < 159; irow++){
506 if (irow > 62) padRegion = 1;
507 if (irow > 126) padRegion = 2;
509 AliTPCclusterMI* cluster = seed->GetClusterPointer(irow);
510 if (!cluster) continue;
511 if (sector[padRegion] == -1) {
512 sector[padRegion] = cluster->GetDetector();
515 if (sector[padRegion] != -1 && sector[padRegion] != cluster->GetDetector()) isNotSplit[padRegion] = kFALSE;
519 // MIP, sect, pad, run
521 Double_t vecMip[5] = {mipSignalShort/corrFactorMip, mipSignalMed/corrFactorMip, mipSignalLong/corrFactorMip, signal/corrFactorMip, mipSignalOroc/corrFactorMip};
523 for(Int_t ipad = 0; ipad < 3; ipad++) {
524 // AK. - run Number To be removed - not needed
525 Double_t vecGainSec[5] = {vecMip[ipad], sector[ipad], static_cast<Double_t>(ipad), static_cast<Double_t>(runNumber), static_cast<Double_t>(ncl[ipad])};
526 if (isNotSplit[ipad]) fHistGainSector->Fill(vecGainSec);
536 void AliTPCcalibGainMult::MakeLookup(THnSparse * /*hist*/, Char_t * /*outputFile*/) {
538 // Not yet implemented
543 void AliTPCcalibGainMult::Analyze() {
553 Long64_t AliTPCcalibGainMult::Merge(TCollection *li) {
555 // merging of the component
558 const Int_t kMaxEntriesSparse=2000000; // MI- temporary - restrict the THnSparse size
559 TIterator* iter = li->MakeIterator();
560 AliTPCcalibGainMult* cal = 0;
562 while ((cal = (AliTPCcalibGainMult*)iter->Next())) {
563 if (!cal->InheritsFrom(AliTPCcalibGainMult::Class())) {
564 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
568 if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
569 if (cal->GetHistClusterShape()) fHistClusterShape->Add(cal->GetHistClusterShape());
570 if (cal->GetHistQA()) fHistQA->Add(cal->GetHistQA());
571 if (cal->GetHistGainSector() && fHistGainSector )
573 if ((fHistGainSector->GetEntries()+cal->GetHistGainSector()->GetEntries()) < fgMergeEntriesCut)
575 //AliInfo(Form("fHistGainSector has %.0f tracks, going to add %.0f\n",fHistGainSector->GetEntries(),cal->GetHistGainSector()->GetEntries()));
576 fHistGainSector->Add(cal->GetHistGainSector());
580 AliInfo(Form("fHistGainSector full (has %.0f entries, trying to add %.0f., max allowed: %.0f)",
581 fHistGainSector->GetEntries(),cal->GetHistGainSector()->GetEntries(),fgMergeEntriesCut));
584 if (cal->GetHistPadEqual()) fHistPadEqual->Add(cal->GetHistPadEqual());
585 if (cal->GetHistGainMult()) {
586 if (fHistGainMult->GetEntries()<kMaxEntriesSparse) fHistGainMult->Add(cal->GetHistGainMult());
588 if (cal->GetHistTopology()) {
589 fHistTopology->Add(cal->GetHistTopology());
592 if (cal->fHistdEdxMap){
593 if (fHistdEdxMap) fHistdEdxMap->Add(cal->fHistdEdxMap);
595 if (cal->fHistdEdxMax){
596 if (fHistdEdxMax) fHistdEdxMax->Add(cal->fHistdEdxMax);
598 if (cal->fHistdEdxTot){
599 if (fHistdEdxTot) fHistdEdxTot->Add(cal->fHistdEdxTot);
602 // Originally we tireied to write the tree to the same file as other calibration components
603 // We failed in merging => therefore this optio was disabled
605 // if (cal->fdEdxTree && cal->fdEdxTree->GetEntries()>0) {
607 // const Int_t kMax=100000;
608 // Int_t entriesSum = (Int_t)fdEdxTree->GetEntries();
609 // Int_t entriesCurrent = (Int_t)cal->fdEdxTree->GetEntries();
610 // Int_t entriesCp=TMath::Min((Int_t) entriesCurrent*(kMax*entriesSum),entriesCurrent);
611 // // cal->fdEdxTree->SetBranchStatus("track*",0);
612 // // cal->fdEdxTree->SetBranchStatus("vertex*",0);
613 // // cal->fdEdxTree->SetBranchStatus("tpcOut*",0);
614 // // cal->fdEdxTree->SetBranchStatus("vec*",0);
615 // // fdEdxTree->SetBranchStatus("track*",0);
616 // // fdEdxTree->SetBranchStatus("vertex*",0);
617 // // fdEdxTree->SetBranchStatus("tpcOut*",0);
618 // // fdEdxTree->SetBranchStatus("vec*",0);
619 // fdEdxTree->Print();
620 // fdEdxTree->Dump();
621 // fdEdxTree->GetEntry(0);
622 // for (Int_t i=0; i<entriesCurrent; i++){
623 // cal->fdEdxTree->CopyAddresses(fdEdxTree);
624 // cal->fdEdxTree->GetEntry(i);
625 // fdEdxTree->Fill();
627 // TObjArray *brarray = cal->fdEdxTree->GetListOfBranches();
628 // for (Int_t i=0; i<brarray->GetEntries(); i++) {TBranch * br = (TBranch *)brarray->At(i); br->SetAddress(0); }
631 // fdEdxTree = (TTree*)(cal->fdEdxTree->Clone());
632 // TObjArray *brarray = fdEdxTree->GetListOfBranches();
633 // for (Int_t i=0; i<brarray->GetEntries(); i++) {TBranch * br = (TBranch *)brarray->At(i); br->SetAddress(0);}
647 void AliTPCcalibGainMult::UpdateGainMap() {
649 // read in the old gain map and scale it appropriately...
652 gSystem->Load("libANALYSIS");
653 gSystem->Load("libTPCcalib");
655 TFile jj("Run0_999999999_v1_s0.root");
656 AliTPCCalPad * pad = AliCDBEntry->GetObject()->Clone();
657 TFile hh("output.root");
658 AliTPCcalibGainMult * gain = calibTracksGain;
659 TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,1);
662 histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
663 TH1D * meanGainSec = arr->At(1);
664 Double_t gainsIROC[36];
665 Double_t gainsOROC[36];
668 for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) {
669 cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl;
670 gains[isec-1] = meanGainSec->GetBinContent(isec);
672 gainsIROC[isec-1] = meanGainSec->GetBinContent(isec);
674 gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec);
677 Double_t meanIroc = TMath::Mean(36, gainsIROC);
678 Double_t meanOroc = TMath::Mean(36, gainsIROC);
679 for(Int_t i = 0; i < 36; i++) gains[i] /= meanIroc;
680 for(Int_t i = 36; i < 72; i++) gains[i] /= meanOroc;
682 for(Int_t i = 0; i < 72; i++) {
683 AliTPCCalROC * chamber = pad->GetCalROC(i);
684 chamber->Multiply(gains[i]);
685 cout << i << " "<< chamber->GetMean() << endl;
690 AliCDBMetaData *metaData= new AliCDBMetaData();
691 metaData->SetObjectClassName("AliTPCCalPad");
692 metaData->SetResponsible("Alexander Kalweit");
693 metaData->SetBeamPeriod(1);
694 metaData->SetAliRootVersion("04-19-05"); //root version
695 metaData->SetComment("New gain map for 1600V OROC gain increase and equalization. Valid for runs starting after technical stop beginning of September.");
696 AliCDBId id1("TPC/Calib/GainFactorDedx", 131541, AliCDBRunRange::Infinity()); // important: new gain runs here..
697 AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///d/alice05/akalweit/projects/OCDBupdate/HighGain_2010-09-03/OCDB/");
698 gStorage->Put(pad, id1, metaData);
703 void AliTPCcalibGainMult::UpdateClusterParam() {
708 gSystem->Load("libANALYSIS");
709 gSystem->Load("libTPCcalib");
711 TFile ff("OldClsParam.root");
712 AliTPCClusterParam * param = AliCDBEntry->GetObject()->Clone();
714 TFile hh("output.root");
715 AliTPCcalibGainMult * gain = calibGainMult;
716 TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,2);
718 histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
719 histGainSec->Draw("colz");
720 TH1D * fitVal = arr.At(1);
721 fitVal->Draw("same");
722 param->GetQnormCorrMatrix()->Print();
723 param->GetQnormCorrMatrix()(0,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qtot
724 param->GetQnormCorrMatrix()(1,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qtot
725 param->GetQnormCorrMatrix()(2,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qtot
727 param->GetQnormCorrMatrix()(3,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qmax -> scaling assumed
728 param->GetQnormCorrMatrix()(4,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qmax -> scaling assumed
729 param->GetQnormCorrMatrix()(5,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qmax -> scaling assumed
731 TFile jj("newClusterParam.root","RECREATE");
733 param->GetQnormCorrMatrix()->Print();
737 AliCDBMetaData *metaData= new AliCDBMetaData();
738 metaData->SetObjectClassName("AliTPCClusterParam");
739 metaData->SetResponsible("Alexander Kalweit");
740 metaData->SetBeamPeriod(1);
741 metaData->SetAliRootVersion("04-19-04"); //root version
742 metaData->SetComment("1600V OROC / hard thres. / new algorithm");
743 AliCDBId id1("TPC/Calib/ClusterParam", 0, AliCDBRunRange::Infinity()); // important: new gain runs here..
744 AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///lustre/alice/akalweit/baseline/CalibrationEntries/OldThres_NewAlgo_PP");
745 gStorage->Put(param, id1, metaData);
752 void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftrack, AliTPCseed * seed, Int_t index){
754 // dump interesting tracks
755 // 1. track at MIP region
756 // 2. highly ionizing protons
757 // pidType: 0 - unselected
763 const Int_t kMax=10000;
764 const Int_t kMinRows=80;
765 const Double_t kDCAcut=30;
767 // Bethe Bloch paramerization
769 Double_t kp1= 0.0851148;
770 Double_t kp2= 9.25771;
771 Double_t kp3= 2.6558e-05;
772 Double_t kp4= 2.32742;
773 Double_t kp5= 1.83039;
782 static const AliTPCROC *roc = AliTPCROC::Instance();
783 static const TDatabasePDG *pdg = TDatabasePDG::Instance();
785 Int_t nclITS = track->GetNcls(0);
786 Int_t ncl = track->GetTPCncls();
787 Double_t ncl21 = track->GetTPCClusterInfo(3,1);
788 Double_t ncl20 = track->GetTPCClusterInfo(3,0);
791 if (ncl21<kMinRows) return;
792 static Int_t counter=0;
793 static Int_t counterHPT=0;
795 static TH1F *hisBB=new TH1F("hisBB","hisBB",20,0.1,1.00); // bethe bloch histogram =
796 // used to cover more homegenously differnt dEdx regions
797 static Double_t massPi = pdg->GetParticle("pi-")->Mass(); //
798 static Double_t massK = pdg->GetParticle("K-")->Mass();
799 static Double_t massP = pdg->GetParticle("proton")->Mass();
800 static Double_t massE = pdg->GetParticle("e-")->Mass();
801 static Double_t massMuon = pdg->GetParticle("mu-")->Mass();
802 static Double_t radius0= roc->GetPadRowRadiiLow(roc->GetNRows(0)/2);
803 static Double_t radius1= roc->GetPadRowRadiiUp(30);
804 static Double_t radius2= roc->GetPadRowRadiiUp(roc->GetNRows(36)-15);
806 AliESDVertex *vertex= (AliESDVertex *)fCurrentEvent->GetPrimaryVertex();
808 // Estimate current MIP position -
810 static Double_t mipArray[kMax]; // mip array
811 static Int_t counterMIP0=0;
812 static Double_t medianMIP0=100000; // current MIP median position - estimated after some mimnimum number of MIP tracks
814 if (TMath::Abs(track->GetP()-0.5)<0.1&&track->GetTPCsignal()/medianMIP0<1.5){
815 mipArray[counterMIP0%kMax]= track->GetTPCsignal();
817 if (counterMIP0>10) medianMIP0=TMath::Median(counterMIP0%kMax, mipArray);
819 // the PID as defiend from the external sources
821 Int_t isElectron = TMath::Nint((*fPIDMatrix)(index,0));
822 Int_t isMuon = TMath::Nint((*fPIDMatrix)(index,1));
823 Int_t isPion = TMath::Nint((*fPIDMatrix)(index,2));
824 Int_t isKaon = TMath::Nint((*fPIDMatrix)(index,3));
825 Int_t isProton = TMath::Nint((*fPIDMatrix)(index,4));
827 track->GetImpactParameters(dca[0],dca[1]);
829 if ( (isMuon==0 && isElectron==0) && (TMath::Sqrt(dca[0]*dca[0]+dca[1]*dca[1])>kDCAcut) ) return;
830 Double_t normdEdx= track->GetTPCsignal()/(medianMIP0); // TPC signal normalized to the MIP
832 AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam();
833 AliExternalTrackParam * trackOut = (AliExternalTrackParam *)track->GetOuterParam();
834 AliExternalTrackParam * tpcOut = (AliExternalTrackParam *)ftrack->GetTPCOut();
835 if (!trackIn) return;
836 if (!trackOut) return;
838 if (trackIn->GetZ()*trackOut->GetZ()<0) return; // remove crossing tracks
840 // calculate local and global angle
841 Int_t side = (trackIn->GetZ()>0)? 1:-1;
842 Double_t tgl=trackIn->GetTgl();
843 Double_t gangle[3]={0,0,0};
844 Double_t langle[3]={0,0,0};
845 Double_t length[3]={0,0,0};
846 Double_t pxpypz[3]={0,0,0};
848 trackIn->GetXYZAt(radius0,bz,pxpypz); // get the global position at the middle of the IROC
849 gangle[0]=TMath::ATan2(pxpypz[1],pxpypz[0]); // global angle IROC
850 trackIn->GetXYZAt(radius1,bz,pxpypz); // get the global position at the middle of the OROC - medium pads
851 gangle[1]=TMath::ATan2(pxpypz[1],pxpypz[0]); // global angle OROC
852 trackOut->GetXYZAt(radius2,bz,pxpypz); // get the global position at the middle of OROC - long pads
853 gangle[2]=TMath::ATan2(pxpypz[1],pxpypz[0]);
855 trackIn->GetPxPyPzAt(radius0,bz,pxpypz); //get momentum vector
856 langle[0]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[0]; //local angle between padrow and track IROC
857 trackIn->GetPxPyPzAt(radius1,bz,pxpypz);
858 langle[1]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[1];
859 trackOut->GetPxPyPzAt(radius2,bz,pxpypz); // OROC medium
860 langle[2]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[2];
861 for (Int_t i=0; i<3; i++){
862 if (langle[i]>TMath::Pi()) langle[i]-=TMath::TwoPi();
863 if (langle[i]<-TMath::Pi()) langle[i]+=TMath::TwoPi();
864 length[i]=TMath::Sqrt(1+langle[i]*langle[i]+tgl*tgl); // the tracklet length
867 // Select the kaons and Protons which are "isolated" in TPC dedx curve
870 Double_t dedxP = AliExternalTrackParam::BetheBlochAleph(track->GetInnerParam()->GetP()/massP,kp1,kp2,kp3,kp4,kp5);
871 Double_t dedxK = AliExternalTrackParam::BetheBlochAleph(track->GetInnerParam()->GetP()/massK,kp1,kp2,kp3,kp4,kp5);
872 if (dedxP>2 || dedxK>2){
873 if (track->GetP()<1.2 && normdEdx>1.8&&counterMIP0>10){ // not enough from TOF and V0s triggered by high dedx
874 // signing the Proton and kaon - using the "bitmask" bit 1 and 2 is dedicated for V0s and TOF selected
875 if ( TMath::Abs(normdEdx/dedxP-1)<0.3) isProton+=4;
876 if ( TMath::Abs(normdEdx/dedxK-1)<0.3) isKaon+=4;
877 if (normdEdx/dedxK>1.3) isProton+=8;
878 if (normdEdx/dedxP<0.7) isKaon+=8;
885 Bool_t isHighPt = ((TMath::Power(1/track->Pt(),4)*gRandom->Rndm())<0.005); // rnadomly selected HPT tracks
886 // there are selected for the QA of the obtained calibration
887 Bool_t isMIP = TMath::Abs(track->GetInnerParam()->P()-0.4)<0.005&&(counter<kMax); //
888 // REMINDER - it is not exactly MIP - we select the regtion where the Kaon and Electrons are well separated
890 if (isElectron>0) mass = massE;
891 if (isProton>0) mass = massP;
892 if (isKaon>0) mass = massK;
893 if (isMuon>0) mass = massMuon;
894 if (isPion>0) mass = massPi;
895 if (isHighPt) mass = massPi; //assign mass of pions
896 if (isMIP&&track->GetTPCsignal()/medianMIP0<1.5) mass = massPi; //assign mass of pions
899 // calculate expected dEdx
901 Double_t dedxDefPion= 0,dedxDefProton=0, dedxDefKaon=0;
902 Double_t pin=trackIn->GetP();
903 Double_t pout=trackOut->GetP();
904 Double_t p=(pin+pout)*0.5; // momenta as the mean between tpc momenta at inner and outer wall of the TPC
905 if (mass>0) dedxDef = AliExternalTrackParam::BetheBlochAleph(p/mass,kp1,kp2,kp3,kp4,kp5);
906 dedxDefPion = AliExternalTrackParam::BetheBlochAleph(p/massPi,kp1,kp2,kp3,kp4,kp5);
907 dedxDefProton = AliExternalTrackParam::BetheBlochAleph(p/massP,kp1,kp2,kp3,kp4,kp5);
908 dedxDefKaon = AliExternalTrackParam::BetheBlochAleph(p/massK,kp1,kp2,kp3,kp4,kp5);
910 // dEdx Truncated mean vectros with differnt tuncation
911 // 11 truncations array - 0-10 - 0~50% 11=100%
912 // 3 Regions - IROC,OROC0, OROC1
913 // 2 Q - total charge and max charge
914 // Log - Logarithmic mean used
915 // Up/Dwon - Upper half or lower half of truncation used
916 // RMS - rms of the distribction (otherwise truncated mean)
917 // M2 suffix - second moment ( truncated)
918 TVectorF truncUp(11);
919 TVectorF truncDown(11);
920 TVectorF vecAllMax(11);
921 TVectorF vecIROCMax(11);
922 TVectorF vecOROCMax(11);
923 TVectorF vecOROC0Max(11);
924 TVectorF vecOROC1Max(11);
926 TVectorF vecAllTot(11);
927 TVectorF vecIROCTot(11);
928 TVectorF vecOROCTot(11);
929 TVectorF vecOROC0Tot(11);
930 TVectorF vecOROC1Tot(11);
932 TVectorF vecAllTotLog(11);
933 TVectorF vecIROCTotLog(11);
934 TVectorF vecOROCTotLog(11);
935 TVectorF vecOROC0TotLog(11);
936 TVectorF vecOROC1TotLog(11);
938 TVectorF vecAllTotUp(11);
939 TVectorF vecIROCTotUp(11);
940 TVectorF vecOROCTotUp(11);
941 TVectorF vecOROC0TotUp(11);
942 TVectorF vecOROC1TotUp(11);
944 TVectorF vecAllTotDown(11);
945 TVectorF vecIROCTotDown(11);
946 TVectorF vecOROCTotDown(11);
947 TVectorF vecOROC0TotDown(11);
948 TVectorF vecOROC1TotDown(11);
950 TVectorF vecAllTotRMS(11);
951 TVectorF vecIROCTotRMS(11);
952 TVectorF vecOROCTotRMS(11);
953 TVectorF vecOROC0TotRMS(11);
954 TVectorF vecOROC1TotRMS(11);
956 TVectorF vecAllTotM2(11);
957 TVectorF vecIROCTotM2(11);
958 TVectorF vecOROCTotM2(11);
959 TVectorF vecOROC0TotM2(11);
960 TVectorF vecOROC1TotM2(11);
962 TVectorF vecAllTotMS(11);
963 TVectorF vecIROCTotMS(11);
964 TVectorF vecOROCTotMS(11);
965 TVectorF vecOROC0TotMS(11);
966 TVectorF vecOROC1TotMS(11);
968 // Differnt number of clusters definitions - in separate regions of the TPC
969 // 20 - ratio - found/findabel
970 // 21 - number of clusters used for given dEdx calculation
972 // suffix - 3 or 4 - number of padrows before and after given row to define findable row
974 Double_t ncl20All = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3);
975 Double_t ncl20IROC = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3);
976 Double_t ncl20OROC = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3);
977 Double_t ncl20OROC0= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3);
978 Double_t ncl20OROC1= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3);
980 Double_t ncl20All4 = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3,4);
981 Double_t ncl20IROC4 = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3,4);
982 Double_t ncl20OROC4 = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3,4);
983 Double_t ncl20OROC04= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3,4);
984 Double_t ncl20OROC14= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3,4);
986 Double_t ncl20All3 = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3,3);
987 Double_t ncl20IROC3 = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3,3);
988 Double_t ncl20OROC3 = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3,3);
989 Double_t ncl20OROC03= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3,3);
990 Double_t ncl20OROC13= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3,3);
992 Double_t ncl21All = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,2);
993 Double_t ncl21IROC = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,2);
994 Double_t ncl21OROC = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,2);
995 Double_t ncl21OROC0= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,2);
996 Double_t ncl21OROC1= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,2);
997 // calculate truncated dEdx - mean rms M2 ...
999 for (Int_t ifracDown=0; ifracDown<1; ifracDown++){
1000 for (Int_t ifracUp=0; ifracUp<11; ifracUp++){
1001 Double_t fracDown = 0.0+Double_t(ifracDown)*0.05;
1002 Double_t fracUp = 0.5+Double_t(ifracUp)*0.05;
1003 vecAllMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,0,159,0);
1004 vecIROCMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,0,63,0);
1005 vecOROCMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,64,159,0);
1006 vecOROC0Max[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,64,128,0);
1007 vecOROC1Max[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,129,159,0);
1009 vecAllTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,0);
1010 vecIROCTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,0);
1011 vecOROCTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,0);
1012 vecOROC0Tot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,0);
1013 vecOROC1Tot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,0);
1015 vecAllTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,0,2,1);
1016 vecIROCTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,0,2,1);
1017 vecOROCTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,0,2,1);
1018 vecOROC0TotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,0,2,1);
1019 vecOROC1TotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,0,2,1);
1021 vecAllTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,4,2,1);
1022 vecIROCTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,4,2,1);
1023 vecOROCTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,4,2,1);
1024 vecOROC0TotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,4,2,1);
1025 vecOROC1TotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,4,2,1);
1027 vecAllTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,5,2,1);
1028 vecIROCTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,5,2,1);
1029 vecOROCTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,5,2,1);
1030 vecOROC0TotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,5,2,1);
1031 vecOROC1TotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,5,2,1);
1033 vecAllTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,1,2,0);
1034 vecIROCTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,1,2,0);
1035 vecOROCTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,1,2,0);
1036 vecOROC0TotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,1,2,0);
1037 vecOROC1TotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,1,2,0);
1039 vecAllTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,6,2,1);
1040 vecIROCTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,6,2,1);
1041 vecOROCTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,6,2,1);
1042 vecOROC0TotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,6,2,1);
1043 vecOROC1TotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,6,2,1);
1045 vecAllTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,8,2,1);
1046 vecIROCTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,8,2,1);
1047 vecOROCTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,8,2,1);
1048 vecOROC0TotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,8,2,1);
1049 vecOROC1TotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,8,2,1);
1050 truncUp[ifrac]=fracUp;
1051 truncDown[ifrac]=fracDown;
1058 if (((isKaon||isProton||isPion||isElectron||isMIP||isMuon)&&(!isHighPt)) && dedxDef>0) {
1060 Int_t ncont = vertex->GetNContributors();
1061 for (Int_t ipad=0; ipad<3; ipad++){
1062 // histogram with enahanced phi granularity - to make gain phi maps
1063 Double_t xxx[4]={0,gangle[ipad],1./dedxDef, static_cast<Double_t>(ipad*2+((side>0)?0:1))};
1065 if (ipad==0) {xxx[0]=vecIROCTot[4]/medianMIP0; nclR=ncl21IROC/63.;}
1066 if (ipad==1) {xxx[0]=vecOROC0Tot[4]/medianMIP0;nclR=ncl21OROC0/63.;}
1067 if (ipad==2) {xxx[0]=vecOROC1Tot[4]/medianMIP0;nclR=ncl21OROC1/32.;}
1069 if (xxx[0]>0) xxx[0]=1./xxx[0];
1070 if (TMath::Abs(langle[ipad])<0.25&&nclR>0.4) fHistdEdxMap->Fill(xxx);
1072 for (Int_t ipad=0; ipad<3; ipad++){
1074 // this are histogram to define overall main gain correction
1075 // Maybe dead end - we can not put all info which can be used into the THnSparse
1076 // It is keeped there for educational point of view
1078 Double_t xxx[6]={0,1./dedxDef, TMath::Abs(langle[ipad]), TMath::Abs(tgl), static_cast<Double_t>(ncont), static_cast<Double_t>(ipad) };
1079 if (ipad==0) {xxx[0]=vecIROCTot[4]/medianMIP0;}
1080 if (ipad==1) {xxx[0]=vecOROC0Tot[4]/medianMIP0;}
1081 if (ipad==2) {xxx[0]=vecOROC1Tot[4]/medianMIP0;}
1083 if (xxx[0]>0) xxx[0]=1./xxx[0];
1084 if (xxx[0]>0) fHistdEdxTot->Fill(xxx);
1085 if (ipad==0) {xxx[0]=vecIROCMax[4]/medianMIP0;}
1086 if (ipad==1) {xxx[0]=vecOROC0Max[4]/medianMIP0;}
1087 if (ipad==2) {xxx[0]=vecOROC1Max[4]/medianMIP0;}
1089 if (xxx[0]>0) xxx[0]=1./xxx[0];
1090 fHistdEdxMax->Fill(xxx);
1094 // Downscale selected tracks before filling the tree
1096 Bool_t isSelected = kFALSE;
1098 if (isKaon||isProton||isPion||isElectron||isMIP||isMuon) isSelected=kTRUE;
1100 if (!isSelected) isHighPt = ((TMath::Power(1/track->Pt(),4)*gRandom->Rndm())<0.005);
1101 //if (counter>kMax && ((1/track->Pt()*gRandom->Rndm())>kMax/counter)) return;
1102 isSelected|=isHighPt;
1106 // Equalize statistic in BB bins - special care about pions
1107 Int_t entriesBB = (Int_t)hisBB->GetEntries();
1108 if ((isElectron==0 &&isMuon==0 && p<2.) && entriesBB>20 &&dedxDef>0.01){
1109 Int_t bin = hisBB->GetXaxis()->FindBin(1./dedxDef);
1110 Double_t cont = hisBB->GetBinContent(bin);
1111 Double_t mean =(entriesBB)/20.;
1112 if ((isPion>0) && gRandom->Rndm()*cont > 0.1*mean) return;
1113 if ((isPion==0) && gRandom->Rndm()*cont > 0.25*mean) return;
1115 if (!isSelected) return;
1116 if (dedxDef>0.01) hisBB->Fill(1./dedxDef);
1118 if (isHighPt) counterHPT++;
1121 TTreeSRedirector * pcstream = GetDebugStreamer();
1122 Double_t ptrel0 = AliTPCcalibDB::GetPTRelative(fTime,fRun,0);
1123 Double_t ptrel1 = AliTPCcalibDB::GetPTRelative(fTime,fRun,1);
1124 Int_t sectorIn = Int_t(18+9*(trackIn->GetAlpha()/TMath::Pi()))%18;
1125 Int_t sectorOut = Int_t(18+9*(trackOut->GetAlpha()/TMath::Pi()))%18;
1128 (*pcstream)<<"dump"<<
1129 "vertex.="<<vertex<<
1133 "sectorIn="<<sectorIn<<
1134 "sectorOut="<<sectorOut<<
1138 "isProton="<<isProton<<
1141 "isElectron="<<isElectron<<
1143 "isHighPt="<<isHighPt<<
1145 "dedxDef="<<dedxDef<<
1146 "dedxDefPion="<<dedxDefPion<<
1147 "dedxDefKaon="<<dedxDefKaon<<
1148 "dedxDefProton="<<dedxDefProton<<
1155 "ncl20All="<<ncl20All<<
1156 "ncl20OROC="<<ncl20OROC<<
1157 "ncl20IROC="<<ncl20IROC<<
1158 "ncl20OROC0="<<ncl20OROC0<<
1159 "ncl20OROC1="<<ncl20OROC1<<
1161 "ncl20All4="<<ncl20All4<<
1162 "ncl20OROC4="<<ncl20OROC4<<
1163 "ncl20IROC4="<<ncl20IROC4<<
1164 "ncl20OROC04="<<ncl20OROC04<<
1165 "ncl20OROC14="<<ncl20OROC14<<
1167 "ncl20All3="<<ncl20All3<<
1168 "ncl20OROC3="<<ncl20OROC3<<
1169 "ncl20IROC3="<<ncl20IROC3<<
1170 "ncl20OROC03="<<ncl20OROC03<<
1171 "ncl20OROC13="<<ncl20OROC13<<
1173 "ncl21All="<<ncl21All<<
1174 "ncl21OROC="<<ncl21OROC<<
1175 "ncl21IROC="<<ncl21IROC<<
1176 "ncl21OROC0="<<ncl21OROC0<<
1177 "ncl21OROC1="<<ncl21OROC1<<
1179 "langle0="<<langle[0]<<
1180 "langle1="<<langle[1]<<
1181 "langle2="<<langle[2]<<
1182 "gangle0="<<gangle[0]<< //global angle phi IROC
1183 "gangle1="<<gangle[1]<< // OROC medium
1184 "gangle2="<<gangle[2]<< // OROC long
1193 //"trackIn.="<<trackIn<<
1194 //"trackOut.="<<trackOut<<
1195 //"tpcOut.="<<tpcOut<<
1197 "medianMIP0="<<medianMIP0<< // median MIP position as estimated from the array of (not only) "MIPS"
1199 "truncUp.="<<&truncUp<<
1200 "truncDown.="<<&truncDown<<
1201 "vecAllMax.="<<&vecAllMax<<
1202 "vecIROCMax.="<<&vecIROCMax<<
1203 "vecOROCMax.="<<&vecOROCMax<<
1204 "vecOROC0Max.="<<&vecOROC0Max<<
1205 "vecOROC1Max.="<<&vecOROC1Max<<
1207 "vecAllTot.="<<&vecAllTot<<
1208 "vecIROCTot.="<<&vecIROCTot<<
1209 "vecOROCTot.="<<&vecOROCTot<<
1210 "vecOROC0Tot.="<<&vecOROC0Tot<<
1211 "vecOROC1Tot.="<<&vecOROC1Tot<<
1213 "vecAllTotLog.="<<&vecAllTotLog<<
1214 "vecIROCTotLog.="<<&vecIROCTotLog<<
1215 "vecOROCTotLog.="<<&vecOROCTotLog<<
1216 "vecOROC0TotLog.="<<&vecOROC0TotLog<<
1217 "vecOROC1TotLog.="<<&vecOROC1TotLog<<
1219 "vecAllTotUp.="<<&vecAllTotUp<<
1220 "vecIROCTotUp.="<<&vecIROCTotUp<<
1221 "vecOROCTotUp.="<<&vecOROCTotUp<<
1222 "vecOROC0TotUp.="<<&vecOROC0TotUp<<
1223 "vecOROC1TotUp.="<<&vecOROC1TotUp<<
1225 "vecAllTotDown.="<<&vecAllTotDown<<
1226 "vecIROCTotDown.="<<&vecIROCTotDown<<
1227 "vecOROCTotDown.="<<&vecOROCTotDown<<
1228 "vecOROC0TotDown.="<<&vecOROC0TotDown<<
1229 "vecOROC1TotDown.="<<&vecOROC1TotDown<<
1231 "vecAllTotRMS.="<<&vecAllTotRMS<<
1232 "vecIROCTotRMS.="<<&vecIROCTotRMS<<
1233 "vecOROCTotRMS.="<<&vecOROCTotRMS<<
1234 "vecOROC0TotRMS.="<<&vecOROC0TotRMS<<
1235 "vecOROC1TotRMS.="<<&vecOROC1TotRMS<<
1237 "vecAllTotM2.="<<&vecAllTotM2<<
1238 "vecIROCTotM2.="<<&vecIROCTotM2<<
1239 "vecOROCTotM2.="<<&vecOROCTotM2<<
1240 "vecOROC0TotM2.="<<&vecOROC0TotM2<<
1241 "vecOROC1TotM2.="<<&vecOROC1TotM2<<
1243 "vecAllTotMS.="<<&vecAllTotMS<<
1244 "vecIROCTotMS.="<<&vecIROCTotMS<<
1245 "vecOROCTotMS.="<<&vecOROCTotMS<<
1246 "vecOROC0TotMS.="<<&vecOROC0TotMS<<
1247 "vecOROC1TotMS.="<<&vecOROC1TotMS<<
1255 void AliTPCcalibGainMult::ProcessV0s(AliESDEvent * event){
1257 // Select the K0s and gamma - and sign daughter products
1259 TTreeSRedirector * pcstream = GetDebugStreamer();
1260 AliKFParticle::SetField(event->GetMagneticField());
1261 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1263 //Printf("ERROR: esdFriend not available");
1266 if (esdFriend->TestSkipBit()) return;
1269 static const TDatabasePDG *pdg = TDatabasePDG::Instance();
1270 const Double_t kChi2Cut=5;
1271 const Double_t kMinR=2;
1272 const Int_t kMinNcl=110;
1273 const Double_t kMinREl=5;
1274 const Double_t kMaxREl=70;
1276 Int_t nv0 = event->GetNumberOfV0s();
1277 AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex();
1278 AliKFVertex kfvertex=*vertex;
1280 for (Int_t iv0=0;iv0<nv0;iv0++){
1281 AliESDv0 *v0 = event->GetV0(iv0);
1283 if (v0->GetOnFlyStatus()<0.5) continue;
1284 if (v0->GetPindex()<0) continue;
1285 if (v0->GetNindex()<0) continue;
1286 if (TMath::Max(v0->GetPindex(), v0->GetNindex())>event->GetNumberOfTracks()) continue;
1289 AliExternalTrackParam pp=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamP())):(*(v0->GetParamN()));
1290 AliExternalTrackParam pn=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamN())):(*(v0->GetParamP()));
1291 AliKFParticle kfp1( pp, 211 );
1292 AliKFParticle kfp2( pn, -211 );
1294 AliKFParticle *v0KFK0 = new AliKFParticle(kfp1,kfp2);
1295 AliKFParticle *v0KFK0CV = new AliKFParticle(*v0KFK0);
1296 v0KFK0CV->SetProductionVertex(kfvertex);
1297 v0KFK0CV->TransportToProductionVertex();
1298 Double_t chi2K0 = v0KFK0CV->GetChi2();
1299 if (chi2K0>kChi2Cut) continue;
1300 if (v0->GetRr()<kMinR) continue;
1301 Bool_t isOKC=TMath::Max(v0->GetCausalityP()[0],v0->GetCausalityP()[1])<0.7&&TMath::Min(v0->GetCausalityP()[2],v0->GetCausalityP()[3])>0.2;
1303 Double_t effMass22=v0->GetEffMass(2,2);
1304 Double_t effMass42=v0->GetEffMass(4,2);
1305 Double_t effMass24=v0->GetEffMass(2,4);
1306 Double_t effMass00=v0->GetEffMass(0,0);
1307 AliKFParticle *v0KFK0CVM = new AliKFParticle(*v0KFK0CV);
1308 v0KFK0CVM->SetMassConstraint(pdg->GetParticle("K_S0")->Mass());
1309 Bool_t isV0= kFALSE;
1311 Double_t d22 = TMath::Abs(effMass22-pdg->GetParticle("K_S0")->Mass());
1312 Double_t d42 = TMath::Abs(effMass42-pdg->GetParticle("Lambda0")->Mass());
1313 Double_t d24 = TMath::Abs(effMass24-pdg->GetParticle("Lambda0")->Mass());
1314 Double_t d00 = TMath::Abs(effMass00);
1316 Bool_t isKaon = d22<0.01 && d22< 0.3 * TMath::Min(TMath::Min(d42,d24),d00);
1317 Bool_t isLambda = d42<0.01 && d42< 0.3 * TMath::Min(TMath::Min(d22,d24),d00);
1318 Bool_t isAntiLambda= d24<0.01 && d24< 0.3 * TMath::Min(TMath::Min(d22,d42),d00);
1319 Bool_t isGamma = d00<0.02 && d00< 0.3 * TMath::Min(TMath::Min(d42,d24),d22);
1321 if (isGamma && (isKaon||isLambda||isAntiLambda)) continue;
1322 if (isLambda && (isKaon||isGamma||isAntiLambda)) continue;
1323 if (isKaon && (isLambda||isGamma||isAntiLambda)) continue;
1324 Double_t sign= v0->GetParamP()->GetSign()* v0->GetParamN()->GetSign();
1325 if (sign>0) continue;
1326 isV0=isKaon||isLambda||isAntiLambda||isGamma;
1327 if (!(isV0)) continue;
1328 if (isGamma&&v0->GetRr()<kMinREl) continue;
1329 if (isGamma&&v0->GetRr()>kMaxREl) continue;
1330 if (!isOKC) continue;
1332 Int_t pindex = (v0->GetParamP()->GetSign()>0) ? v0->GetPindex() : v0->GetNindex();
1333 Int_t nindex = (v0->GetParamP()->GetSign()>0) ? v0->GetNindex() : v0->GetPindex();
1334 AliESDtrack * trackP = event->GetTrack(pindex);
1335 AliESDtrack * trackN = event->GetTrack(nindex);
1336 if (!trackN) continue;
1337 if (!trackP) continue;
1338 Int_t nclP= (Int_t)trackP->GetTPCClusterInfo(2,1);
1339 Int_t nclN= (Int_t)trackN->GetTPCClusterInfo(2,1);
1340 if (TMath::Min(nclP,nclN)<kMinNcl) continue;
1341 Double_t eta = TMath::Max(TMath::Abs(trackP->Eta()), TMath::Abs(trackN->Eta()));
1342 if (TMath::Abs(eta)>1) continue;
1345 AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(pindex);
1346 AliESDfriendTrack *friendTrackN = esdFriend->GetTrack(nindex);
1347 if (!friendTrackP) continue;
1348 if (!friendTrackN) continue;
1349 TObject *calibObject;
1350 AliTPCseed *seedP = 0;
1351 AliTPCseed *seedN = 0;
1352 for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
1353 if ((seedP=dynamic_cast<AliTPCseed*>(calibObject))) break;
1355 for (Int_t l=0;(calibObject=friendTrackN->GetCalibObject(l));++l) {
1356 if ((seedN=dynamic_cast<AliTPCseed*>(calibObject))) break;
1359 if ( TMath::Abs((trackP->GetTPCsignal()/(trackN->GetTPCsignal()+0.0001)-1)>0.3)) continue;
1361 if (isGamma) (*fPIDMatrix)(pindex, 0)+=2;
1362 if (isGamma) (*fPIDMatrix)(nindex, 0)+=2;
1364 if (isKaon) (*fPIDMatrix)(pindex, 2)+=2;
1365 if (isKaon) (*fPIDMatrix)(nindex, 2)+=2;
1369 (*pcstream)<<"v0s"<<
1370 "isGamma="<<isGamma<<
1372 "isLambda="<<isLambda<<
1373 "isAntiLambda="<<isAntiLambda<<
1375 "trackP.="<<trackP<<
1376 "trackN.="<<trackN<<
1386 void AliTPCcalibGainMult::ProcessCosmic(const AliESDEvent * event) {
1388 // Find cosmic pairs trigger by random trigger
1391 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1392 AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
1394 AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
1395 AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
1396 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1397 const Double_t kMinPt=4;
1398 const Double_t kMinPtMax=0.8;
1399 const Double_t kMinNcl=159*0.5;
1400 const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
1401 Int_t ntracks=event->GetNumberOfTracks();
1402 const Double_t kMaxImpact=80;
1403 // Float_t dcaTPC[2]={0,0};
1404 // Float_t covTPC[3]={0,0,0};
1406 UInt_t specie = event->GetEventSpecie(); // skip laser events
1407 if (specie==AliRecoParam::kCalib) return;
1410 for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
1411 AliESDtrack *track0 = event->GetTrack(itrack0);
1412 if (!track0) continue;
1413 if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
1415 if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue;
1416 if (track0->GetTPCncls()<kMinNcl) continue;
1417 if (TMath::Abs(track0->GetY())<2*kMaxDelta[0]) continue;
1418 if (TMath::Abs(track0->GetY())>kMaxImpact) continue;
1419 if (track0->GetKinkIndex(0)>0) continue;
1420 const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
1423 for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
1424 AliESDtrack *track1 = event->GetTrack(itrack1);
1425 if (!track1) continue;
1426 if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
1427 if (track1->GetKinkIndex(0)>0) continue;
1428 if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue;
1429 if (track1->GetTPCncls()<kMinNcl) continue;
1430 if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
1431 if (TMath::Abs(track1->GetY())<2*kMaxDelta[0]) continue;
1432 if (TMath::Abs(track1->GetY())>kMaxImpact) continue;
1434 const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
1436 Bool_t isPair=kTRUE;
1437 for (Int_t ipar=0; ipar<5; ipar++){
1438 if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
1439 if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
1441 if (!isPair) continue;
1442 if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
1443 //delta with correct sign
1444 if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
1445 if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
1446 if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
1447 if (!isPair) continue;
1448 TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1449 Int_t eventNumber = event->GetEventNumberInFile();
1450 Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(itrack0)!=0):0;
1451 Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);
1452 AliInfo(Form("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS));
1455 TTreeSRedirector * pcstream = GetDebugStreamer();
1456 Int_t ntracksSPD = vertexSPD->GetNContributors();
1457 Int_t ntracksTPC = vertexTPC->GetNContributors();
1460 (*pcstream)<<"cosmicPairsAll"<<
1461 "run="<<fRun<< // run number
1462 "event="<<fEvent<< // event number
1463 "time="<<fTime<< // time stamp of event
1464 "trigger="<<fTrigger<< // trigger
1465 "triggerClass="<<&fTriggerClass<< // trigger
1466 "bz="<<fMagF<< // magnetic field
1468 "nSPD="<<ntracksSPD<<
1469 "nTPC="<<ntracksTPC<<
1470 "vSPD.="<<vertexSPD<< //primary vertex -SPD
1471 "vTPC.="<<vertexTPC<< //primary vertex -TPC
1472 "t0.="<<track0<< //track0
1473 "t1.="<<track1<< //track1
1477 AliESDfriendTrack *friendTrack0 = esdFriend->GetTrack(itrack0);
1478 if (!friendTrack0) continue;
1479 AliESDfriendTrack *friendTrack1 = esdFriend->GetTrack(itrack1);
1480 if (!friendTrack1) continue;
1481 TObject *calibObject;
1482 AliTPCseed *seed0 = 0;
1483 AliTPCseed *seed1 = 0;
1485 for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) {
1486 if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1488 for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) {
1489 if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1493 (*pcstream)<<"cosmicPairs"<<
1494 "run="<<fRun<< // run number
1495 "event="<<fEvent<< // event number
1496 "time="<<fTime<< // time stamp of event
1497 "trigger="<<fTrigger<< // trigger
1498 "triggerClass="<<&fTriggerClass<< // trigger
1499 "bz="<<fMagF<< // magnetic field
1501 "nSPD="<<ntracksSPD<<
1502 "nTPC="<<ntracksTPC<<
1503 "vSPD.="<<vertexSPD<< //primary vertex -SPD
1504 "vTPC.="<<vertexTPC<< //primary vertex -TPC
1505 "t0.="<<track0<< //track0
1506 "t1.="<<track1<< //track1
1507 "ft0.="<<friendTrack0<< //track0
1508 "ft1.="<<friendTrack1<< //track1
1509 "s0.="<<seed0<< //track0
1510 "s1.="<<seed1<< //track1
1513 if (!seed0) continue;
1514 if (!seed1) continue;
1515 Int_t nclA0=0, nclC0=0; // number of clusters
1516 Int_t nclA1=0, nclC1=0; // number of clusters
1518 for (Int_t irow=0; irow<159; irow++){
1519 AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1520 AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1522 if (cluster0->GetQ()>0 && cluster0->GetDetector()%36<18) nclA0++;
1523 if (cluster0->GetQ()>0 && cluster0->GetDetector()%36>=18) nclC0++;
1526 if (cluster1->GetQ()>0 && cluster1->GetDetector()%36<18) nclA1++;
1527 if (cluster1->GetQ()>0 && cluster1->GetDetector()%36>=18) nclC1++;
1530 Int_t cosmicType=0; // types of cosmic topology
1531 if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
1532 if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
1533 if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
1534 if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
1535 if (cosmicType<2) continue; // use only crossing tracks
1537 Double_t deltaTimeCluster=0;
1538 deltaTimeCluster=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
1539 if (nclA0>nclC0) deltaTimeCluster*=-1; // if A side track
1541 for (Int_t irow=0; irow<159; irow++){
1542 AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1543 if (cluster0 &&cluster0->GetX()>10){
1544 Double_t x0[3]={ static_cast<Double_t>(cluster0->GetRow()),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
1545 Int_t index0[1]={cluster0->GetDetector()};
1546 transform->Transform(x0,index0,0,1);
1547 cluster0->SetX(x0[0]);
1548 cluster0->SetY(x0[1]);
1549 cluster0->SetZ(x0[2]);
1552 AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1553 if (cluster1&&cluster1->GetX()>10){
1554 Double_t x1[3]={ static_cast<Double_t>(cluster1->GetRow()),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
1555 Int_t index1[1]={cluster1->GetDetector()};
1556 transform->Transform(x1,index1,0,1);
1557 cluster1->SetX(x1[0]);
1558 cluster1->SetY(x1[1]);
1559 cluster1->SetZ(x1[2]);
1565 (*fPIDMatrix)(itrack0,1)+=4; //
1566 (*fPIDMatrix)(itrack1,1)+=4; //
1574 void AliTPCcalibGainMult::ProcessKinks(const AliESDEvent * event){
1578 AliKFParticle::SetField(event->GetMagneticField());
1579 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1581 //Printf("ERROR: esdFriend not available");
1584 // if (esdFriend->TestSkipBit()) return;
1587 const Double_t kChi2Cut=10;
1588 const Double_t kMinR=100;
1589 const Double_t kMaxR=230;
1590 const Int_t kMinNcl=110;
1592 Int_t nkinks = event->GetNumberOfKinks();
1593 AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex();
1594 AliKFVertex kfvertex=*vertex;
1595 TTreeSRedirector * pcstream = GetDebugStreamer();
1597 for (Int_t ikink=0;ikink<nkinks;ikink++){
1598 AliESDkink *kink = event->GetKink(ikink);
1599 if (!kink) continue;
1600 if (kink->GetIndex(0)<0) continue;
1601 if (kink->GetIndex(1)<0) continue;
1602 if (TMath::Max(kink->GetIndex(1), kink->GetIndex(0))>event->GetNumberOfTracks()) continue;
1605 AliExternalTrackParam pd=kink->RefParamDaughter();
1606 AliExternalTrackParam pm=kink->RefParamMother();
1607 AliKFParticle kfpd( pd, 211 );
1608 AliKFParticle kfpm( pm, -13 );
1610 AliKFParticle *v0KF = new AliKFParticle(kfpm,kfpd);
1611 v0KF->SetVtxGuess(kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]);
1612 Double_t chi2 = v0KF->GetChi2();
1613 AliESDtrack * trackM = event->GetTrack(kink->GetIndex(0));
1614 AliESDtrack * trackD = event->GetTrack(kink->GetIndex(1));
1615 if (!trackM) continue;
1616 if (!trackD) continue;
1617 Int_t nclM= (Int_t)trackM->GetTPCClusterInfo(2,1);
1618 Int_t nclD= (Int_t)trackD->GetTPCClusterInfo(2,1);
1619 Double_t eta = TMath::Max(TMath::Abs(trackM->Eta()), TMath::Abs(trackD->Eta()));
1620 Double_t kx= v0KF->GetX();
1621 Double_t ky= v0KF->GetY();
1622 Double_t kz= v0KF->GetZ();
1623 Double_t ex= v0KF->GetErrX();
1624 Double_t ey= v0KF->GetErrY();
1625 Double_t ez= v0KF->GetErrZ();
1627 Double_t radius=TMath::Sqrt(kx*kx+ky*ky);
1628 Double_t alpha=TMath::ATan2(ky,kx);
1629 if (!pd.Rotate(alpha)) continue;
1630 if (!pm.Rotate(alpha)) continue;
1631 if (!pd.PropagateTo(radius,event->GetMagneticField())) continue;
1632 if (!pm.PropagateTo(radius,event->GetMagneticField())) continue;
1633 Double_t pos[2]={0,kz};
1634 Double_t cov[3]={ex*ex+ey*ey,0,ez*ez};
1639 (*pcstream)<<"kinks"<<
1644 "trackM.="<<trackM<<
1645 "trackD.="<<trackD<<
1646 "pm.="<<&pm<< //updated parameters
1647 "pd.="<<&pd<< // updated parameters
1653 TCut cutQ="chi2<10&&kink.fRr>90&&kink.fRr<220";
1654 TCut cutRD="20*sqrt(pd.fC[14])<abs(pd.fP[4])&&trackD.fTPCsignal>10&&trackD.fTPCsignalN>50";
1658 if (chi2>kChi2Cut) continue;
1659 if (kink->GetR()<kMinR) continue;
1660 if (kink->GetR()>kMaxR) continue;
1661 if ((nclM+nclD)<kMinNcl) continue;
1662 if (TMath::Abs(eta)>1) continue;
1665 AliESDfriendTrack *friendTrackM = esdFriend->GetTrack(kink->GetIndex(0));
1666 AliESDfriendTrack *friendTrackD = esdFriend->GetTrack(kink->GetIndex(1));
1667 if (!friendTrackM) continue;
1668 if (!friendTrackD) continue;
1669 TObject *calibObject;
1670 AliTPCseed *seedM = 0;
1671 AliTPCseed *seedD = 0;
1672 for (Int_t l=0;(calibObject=friendTrackM->GetCalibObject(l));++l) {
1673 if ((seedM=dynamic_cast<AliTPCseed*>(calibObject))) break;
1675 for (Int_t l=0;(calibObject=friendTrackD->GetCalibObject(l));++l) {
1676 if ((seedD=dynamic_cast<AliTPCseed*>(calibObject))) break;
1681 void AliTPCcalibGainMult::DumpHPT(const AliESDEvent * event){
1683 // Function to select the HPT tracks and events
1684 // It is used to select event with HPT - list used later for the raw data downloading
1685 // - and reconstruction
1686 // Not actualy used for the calibration of the data
1688 TTreeSRedirector * pcstream = GetDebugStreamer();
1689 AliKFParticle::SetField(event->GetMagneticField());
1690 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1692 //Printf("ERROR: esdFriend not available");
1695 if (esdFriend->TestSkipBit()) return;
1697 Int_t ntracks=event->GetNumberOfTracks();
1699 for (Int_t i=0;i<ntracks;++i) {
1701 AliESDtrack *track = event->GetTrack(i);
1702 if (!track) continue;
1703 if (track->Pt()<4) continue;
1704 UInt_t status = track->GetStatus();
1706 AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam();
1707 if (!trackIn) continue;
1708 if ((status&AliESDtrack::kTPCrefit)==0) continue;
1709 if ((status&AliESDtrack::kITSrefit)==0) continue;
1710 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
1711 if (!friendTrack) continue;
1712 AliExternalTrackParam * itsOut = (AliExternalTrackParam *)(friendTrack->GetITSOut());
1713 if (!itsOut) continue;
1714 AliExternalTrackParam * itsOut2 = (AliExternalTrackParam *)(friendTrack->GetITSOut()->Clone());
1715 AliExternalTrackParam * tpcIn2 = (AliExternalTrackParam *)(trackIn->Clone());
1716 if (!itsOut2->Rotate(trackIn->GetAlpha())) continue;
1717 //Double_t xmiddle=0.5*(itsOut2->GetX()+tpcIn2->GetX());
1718 Double_t xmiddle=(itsOut2->GetX());
1719 if (!itsOut2->PropagateTo(xmiddle,event->GetMagneticField())) continue;
1720 if (!tpcIn2->PropagateTo(xmiddle,event->GetMagneticField())) continue;
1722 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
1723 if (!tpcInner) continue;
1724 tpcInner->Rotate(track->GetAlpha());
1725 tpcInner->PropagateTo(track->GetX(),event->GetMagneticField());
1729 AliExternalTrackParam * tpcInnerC = (AliExternalTrackParam *)(track->GetTPCInnerParam()->Clone());
1730 if (!tpcInnerC) continue;
1731 tpcInnerC->Rotate(track->GetAlpha());
1732 tpcInnerC->PropagateTo(track->GetX(),event->GetMagneticField());
1733 Double_t dz[2],cov[3];
1734 AliESDVertex *vtx= (AliESDVertex *)event->GetPrimaryVertex();
1736 if (!tpcInnerC->PropagateToDCA(vtx, event->GetMagneticField(), 3, dz, cov)) continue;
1737 Double_t covar[6]; vtx->GetCovMatrix(covar);
1738 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
1739 Double_t c[3]={covar[2],0.,covar[5]};
1741 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
1742 tpcInnerC->Update(p,c);
1745 (*pcstream)<<"hpt"<<
1751 "tpcInner.="<<tpcInner<<
1752 "tpcInnerC.="<<tpcInnerC<<
1766 void AliTPCcalibGainMult::ProcessTOF(const AliESDEvent * event){
1768 // 1. Loop over tracks
1770 // 3. Sign positivelly identified tracks
1772 const Double_t kMaxDelta=1000;
1773 const Double_t kOrbit=50000; // distance in the time beween two orbits in the TOF time units - 50000=50 ns
1774 const Double_t kMaxD=20;
1775 const Double_t kRMS0=200;
1776 const Double_t kMaxDCAZ=10;
1777 AliESDVertex *vtx= (AliESDVertex *)event->GetPrimaryVertex();
1779 TTreeSRedirector * pcstream = GetDebugStreamer();
1780 AliKFParticle::SetField(event->GetMagneticField());
1781 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1783 //Printf("ERROR: esdFriend not available");
1786 //if (esdFriend->TestSkipBit()) return;
1788 Int_t ntracks=event->GetNumberOfTracks();
1790 Double_t deltaTPion[10000];
1791 Double_t medianT0=0;
1796 // Get Median time for pion hypothesy
1798 for (Int_t iter=0; iter<3; iter++){
1800 for (Int_t i=0;i<ntracks;++i) {
1802 AliESDtrack *track = event->GetTrack(i);
1803 if (!track) continue;
1804 if (!track->IsOn(AliESDtrack::kTIME)) continue;
1805 if (TMath::Abs(track->GetZ())>kMaxDCAZ) continue; // remove overlaped events
1806 if (TMath::Abs(track->GetTOFsignalDz())>kMaxD) continue;
1807 Double_t times[1000];
1808 track->GetIntegratedTimes(times);
1809 Int_t norbit=TMath::Nint((track->GetTOFsignal()-times[2])/kOrbit);
1810 Double_t torbit=norbit*kOrbit;
1811 if (iter==1 &&TMath::Abs(times[2]-times[3])<3*rms) continue; // skip umbigous points - kaon pion
1815 for (Int_t j=3; j<5; j++)
1816 if (TMath::Abs(track->GetTOFsignal()-times[j]-torbit-medianT0)<TMath::Abs(track->GetTOFsignal()-times[j]-torbit-medianT0)) indexBest=j;
1819 if (iter>0) if (TMath::Abs(track->GetTOFsignal()-times[indexBest]-torbit-medianT0)>3*(kRMS0+rms)) continue;
1820 if (iter>0) if (TMath::Abs(track->GetTOFsignal()-times[indexBest]-torbit-medianT0)>kMaxDelta) continue;
1821 deltaTPion[counter]=track->GetTOFsignal()-times[indexBest]-torbit;
1824 if (counter<2) return;
1825 medianT0=TMath::Median(counter,deltaTPion);
1826 meanT0=TMath::Median(counter,deltaTPion);
1827 rms=TMath::RMS(counter,deltaTPion);
1829 if (counter<3) return;
1833 for (Int_t i=0;i<ntracks;++i) {
1835 AliESDtrack *track = event->GetTrack(i);
1836 if (!track) continue;
1837 if (!track->IsOn(AliESDtrack::kTIME)) continue;
1838 if (TMath::Abs(track->GetZ())>kMaxDCAZ) continue; //remove overlapped events
1839 if (TMath::Abs(track->GetTOFsignalDz())>kMaxD) continue;
1840 Double_t times[1000];
1841 track->GetIntegratedTimes(times);
1842 Int_t norbit=TMath::Nint((track->GetTOFsignal()-times[2])/kOrbit);
1843 Double_t torbit=norbit*kOrbit;
1844 if (rms<=0) continue;
1846 Double_t tPion = (track->GetTOFsignal()-times[2]-medianT0-torbit);
1847 Double_t tKaon = (track->GetTOFsignal()-times[3]-medianT0-torbit);
1848 Double_t tProton= (track->GetTOFsignal()-times[4]-medianT0-torbit);
1849 Double_t tElectron= (track->GetTOFsignal()-times[0]-medianT0-torbit);
1851 Bool_t isPion = (TMath::Abs(tPion/rms)<6) && TMath::Abs(tPion)<(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tProton))-rms);
1852 Bool_t isKaon = (TMath::Abs(tKaon/rms)<3) && TMath::Abs(tKaon)<0.2*(TMath::Min(TMath::Abs(tPion), TMath::Abs(tProton))-3*rms);
1853 Bool_t isProton = (TMath::Abs(tProton/rms)<6) && TMath::Abs(tProton)<0.5*(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tPion))-rms);
1854 Bool_t isElectron = (TMath::Abs(tElectron/rms)<6) && TMath::Abs(tElectron)<0.2*(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tPion))-rms) &&TMath::Abs(tElectron)<0.5*(TMath::Abs(tPion)-rms);
1856 if (isPion) (*fPIDMatrix)(i,2)+=1;
1857 if (isKaon) (*fPIDMatrix)(i,3)+=1;
1858 if (isProton) (*fPIDMatrix)(i,4)+=1;
1859 // if (isElectron) (*fPIDMatrix)(i,0)+=1;
1862 // debug streamer to dump the information
1863 (*pcstream)<<"tof"<<
1866 "isProton="<<isProton<<
1867 "isElectron="<<isElectron<<
1869 "counter="<<counter<<
1872 "medianT0="<<medianT0<<
1882 tof->SetAlias("isProton","(abs(track.fTOFsignal-track.fTrackTime[4]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)-rmsT0))");
1883 tof->SetAlias("isPion","(abs(track.fTOFsignal-track.fTrackTime[2]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)-rmsT0))");
1884 tof->SetAlias("isKaon","(abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[2]-medianT0-torbit)-rmsT0))&&(abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[4]-medianT0-torbit)-rmsT0))");
1891 TGraphErrors* AliTPCcalibGainMult::GetGainPerChamber(Int_t padRegion/*=1*/, Bool_t plotQA/*=kFALSE*/)
1894 // Extract gain variations per chamger for 'padRegion'
1897 if (padRegion<0||padRegion>2) return 0x0;
1899 if (!fHistGainSector) return NULL;
1900 if (!fHistGainSector->GetAxis(2)) return NULL;
1901 fHistGainSector->GetAxis(2)->SetRangeUser(padRegion,padRegion);
1902 TH2D * histGainSec = fHistGainSector->Projection(0,1);
1903 // TH1D* proj=fHistGainSector->Projection(0);
1904 // Double_t max=proj->GetBinCenter(proj->GetMaximumBin());
1908 // TF1 fg("gaus","gaus",histGainSec->GetYaxis()->GetXmin()+1,histGainSec->GetYaxis()->GetXmax()-1);
1909 // histGainSec->FitSlicesY(&fg, 0, -1, 0, "QNR", &arr);
1910 histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
1911 TH1D * meanGainSec = (TH1D*)arr.At(1);
1912 Double_t gainsIROC[36]={0.};
1913 Double_t gainsOROC[36]={0.};
1914 Double_t gains[72]={0.};
1915 Double_t gainsErr[72]={0.};
1916 TGraphErrors *gr=new TGraphErrors(36);
1918 for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) {
1919 TH1D *slice=histGainSec->ProjectionY("_py",isec,isec);
1920 Double_t max=slice->GetBinCenter(slice->GetMaximumBin());
1921 TF1 fg("gaus","gaus",max-30,max+30);
1922 slice->Fit(&fg,"QNR");
1923 meanGainSec->SetBinContent(isec,fg.GetParameter(1));
1924 meanGainSec->SetBinError(isec,fg.GetParError(1));
1926 // cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl;
1927 gains[isec-1] = meanGainSec->GetBinContent(isec);
1929 gainsIROC[isec-1] = meanGainSec->GetBinContent(isec);
1931 gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec);
1933 gainsErr[isec-1]=meanGainSec->GetBinError(isec);
1936 Double_t meanIroc = TMath::Median(36, gainsIROC);
1937 Double_t meanOroc = TMath::Median(36, gainsOROC);
1938 if (TMath::Abs(meanIroc)<1e-30) meanIroc=1.;
1939 if (TMath::Abs(meanOroc)<1e-30) meanOroc=1.;
1940 for(Int_t i = 0; i < 36; i++) {
1941 gains[i] /= meanIroc;
1942 gainsErr[i] /= meanIroc;
1945 for(Int_t i = 36; i < 72; i++){
1946 gains[i] /= meanOroc;
1947 gainsErr[i] /= meanOroc;
1951 for(Int_t i = 0; i < 72; i++) {
1952 if (padRegion==0 && i>35) continue;
1953 if ( (padRegion==1 || padRegion==2) && i<36) continue;
1955 if (gains[i]<1e-20 || gainsErr[i]/gains[i]>.2){
1956 AliWarning(Form("Invalid chamber gain in ROC/region: %d / %d", i, padRegion));
1961 gr->SetPoint(ipoint,i,gains[i]);
1962 gr->SetPointError(ipoint,0,gainsErr[i]);
1966 const char* names[3]={"SHORT","MEDIUM","LONG"};
1967 gr->SetNameTitle(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]),Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]));
1970 //=====================================
1971 // Do QA plotting if requested
1973 TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("cQA");
1974 if (!c) c=new TCanvas("cQA","cQA");
1978 histGainSec->DrawCopy("colz");
1979 meanGainSec->DrawCopy("same");
1980 gr->SetMarkerStyle(20);
1981 gr->SetMarkerSize(.5);
1989 TGraphErrors* AliTPCcalibGainMult::GetGainPerChamberRobust(Int_t padRegion/*=1*/, Bool_t /*plotQA=kFALSE*/)
1992 // Extract gain variations per chamger for 'padRegion'
1993 // Use Robust mean - LTM with 60 % 0 should be similar to the truncated mean 60 %
1994 if (padRegion<0||padRegion>2) return 0x0;
1995 const Int_t colors[10]={1,2,4,6};
1996 const Int_t markers[10]={21,25,22,20};
1998 if (!fHistGainSector) return NULL;
1999 if (!fHistGainSector->GetAxis(2)) return NULL;
2000 fHistGainSector->GetAxis(2)->SetRangeUser(padRegion,padRegion);
2001 if (padRegion==0) fHistGainSector->GetAxis(1)->SetRangeUser(0.0,35.);
2002 if (padRegion>0) fHistGainSector->GetAxis(1)->SetRangeUser(36.,71.);
2004 TH2D * histGainSec = fHistGainSector->Projection(0,1);
2005 TGraphErrors * gr = TStatToolkit::MakeStat1D(histGainSec, 0, 0.6,2,markers[padRegion],colors[padRegion]);
2007 Double_t median = TMath::Median(gr->GetN(),gr->GetY());
2009 for (Int_t i=0; i<gr->GetN();i++) {
2010 gr->GetY()[i]/=median;
2011 gr->SetPointError(i,gr->GetErrorX(i),gr->GetErrorY(i)/median);
2014 const char* names[3]={"SHORT","MEDIUM","LONG"};
2015 gr->SetNameTitle(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]),Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[padRegion]));
2019 // void AliTPCcalibGainMult::Terminate(){
2021 // // Terminate function
2022 // // call base terminate + Eval of fitters
2024 // Info("AliTPCcalibGainMult","Terminate");
2025 // TTreeSRedirector *pcstream = GetDebugStreamer();
2027 // TTreeStream &stream = (*pcstream)<<"dump";
2028 // TTree* tree = stream.GetTree();
2029 // if (tree) if ( tree->GetEntries()>0){
2030 // TObjArray *array = tree->GetListOfBranches();
2031 // for (Int_t i=0; i<array->GetEntries(); i++) {TBranch * br = (TBranch *)array->At(i); br->SetAddress(0);}
2032 // gDirectory=gROOT;
2033 // fdEdxTree=tree->CloneTree(10000);
2036 // AliTPCcalibBase::Terminate();