[u/mrichter/AliRoot.git] / ITS / AliITSCompareHitsRecPoints.C

Bool_t AliITSCompareHitsRecPoints(Char_t *rfn="galice.root"){
    // This macro compares the average location of the hits associated
    // with a specific track and the best match RecPoint location.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //   kTRUE if no errors were encountered, otherwise kFALSE.
    TProfile *pxg,*pyg,*pzg,*pxl,*pyl,*pzl;
    Double_t hg[3],hl[3],tof,rg[3],rl[3];
    Float_t rgf[3];

    // Dynamically link some shared libs
    if (gClassTable->GetID("AliRun") < 0) {
        gROOT->LoadMacro("loadlibs.C");
        loadlibs();
    } else {
        if(gAlice){
            delete AliRunLoader::Instance();
            delete gAlice;
            gAlice=0;
        } // end if gAlice
    }// end if gClassTable->GetID()

    gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");

    AliRunLoader *runl = AccessFile(rfn); // Set up to read in Data
    Int_t retval = runl->LoadHeader();
    if (retval){
        cerr<<"AliITSPrintRecPoints.C : LoadHeader returned error"<<endl;
        return kFALSE;
    } // end if retval

    AliITSLoader* ITSloader =  (AliITSLoader*) runl->GetLoader("ITSLoader");
    if(!ITSloader){
        cerr<<"AliITSPrintRecPoints.C :  ITS loader not found"<<endl;
        return kFALSE;
    } // end if !ITSloader

    if(!gGeoManager){
        gGeoManger = new TGeoManager();
        gGeoManager->Import("geometry.root","");
    } // end if !gGeoManager
    if(ITSloader->LoadHits("read")!=0){
        cerr<<"Error Loading Hits"<<endl;
        return kFALSE;
    }// end if ITSloader->LoadHits
    if(ITSloader->LoadRecPoints("read")!=0){
        cerr<<"Error Loading RecPoints"<<endl;
        return kFALSE;
    }// end if ITSloader->LoadRecPoints
    AliITS *ITS = 0;
    ITS = (AliITS*)(gAlice->GetDetector("ITS"));
    if(!ITS){
        cout << "Error: no ITS found. Aborting"<<endl;
        return kFALSE;
    } // end if !ITS
    if(!(ITS->GetDetTypeSim())){
        cout <<"No AliITSDetTypeSim object found in ITS"<<endl;
        return kFALSE;
    } // end if
    AliITSmodule     *m   = 0;
    AliITShit        *h   = 0;
    AliITSDetTypeSim *sim = ITS->GetDetTypeSim();
    AliITSgeom *gm=0;
    gm = ITS->GetITSgeom();
    if(!gm){
        cout <<"No AliITSgeom object found in ITS"<<endl;
        if(!gGeoManager){
            cout <<"No gGeoManger. Aborting"<<endl;
            return kFALSE;
        }else{
            ITS->UpdateInternalGeometry(gGeoManager);
        } // end if
    } // end if !AliITSgeom
    //
    Int_t nMods= gm->GetIndexMax(),nEvents=AliRunLoader::GetNumberOfEvents();
    Int_t mod=0,evnt=0,size=-1,irp=0,ih=0,trkindexOld=-1;
    Double_t xmod,nHitPerTrack;
    TTree            *rpt = 0;
    TClonesArray     *rpa = 0;
    AliITSRecPoint   *rp  = 0;
    AliITSDetTypeRec *rec = new AliITSDetTypeRec(ITSloader);
    rec->SetDefaults();
    // We are going to use the following
    // <hit_x[i]> == { Sum_j hit_x[i][j] }/N_hits per track
    // d_x[i] == <hit_x[i]>-recpoint_x[i]
    // <d_x>  == {sum_i d_x[i]}/N_RecPoints   (excluding noise and no merging)
    //         = {sum_i (<hit_x[i]>-recpoint_x[i])}/N_RecPoints
    //         = {sum_i ((sum_j hit_x[i][j])/N_hits_per_track -recpoinnt_x[i])}/N_Recpoints
    //         = sum_i sum_j {hit_x[i][j]/N_hits_per_track}/N_Recpoints -
    //           sum_i {recpoint_x[i]}/N_Recpoints
    pxg = new TProfile("XdiffGlobal","Mean displacement in gX direction",nMods,
                       -0.5,(Double_t)(nMods)+0.5," ");
    pxg->SetXTitle("module");
    pxg->SetYTitle("Global X [cm]");
    pyg = new TProfile("YdiffGlobal","Mean displacement in gY direction",nMods,
                       -0.5,(Double_t)(nMods)+0.5," ");
    pyg->SetXTitle("module");
    pyg->SetYTitle("Global Y [cm]");
    pzg = new TProfile("ZdiffGlobal","Mean displacement in gZ direction",nMods,
                       -0.5,(Double_t)(nMods)+0.5," ");
    pzg->SetXTitle("module");
    pzg->SetYTitle("Global Z [cm]");
    pxl = new TProfile("XdiffLocal","Mean displacement in lX direction",nMods,
                       -0.5,(Double_t)(nMods)+0.5," ");
    pxl->SetXTitle("module");
    pxl->SetYTitle("local X [cm]");
    pyl = new TProfile("YdiffLocal","Mean displacement in lY direction",nMods,
                       -0.5,(Double_t)(nMods)+0.5," ");
    pyl->SetXTitle("module");
    pyl->SetYTitle("local Y [cm]");
    pzl = new TProfile("ZdiffLocal","Mean displacement in lZ direction",nMods,
                       -0.5,(Double_t)(nMods)+0.5," ");
    pzl->SetXTitle("module");
    pzl->SetYTitle("local Z [cm]");
    //
    for(evnt=0;evnt<nEvents;evnt++){
        runl->GetEvent(evnt);
        ITS->InitModules(size,nMods);
        ITS->FillModules(evnt,0,-1," "," ");
        rec->SetTreeAddress();
        for(mod=0;mod<nMods;mod++){
            xmod = (Double_t) mod;
            m = ITS->GetModule(mod);
            rec->ResetRecPoints();
            rpt = ITSloader->TreeR();
            rpt->GetEvent(mod);
            rpa = rec->RecPoints();
            for(irp=0;irp<rpa->GetEntriesFast();irp++){
                rp = (AliITSRecPoint*)(rpa->At(irp));
                rp->GetGlobalXYZ(rgf);rg[0]=rgf[0];rg[1]=rgf[1];rg[2]=rgf[2];
                rl[0] = rp->GetDetLocalX();
                rl[1] = 0.0;
                rl[2] = rp->GetDetLocalZ();
                pxg->Fill(xmod,-rg[0],0.5);
                pyg->Fill(xmod,-rg[1],0.5);
                pzg->Fill(xmod,-rg[2],0.5);
                pxl->Fill(xmod,-rl[0],0.5);
                //pyl->Fill(xmod,-rl[1],0.5); // assumed to be zero always.
                pzl->Fill(xmod,-rl[2],0.5);
            }// rnf got itp
            trkindexOld=-1;
            nHitPerTrack = 0.0;
            for(ih=0;ih<m->GetNhits();ih++){ // We want the median hit location
                // for each track.
                h = m->GetHit(ih);
                if(m->GetHitTrackIndex(ih)!=trkindexOld){// Enterence location
                    trkindexOld = m->GetHitTrackIndex(ih);
                    nHitPerTrack = 1.0;
                    do{
                        nHitPerTrack += 1.0;
                    }while (m->GetHitTrackIndex(ih+nHitPerTrack-1)==trkindexOld);
                    h->GetPositionG0(hg[0],hg[1],hg[2],tof);
                    h->GetPositionL0(hl[0],hl[1],hl[2],tof);
                    pxg->Fill(xmod,hg[0],0.5/nHitPerTrack);
                    pyg->Fill(xmod,hg[1],0.5/nHitPerTrack);
                    pzg->Fill(xmod,hg[2],0.5/nHitPerTrack);
                    pxl->Fill(xmod,hl[0],0.5/nHitPerTrack);
                    pyl->Fill(xmod,hl[1],1.0/nHitPerTrack); // rl[1]=0 always.
                    pzl->Fill(xmod,hl[2],0.5/nHitPerTrack);
                } // end if
                if(nHitPerTrack==0.0) continue;
                h->GetPositionG(hg[0],hg[1],hg[2],tof);
                h->GetPositionL(hl[0],hl[1],hl[2],tof);
                pxg->Fill(xmod,hg[0],0.5/nHitPerTrack);
                pyg->Fill(xmod,hg[1],0.5/nHitPerTrack);
                pzg->Fill(xmod,hg[2],0.5/nHitPerTrack);
                pxl->Fill(xmod,hl[0],0.5/nHitPerTrack);
                pyl->Fill(xmod,hl[1],1.0/nHitPerTrack); // rl[1]=0 always.
                pzl->Fill(xmod,hl[2],0.5/nHitPerTrack);
            } // end for ih
        } // end for mod
        ITS->ClearModules();
    } // end for evnt
    //
    Int_t wh=800;
    TCanvas *c0 = new TCanvas("c0","Average displacements between hits and RecPoints"
                              ,1.5*wh,wh);
    c0->Divide(3,2);
    c0->cd(1);
    pxg->Draw();
    c0->cd(2);
    pyg->Draw();
    c0->cd(3);
    pzg->Draw();
    c0->cd(4);
    pxl->Draw();
    c0->cd(5);
    pyl->Draw();
    c0->cd(6);
    pzl->Draw();
    c0->Update();
    return kTRUE;
}
Commit	Line	Data
012f0f4c	1	Bool_t AliITSCompareHitsRecPoints(Char_t *rfn="galice.root"){
	2	// This macro compares the average location of the hits associated
	3	// with a specific track and the best match RecPoint location.
	4	// Inputs:
	5	// none.
	6	// Outputs:
	7	// none.
	8	// Return:
	9	// kTRUE if no errors were encountered, otherwise kFALSE.
	10	TProfile pxg,pyg,pzg,pxl,pyl,pzl;
	11	Double_t hg[3],hl[3],tof,rg[3],rl[3];
	12	Float_t rgf[3];
	13
	14	// Dynamically link some shared libs
	15	if (gClassTable->GetID("AliRun") < 0) {
	16	gROOT->LoadMacro("loadlibs.C");
	17	loadlibs();
	18	} else {
	19	if(gAlice){
33c3c91a	20	delete AliRunLoader::Instance();
012f0f4c	21	delete gAlice;
	22	gAlice=0;
	23	} // end if gAlice
	24	}// end if gClassTable->GetID()
	25
	26	gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
	27
	28	AliRunLoader *runl = AccessFile(rfn); // Set up to read in Data
	29	Int_t retval = runl->LoadHeader();
	30	if (retval){
	31	cerr<<"AliITSPrintRecPoints.C : LoadHeader returned error"<<endl;
	32	return kFALSE;
	33	} // end if retval
	34
	35	AliITSLoader* ITSloader = (AliITSLoader*) runl->GetLoader("ITSLoader");
	36	if(!ITSloader){
	37	cerr<<"AliITSPrintRecPoints.C : ITS loader not found"<<endl;
	38	return kFALSE;
	39	} // end if !ITSloader
	40
	41	if(!gGeoManager){
	42	gGeoManger = new TGeoManager();
	43	gGeoManager->Import("geometry.root","");
	44	} // end if !gGeoManager
	45	if(ITSloader->LoadHits("read")!=0){
	46	cerr<<"Error Loading Hits"<<endl;
	47	return kFALSE;
	48	}// end if ITSloader->LoadHits
	49	if(ITSloader->LoadRecPoints("read")!=0){
	50	cerr<<"Error Loading RecPoints"<<endl;
	51	return kFALSE;
	52	}// end if ITSloader->LoadRecPoints
	53	AliITS *ITS = 0;
	54	ITS = (AliITS*)(gAlice->GetDetector("ITS"));
	55	if(!ITS){
	56	cout << "Error: no ITS found. Aborting"<<endl;
	57	return kFALSE;
	58	} // end if !ITS
	59	if(!(ITS->GetDetTypeSim())){
	60	cout <<"No AliITSDetTypeSim object found in ITS"<<endl;
	61	return kFALSE;
	62	} // end if
	63	AliITSmodule *m = 0;
	64	AliITShit *h = 0;
	65	AliITSDetTypeSim *sim = ITS->GetDetTypeSim();
	66	AliITSgeom *gm=0;
	67	gm = ITS->GetITSgeom();
	68	if(!gm){
	69	cout <<"No AliITSgeom object found in ITS"<<endl;
	70	if(!gGeoManager){
	71	cout <<"No gGeoManger. Aborting"<<endl;
	72	return kFALSE;
	73	}else{
	74	ITS->UpdateInternalGeometry(gGeoManager);
	75	} // end if
	76	} // end if !AliITSgeom
	77	//
880d6abe	78	Int_t nMods= gm->GetIndexMax(),nEvents=AliRunLoader::GetNumberOfEvents();
012f0f4c	79	Int_t mod=0,evnt=0,size=-1,irp=0,ih=0,trkindexOld=-1;
	80	Double_t xmod,nHitPerTrack;
	81	TTree *rpt = 0;
	82	TClonesArray *rpa = 0;
	83	AliITSRecPoint *rp = 0;
	84	AliITSDetTypeRec *rec = new AliITSDetTypeRec(ITSloader);
	85	rec->SetDefaults();
	86	// We are going to use the following
	87	// <hit_x[i]> == { Sum_j hit_x[i][j] }/N_hits per track
	88	// d_x[i] == <hit_x[i]>-recpoint_x[i]
	89	// <d_x> == {sum_i d_x[i]}/N_RecPoints (excluding noise and no merging)
	90	// = {sum_i (<hit_x[i]>-recpoint_x[i])}/N_RecPoints
	91	// = {sum_i ((sum_j hit_x[i][j])/N_hits_per_track -recpoinnt_x[i])}/N_Recpoints
	92	// = sum_i sum_j {hit_x[i][j]/N_hits_per_track}/N_Recpoints -
	93	// sum_i {recpoint_x[i]}/N_Recpoints
	94	pxg = new TProfile("XdiffGlobal","Mean displacement in gX direction",nMods,
	95	-0.5,(Double_t)(nMods)+0.5," ");
	96	pxg->SetXTitle("module");
	97	pxg->SetYTitle("Global X [cm]");
	98	pyg = new TProfile("YdiffGlobal","Mean displacement in gY direction",nMods,
	99	-0.5,(Double_t)(nMods)+0.5," ");
	100	pyg->SetXTitle("module");
	101	pyg->SetYTitle("Global Y [cm]");
	102	pzg = new TProfile("ZdiffGlobal","Mean displacement in gZ direction",nMods,
	103	-0.5,(Double_t)(nMods)+0.5," ");
	104	pzg->SetXTitle("module");
	105	pzg->SetYTitle("Global Z [cm]");
	106	pxl = new TProfile("XdiffLocal","Mean displacement in lX direction",nMods,
	107	-0.5,(Double_t)(nMods)+0.5," ");
	108	pxl->SetXTitle("module");
	109	pxl->SetYTitle("local X [cm]");
	110	pyl = new TProfile("YdiffLocal","Mean displacement in lY direction",nMods,
	111	-0.5,(Double_t)(nMods)+0.5," ");
	112	pyl->SetXTitle("module");
	113	pyl->SetYTitle("local Y [cm]");
	114	pzl = new TProfile("ZdiffLocal","Mean displacement in lZ direction",nMods,
	115	-0.5,(Double_t)(nMods)+0.5," ");
	116	pzl->SetXTitle("module");
	117	pzl->SetYTitle("local Z [cm]");
	118	//
	119	for(evnt=0;evnt<nEvents;evnt++){
	120	runl->GetEvent(evnt);
	121	ITS->InitModules(size,nMods);
	122	ITS->FillModules(evnt,0,-1," "," ");
	123	rec->SetTreeAddress();
	124	for(mod=0;mod<nMods;mod++){
	125	xmod = (Double_t) mod;
	126	m = ITS->GetModule(mod);
	127	rec->ResetRecPoints();
	128	rpt = ITSloader->TreeR();
	129	rpt->GetEvent(mod);
	130	rpa = rec->RecPoints();
	131	for(irp=0;irp<rpa->GetEntriesFast();irp++){
	132	rp = (AliITSRecPoint*)(rpa->At(irp));
	133	rp->GetGlobalXYZ(rgf);rg[0]=rgf[0];rg[1]=rgf[1];rg[2]=rgf[2];
	134	rl[0] = rp->GetDetLocalX();
	135	rl[1] = 0.0;
	136	rl[2] = rp->GetDetLocalZ();
	137	pxg->Fill(xmod,-rg[0],0.5);
	138	pyg->Fill(xmod,-rg[1],0.5);
	139	pzg->Fill(xmod,-rg[2],0.5);
	140	pxl->Fill(xmod,-rl[0],0.5);
	141	//pyl->Fill(xmod,-rl[1],0.5); // assumed to be zero always.
	142	pzl->Fill(xmod,-rl[2],0.5);
143	}// rnf got itp
144	trkindexOld=-1;
145	nHitPerTrack = 0.0;
146	for(ih=0;ih<m->GetNhits();ih++){ // We want the median hit location
147	// for each track.
148	h = m->GetHit(ih);
149	if(m->GetHitTrackIndex(ih)!=trkindexOld){// Enterence location
150	trkindexOld = m->GetHitTrackIndex(ih);
151	nHitPerTrack = 1.0;
152	do{
153	nHitPerTrack += 1.0;
154	}while (m->GetHitTrackIndex(ih+nHitPerTrack-1)==trkindexOld);
155	h->GetPositionG0(hg[0],hg[1],hg[2],tof);
156	h->GetPositionL0(hl[0],hl[1],hl[2],tof);
157	pxg->Fill(xmod,hg[0],0.5/nHitPerTrack);
158	pyg->Fill(xmod,hg[1],0.5/nHitPerTrack);
159	pzg->Fill(xmod,hg[2],0.5/nHitPerTrack);
160	pxl->Fill(xmod,hl[0],0.5/nHitPerTrack);
161	pyl->Fill(xmod,hl[1],1.0/nHitPerTrack); // rl[1]=0 always.
162	pzl->Fill(xmod,hl[2],0.5/nHitPerTrack);
163	} // end if
164	if(nHitPerTrack==0.0) continue;
165	h->GetPositionG(hg[0],hg[1],hg[2],tof);
166	h->GetPositionL(hl[0],hl[1],hl[2],tof);
167	pxg->Fill(xmod,hg[0],0.5/nHitPerTrack);
168	pyg->Fill(xmod,hg[1],0.5/nHitPerTrack);
169	pzg->Fill(xmod,hg[2],0.5/nHitPerTrack);
170	pxl->Fill(xmod,hl[0],0.5/nHitPerTrack);
171	pyl->Fill(xmod,hl[1],1.0/nHitPerTrack); // rl[1]=0 always.
172	pzl->Fill(xmod,hl[2],0.5/nHitPerTrack);
173	} // end for ih
174	} // end for mod
175	ITS->ClearModules();
176	} // end for evnt
177	//
178	Int_t wh=800;
179	TCanvas *c0 = new TCanvas("c0","Average displacements between hits and RecPoints"
180	,1.5*wh,wh);
181	c0->Divide(3,2);
182	c0->cd(1);
183	pxg->Draw();
184	c0->cd(2);
185	pyg->Draw();
186	c0->cd(3);
187	pzg->Draw();
188	c0->cd(4);
189	pxl->Draw();
190	c0->cd(5);
191	pyl->Draw();
192	c0->cd(6);
193	pzl->Draw();
194	c0->Update();
195	return kTRUE;
196	}